add files for running a simulation within an actor

runS1S2.py

+#!/usr/bin/env python
+
+"""
+Overview
+========
+
+.. code-block:: bash
+
+    cd ${TUTORIALS}/02_EP_tissue/01_basic_usage
+
+This example introduces to the basics of using the openCARP executable 
+for simulating EP at the tissue and organ scale.
+Details on how to inquire input parameters are provided in the :ref:`openCARP <openCARP>` section
+of the manual.
+
+In absence of pacemaker currents electrical stimulus currents must be supplied in order to 
+depolarize tissue above the threshold to elicit a propagated action potential. 
+The simplest way of initiating a propagating action potential is the use of a transmembrane stimulus, :math:`I_{\mathrm{tr}}`.
+According to the definition in openCARP supplying a positive :math:`I_{\mathrm{tr}}` depolarizes the tissue,
+a negative :math:`I_{\mathrm{tr}}` hyperpolarizes the tissue.
+
+
+Experimental setup
+==================
+
+To introduce definition and usage of transmembrane stimulus currents in openCARP a simple experiment 
+using a thin monolayer preparation was defined. The following input parameters are available to steer the experiment:
+
+.. code-block:: bash
+
+  --duration    
+                Duration of simulation (ms)
+
+  --S1-strength 
+                pick transmembrane current stimulus strength in
+                [uA/cm^2] = [pA/pF] considering fixed Cm (default: 20.)
+
+  --S1-dur     
+                pick transmembrane current stimulus duration in [ms
+
+
+.. _tutorial-subthreshold-stimuation:
+
+Modeling subthreshold conditions
+--------------------------------
+
+.. _exp01-ep-ti-basic:
+
+**Experiment exp01 (subthreshold depolarizing stimulus)**
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+We start with a longer stimulus lasting for 15 ms and gradually increase the stimulus strength
+until we achieve capture. This is equivalent to measuring one data point of a `strength-duration curve <https://en.wikipedia.org/wiki/Rheobase>`_
+to determine chronaxie and rheobase of our preparation.
+
+.. code-block:: bash
+
+   ./run.py --duration 20 --S1-strength 20. --S1-dur 15 --visualize
+
+This weak current slowly depolarizes the membrane, but does not reach threshold.
+After the end of stimulation, the membrane returns to its resting state, as shown in :numref:`fig-tutorial-subthreshold-depol-stim`.
+
+.. _fig-tutorial-subthreshold-depol-stim:
+
+.. figure:: /images/02_01_subthreshold_stim_response.png 
+   :width: 50%
+   :align: center
+
+   Membrane response to a subthresold *depolarizing* transmembrane stimulus.
+
+Inspect the variation in membrane depolarization as a function of distance to the stimulus location 
+by varying the vertex index in meshalyzer.
+
+To understand how the stimulus definition is passed on to the simulator we inspect the command line generated by the script.
+For this sake run
+
+.. code-block:: bash
+
+   ./run.py --duration 20 --S1-strength 20. --S1-dur 15 --dry
+
+The stimulus section of the command line is then
+
+.. code-block:: bash
+
+   -num_stim 1 \                   # use one stimulus electrode
+   -stimulus[0].name S1 \
+    \                              # pulse and type definition
+   -stimulus[0].stimtype 0 \
+   -stimulus[0].strength 20.0 \
+   -stimulus[0].duration 15.0 \
+    \                              # electrode definition (block)       
+   -stimulus[0].x0 -5050.0 \
+   -stimulus[0].xd 100.0 \
+   -stimulus[0].y0 -550.0 \
+   -stimulus[0].yd 1100.0 \
+   -stimulus[0].z0 -150.0 \
+   -stimulus[0].zd 300.0
+
+
+
+**Experiment exp02 (subthreshold hyperpolarizing stimulus)**
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+We repeat exp01 with a negative transmembrane stimulus current.
+
+.. code-block:: bash
+
+   ./run.py --duration 20 --S1-strength -20 --S1-dur 15 --visualize
+
+
+As the membrane behaves essentially only passively in the subthreshold regime 
+we observe the same response symmetric to the resting potential, see :numref:`fig-tutorial-subthreshold-hyper-stim`
+
+.. _fig-tutorial-subthreshold-hyper-stim:
+
+.. figure:: /images/02_01_subthreshold_stim_response_hyper.png 
+   :width: 50%
+   :align: center
+
+   Membrane response to a subthresold *hyperpolarizing* transmembrane stimulus.
+
+
+.. _tutorial-suprathreshold-stimuation:
+
+Initiating action potential propagation
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+We repeat exp01 with a stronger transmembrane stimulus current of shorter duration.
+Peak sodium current is known to be around :math:`250 \mu A/cm^2` and last for less than a millisecond.
+We chose a strength of  :math:`250 \mu A/cm^2` of 1 millisecond duration to achieve suprathreshold depolarization.
+Keep in mind that *suprathreshold depolarization* at the stimulus site does not automatically lead to a propagated response.
+The local membrane response at the stimulus location is shown in :numref:`fig-tutorial-suprathreshold-stim`.
+
+.. code-block:: bash
+
+   ./run.py --duration 20 --S1-strength 250. --S1-dur 2 --visualize
+
+
+.. _fig-tutorial-suprathreshold-stim:
+
+.. figure:: /images/02_01_suprathreshold_stim_response.png 
+   :width: 50%
+   :align: center
+
+   Membrane response to a suprathresold *depolarizing* transmembrane stimulus.
+
+"""
+import os
+EXAMPLE_DESCRIPTIVE_NAME = 'Basics of usage '
+EXAMPLE_AUTHOR = 'Gernot Plank <gernot.plank@medunigraz.at>'
+EXAMPLE_DIR = os.path.dirname(__file__)
+GUIinclude = True
+
+from datetime import date
+from carputils import settings
+from carputils import tools
+from carputils import mesh
+from carputils.carpio import txt
+import numpy as np
+from numpy import array as nplist
+
+#Joyce: Can add additional options for arguments here
+def parser():
+    parser = tools.standard_parser()
+    group  = parser.add_argument_group('experiment specific options')
+    group.add_argument('--duration',
+                        type = float,
+                        default = 300.,
+                        help = 'Duration of simulation in [ms] (default: 20.)')
+    group.add_argument('--S1-strength',
+                       type = float,
+                       default = 172.,
+                       help = 'pick S1 transmembrane current stimulus strength in [uA/cm^2] = [pA/pF] considering fixed Cm (default: 20.)')
+    group.add_argument('--S1-dur',
+                        type = float,
+                        default = 2.,
+                        help = 'pick transmembrane current stimulus duration in [ms] (default: 2.)')
+    group.add_argument('--S1-start',
+                        type = float,
+                        default = 0.0,
+                        help = 'pick S1 start time')
+    group.add_argument('--S2-strength',
+                       type = float,
+                       default = 200.,
+                       help = 'pick S2 transmembrane current stimulus strength in [uA/cm^2] = [pA/pF] considering fixed Cm (default: 20.)')
+    group.add_argument('--S2-dur',
+                        type = float,
+                        default = 2.,
+                        help = 'pick transmembrane current stimulus duration in [ms] (default: 2.)')
+    group.add_argument('--S2-start',
+                        type = float,
+                        default = 150.0,
+                        help = 'pick S2 start time')
+    return parser
+
+def jobID(args):
+    """
+    Generate name of top level output directory.
+    """
+    today = date.today()
+    return '{}_basic_{}'.format(today.isoformat(), args.S1_strength)
+
+@tools.carpexample(parser, jobID)
+def run(args, job):
+
+    #===========================================
+    # 1: Define the mesh: Units are mm
+    #===========================================
+    
+    # Create a block with a thin z dimension
+    x = 10.0 #(mm)
+    y = 1.0  #(mm)
+    z = 0.2  #(mm)
+    geom = mesh.Block(size=(x,y,z)) 
+
+    # Set fibre angle to 0, sheet angle to 0
+    geom.set_fibres(0, 0, 90, 90) #endo, epi, sheet endo, sheet epicard
+
+    # Finilaize the creation of the mesh
+    meshname = mesh.generate(geom)
+
+
+    #===========================================
+    # 2: Define the ionic models & conductivities
+    #===========================================
+
+    # Query for element labels
+    _, etags, _ = txt.read(meshname + '.elem')
+    etags = np.unique(etags)
+
+    IntraTags = etags[etags != 0].tolist()
+    ExtraTags = etags.tolist().copy()
+
+    #===========================================
+    # 3: Configure the simulation
+    #===========================================
+
+    # Define the geometry of the stimulus at one end of the block
+    # Units are um
+    stim = ['-num_stim',             2,
+            '-stimulus[0].name',    'S1', #index 0 because only one type of electrode
+            '-stimulus[0].stimtype', 0, #0 means intracellular stim
+            '-stimulus[0].strength', args.S1_strength, #uA/Cm^2
+            '-stimulus[0].duration', args.S1_dur,
+            '-stimulus[0].start',    args.S1_start,
+            '-stimulus[1].name',    'S2',
+            '-stimulus[1].stimtype', 0,
+            '-stimulus[1].strength', args.S2_strength,
+            '-stimulus[1].duration', args.S2_dur,
+            '-stimulus[1].start',    args.S2_start]
+
+    #electrode = mesh.block_boundary_condition(geom, 'stimulus', 0, 'x', True) #True refers to a BC between tissue and bath
+    E1_lower_bound = nplist([ 0.0, 0.0, 0.0])
+    E1_upper_bound = nplist([ 2.0, 0.3, 0.2])
+    electrode = mesh.block_region(geom, 'stimulus', 0, E1_lower_bound, E1_upper_bound, False) # introudce a stimulus in the region defined by the lower and upper bounds
+
+    #===========================================
+    # 4: Defining simulator options:
+    #===========================================
+
+    # Get basic command line, including solver options
+    cmd = tools.carp_cmd(os.path.join(EXAMPLE_DIR, 'basic.par')) #cmd is a carp command
+
+    cmd += ['-simID',    job.ID,
+            '-meshname', meshname,
+            '-dt',       25,
+            '-tend',     args.duration]
+    cmd += tools.gen_physics_opts(ExtraTags=ExtraTags, IntraTags=IntraTags)
+    cmd += stim + electrode
+
+    # Run simulation 
+    job.carp(cmd) #cmd includes the params
+
+
+
+if __name__ == '__main__':
+    run()    

