diff --git a/runS1S2.py b/runS1S2.py
new file mode 100755
index 0000000000000000000000000000000000000000..7504a5691938f10361cd6f60e74a69ec485f3ec6
--- /dev/null
+++ b/runS1S2.py
@@ -0,0 +1,277 @@
+#!/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()
diff --git a/run_with_actors.py b/run_with_actors.py
new file mode 100644
index 0000000000000000000000000000000000000000..74e7245ebb182c2b8d66676c90acab1b8feddc9c
--- /dev/null
+++ b/run_with_actors.py
@@ -0,0 +1,94 @@
+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()