Add example of a single dendrite to show influence of cv policy (#1562)

.github/workflows/basic.yml

@@ -201,3 +201,4 @@ jobs:
             python python/example/single_cell_swc.py python/example/single_cell_detailed.swc
             python python/example/single_cell_detailed.py python/example/single_cell_detailed.swc
             python python/example/single_cell_detailed_recipe.py python/example/single_cell_detailed.swc
+            python python/example/single_cell_cable.py
\ No newline at end of file

.github/workflows/sanitize.yml

@@ -127,3 +127,4 @@ jobs:
             python python/example/single_cell_swc.py python/example/single_cell_detailed.swc
             python python/example/single_cell_detailed.py python/example/single_cell_detailed.swc
             python python/example/single_cell_detailed_recipe.py python/example/single_cell_detailed.swc
+            python python/example/single_cell_cable.py
\ No newline at end of file

doc/tutorial/index.rst

@@ -24,5 +24,6 @@ Tutorials
    single_cell_recipe
    single_cell_detailed
    single_cell_detailed_recipe
+   single_cell_cable
    network_ring
    mpi

doc/tutorial/single_cell_cable.rst

+.. _tutorialsinglecellcable:
+
+A simple dendrite
+==============================
+
+In this example, we will set up a single dendrite represented by a passive cable.
+On one end the dendrite is stimulated by a short current pulse.
+We will investigate how this pulse travels along the dendrite and calculate the conduction velocity.
+
+.. Note::
+
+   **Concepts covered in this example:**
+
+   1. Creating a simulation recipe of a single dendrite.
+   2. Placing probes on the morphology.
+   3. Running the simulation and extracting the results.
+   4. Investigating the influence of control volume policies.
+
+The full code
+*************
+
+You can find the full code of the example at ``python/examples/single_cell_cable.py``
+
+Executing the script will run the simulation with default parameters.
+
+
+Walk-through
+**********
+
+We set up a recipe for the simulation of a single dendrite with some parameters.
+
+.. literalinclude:: ../../python/example/single_cell_cable.py
+   :language: python
+   :lines: 11-31
+
+Implementing the ``cell_description`` member function constructs the morphology and sets the properties of the dendrite as well as the current stimulus and the discretization policy.
+
+.. literalinclude:: ../../python/example/single_cell_cable.py
+   :language: python
+   :lines: 63-98
+
+We parse the command line arguments, instantiate the recipe, run the simulation, extract results and plot:
+
+.. literalinclude:: ../../python/example/single_cell_cable.py
+   :language: python
+   :lines: 147-210
+
+The output plot below shows how the current pulse is attenuated and
+gets broader the further it travels along the dendrite from the current
+clamp.
+
+.. figure:: single_cell_cable_result.svg
+    :width: 800
+    :align: center
+
+The conduction velocity in simulation is calculated from the time of the maximum of the membrane potential.
+
+.. literalinclude:: ../../python/example/single_cell_cable.py
+   :language: python
+   :lines: 212-
+
+Please have in mind that the calculated (idealized) conduction velocity is only correct for an infinite cable.
+
+::
+
+	calculated conduction velocity: 0.47 m/s
+	simulated conduction velocity:  0.50 m/s
+
+When we set the control volume policy to cover the full dendrite, all probes will see the same voltage.
+
+.. code-block:: bash
+
+	./single_cell_cable.py --length 1000 --cv_policy_max_extent 1000
+
+.. figure:: single_cell_cable_result_cv_policy.svg
+    :width: 800
+    :align: center

python/example/single_cell_cable.py

+#!/usr/bin/env python3
+
+import arbor
+import argparse
+import numpy as np
+
+import pandas
+import seaborn  # You may have to pip install these.
+
+
+class Cable(arbor.recipe):
+    def __init__(self, probes,
+                 Vm, length, radius, cm, rL, g,
+                 stimulus_start, stimulus_duration, stimulus_amplitude,
+                 cv_policy_max_extent):
+        """
+        probes -- list of probes
+
+        Vm -- membrane leak potential
+        length -- length of cable in μm
+        radius -- radius of cable in μm
+        cm -- membrane capacitance in F/m^2
+        rL -- axial resistivity in Ω·cm
+        g -- membrane conductivity in S/cm^2
+
+        stimulus_start -- start of stimulus in ms
+        stimulus_duration -- duration of stimulus in ms
+        stimulus_amplitude -- amplitude of stimulus in nA
+
+        cv_policy_max_extent -- maximum extent of control volume in μm
+        """
+
+        # (4.1) The base C++ class constructor must be called first, to ensure that
+        # all memory in the C++ class is initialized correctly.
+        arbor.recipe.__init__(self)
+
+        self.the_probes = probes
+
+        self.Vm = Vm
+        self.length = length
+        self.radius = radius
+        self.cm = cm
+        self.rL = rL
+        self.g = g
+
+        self.cv_policy_max_extent = cv_policy_max_extent
+
+        self.the_props = arbor.neuron_cable_properties()
+        self.the_cat = arbor.default_catalogue()
+        self.the_props.register(self.the_cat)
+
+    def num_cells(self):
+        return 1
+
+    def num_sources(self, gid):
+        assert gid == 0
+        return 0
+
+    def cell_kind(self, gid):
+        assert gid == 0
+        return arbor.cell_kind.cable
+
+    def cell_description(self, gid):
+        """A high level description of the cell with global identifier gid.
+
+        For example the morphology, synapses and ion channels required
+        to build a multi-compartment neuron.
+        """
+        assert gid == 0
+
+        tree = arbor.segment_tree()
+
+        tree.append(arbor.mnpos,
+                    arbor.mpoint(0, 0, 0, self.radius),
+                    arbor.mpoint(self.length, 0, 0, self.radius),
+                    tag=1)
+
+        labels = arbor.label_dict({'cable': '(tag 1)',
+                                   'start': '(location 0 0)'})
+
+        decor = arbor.decor()
+        decor.set_property(Vm=self.Vm)
+        decor.set_property(cm=self.cm)
+        decor.set_property(rL=self.rL)
+
+        decor.paint('"cable"',
+                    arbor.mechanism(f'pas/e={self.Vm}', {'g': self.g}))
+
+        decor.place('"start"', arbor.iclamp(args.stimulus_start,
+                                            args.stimulus_duration,
+                                            args.stimulus_amplitude))
+
+        policy = arbor.cv_policy_max_extent(self.cv_policy_max_extent)
+        decor.discretization(policy)
+
+        cell = arbor.cable_cell(tree, labels, decor)
+
+        return cell
+
+    def probes(self, gid):
+        assert gid == 0
+        return self.the_probes
+
+    def global_properties(self, kind):
+        return self.the_props
+
+
+def get_rm(g):
+    """Return membrane resistivity in Ohm*m^2
+    g -- membrane conductivity in S/m^2
+    """
+    return 1/g
+
+
+def get_taum(cm, rm):
+    """Return membrane time constant in seconds
+    cm -- membrane capacitance in F/m^2
+    rm -- membrane resistivity in Ohm*m^2
+    """
+    return cm*rm
+
+
+def get_lambdam(a, rm, rL):
+    """Return electronic length in m
+    a -- cable radius in m
+    rm -- membrane resistivity in Ohm*m^2
+    rL -- axial resistivity in Ohm*m
+    """
+    return np.sqrt(a*rm/(2*rL))
+
+
+def get_vcond(lambdam, taum):
+    """Return conductance velocity in m/s
+    lambda -- electronic length in m
+    taum -- membrane time constant
+    """
+    return 2*lambdam/taum
+
+
+def get_tmax(data):
+    """Return time of maximum of membrane voltage"""
+    return data[np.argmax(data[:, 1])][0]
+
+
+if __name__ == "__main__":
+
+    parser = argparse.ArgumentParser(description='Cable')
+
+    parser.add_argument(
+        '--Vm', help="membrane leak potential in mV", type=float, default=-65)
+    parser.add_argument(
+        '--length', help="cable length in μm", type=float, default=1000)
+    parser.add_argument(
+        '--radius', help="cable radius in μm", type=float, default=1)
+    parser.add_argument(
+        '--cm', help="membrane capacitance in F/m^2", type=float, default=0.01)
+    parser.add_argument(
+        '--rL', help="axial resistivity in Ω·cm", type=float, default=90)
+    parser.add_argument(
+        '--g', help="membrane conductivity in S/cm^2", type=float, default=0.001)
+
+    parser.add_argument('--stimulus_start',
+                        help="start of stimulus in ms", type=float, default=10)
+    parser.add_argument('--stimulus_duration',
+                        help="duration of stimulus in ms", type=float, default=0.1)
+    parser.add_argument('--stimulus_amplitude',
+                        help="amplitude of stimulus in nA", type=float, default=1)
+
+    parser.add_argument('--cv_policy_max_extent',
+                        help="maximum extent of control volume in μm", type=float,
+                        default=10)
+
+    # parse the command line arguments
+    args = parser.parse_args()
+
+    # set up membrane voltage probes equidistantly along the dendrites
+    probes = [arbor.cable_probe_membrane_voltage(
+        f'(location 0 {r})') for r in np.linspace(0, 1, 11)]
+    recipe = Cable(probes, **vars(args))
+
+    # create a default execution context and a default domain decomposition
+    context = arbor.context()
+    domains = arbor.partition_load_balance(recipe, context)
+
+    # configure the simulation and handles for the probes
+    sim = arbor.simulation(recipe, domains, context)
+    dt = 0.001
+    handles = [sim.sample((0, i), arbor.regular_schedule(dt))
+               for i in range(len(probes))]
+
+    # run the simulation for 30 ms
+    sim.run(tfinal=30, dt=dt)
+
+    # retrieve the sampled membrane voltages and convert to a pandas DataFrame
+    data = [sim.samples(handle)[0][0] for handle in handles]
+    data_dict = {str(i): d[:, 1] for i, d in enumerate(data)}
+    data_dict["t"] = data[0][:, 0]
+
+    df = pandas.DataFrame.from_dict(data_dict)
+
+    # plot all probes and save to file
+    for i in range(len(probes)):
+        plot = seaborn.lineplot(data=df, x="t", y=str(
+            i), ci=None, legend="full", label=i)
+
+    plot.set(xlabel="t/ms")
+    plot.set(ylabel="U/mV")
+    plot.set_xlim(9, 14)
+
+    plot.figure.savefig("single_cell_cable_result.svg")
+
+    # calculcate the idealized conduction velocity, cf. cable equation
+    rm = get_rm(args.g*1/(0.01*0.01))
+    taum = get_taum(args.cm, rm)
+    lambdam = get_lambdam(args.radius*1e-6, rm, args.rL*0.01)
+    vcond_ideal = get_vcond(lambdam, taum)
+
+    # take the last and second probe
+    delta_t = (get_tmax(data[-1]) - get_tmax(data[1]))
+
+    # 90% because we took the second probe
+    delta_x = args.length*0.9
+    vcond = delta_x*1e-6/(delta_t*1e-3)
+
+    print(f"calculated conduction velocity: {vcond_ideal:.2f} m/s")
+    print(f"simulated conduction velocity:  {vcond:.2f} m/s")