run_with_actors.py

+from thespian.actors import *
+import runS1S2
+import datetime
+import os
+import re
+from subprocess import Popen, PIPE
+import subprocess
+import numpy as np
+
+
+class output_actor(Actor):
+    def receiveMessage(self, msg, sender):
+        S2_interval = msg[1]
+        # subprocess.run(msg)
+        file = msg[0][-1]
+
+        lines = np.loadtxt(file)
+        desired_interval = lines[S2_interval:]
+        self.send(sender, bool(np.any(desired_interval > 0)))
+
+class Simulation_Manager(Actor):
+    def __init__(self):
+        # ===========================================
+        # Define the initial parameters for the S1S2 protocol
+        # ===========================================
+        self.S1_strength = 155.0 # Usually the S1 threshold (for this problem, we will pretend 155 is the S1 threshold)
+        self.t_final = 300 # duration of the experiment
+        self.S2_min = 150000  # Should be set to at least the S1 threshold. (150,000 is closer to the currently set S1-S2 intervals)
+        self.S2_absolute_max = 151000 # If unsure, set to maximum possible strength before solution blows up (for this problem, =225,000)
+        self.S2_max = self.S2_absolute_max # We want to retain the value of S2_absolute_max, so copy to S2_max to start with
+        self.tol = 1 # Find threshold up to a certain tolerance (typically 1 for these experiments).
+        self.S1S2_interval_min = 200 # User sets the range of desired S1-S2 time intervals
+        self.S1S2_interval_max = 220
+        self.S1S2_interval_increment = 5
+        self.S1S2_interval = self.S1S2_interval_max # S1S2 interval is set to max to start with
+        
+        self.nodes = 0 # Number of nodes in the mesh
+        self.actor_responses = []
+
+    def receiveMessage(self, message, sender):
+        if isinstance(message, str):
+            args = ["--S1-strength", str(self.S1_strength), "--S2-strength", str(self.S2_min - 1), \
+                    "--S2-start", str(self.S1S2_interval), "--duration", str(self.t_final)]
+            # runS1S2.run(args)
+
+            # check if output directory exists
+
+            output_folder = str(datetime.date.today().isoformat()) + "_basic_" + str(self.S1_strength)
+            if not os.path.exists(output_folder):
+                Exception("Output folder does not exist")
+            else:
+                # Find out how many nodes are in the mesh
+                # TODO: This is a hacky way to do this. Should be a better way to do this.
+                output_file = output_folder + "/vm.igb"
+                process = Popen(["/opt/openCARP/_build/bin/igbhead", output_file], stdout=PIPE)
+                (output, err) = process.communicate()
+                exit_code = process.wait()
+                string_output = output.decode("utf-8").split("\n")
+                nodes = re.findall(r"\d+", string_output[0])
+                nodes = int(nodes[0])
+                self.nodes = nodes
+
+            # Inspect the mesh
+            
+            # Make a temporary directory to store the output values
+            temp_dir = output_folder + "/temp"
+            if not os.path.exists(temp_dir):
+                os.makedirs(temp_dir)
+            else:
+                Exception("Temporary directory already exists")
+
+            
+            for i in range(0, self.nodes):
+                data_file = temp_dir + "/out_" + str(i) + ".txt"
+                cmd = (["/opt/openCARP/_build/bin/igbextract", "-o", "ascii", "-l", str(i), output_file, "-O", data_file], self.S1S2_interval)
+                actor = self.createActor(output_actor)
+                self.send(actor, cmd)
+        
+        elif isinstance(message, bool):
+            self.actor_responses.append(message)
+            if len(self.actor_responses) == self.nodes:
+                if not all(self.actor_responses):
+                    print("APs not generated")
+
+
+
+
+def actor_system_start():
+  print("Starting actor system")
+  first = ActorSystem().createActor(Simulation_Manager)
+  ActorSystem().tell(first, "start")
+
+if __name__ == '__main__':
+  actor_system_start()