Example and tutorial of two cells connected via gap-junction (#1771)

doc/tutorial/index.rst

@@ -27,3 +27,4 @@ Tutorials
    single_cell_cable
    network_ring
    network_ring_mpi
+   two_cells_gap_junctions

doc/tutorial/two_cells_gap_junctions.rst

+.. _tutorialtwocellsgapjunction:
+
+Two cells connected via a gap junction
+======================================
+
+In this example, we will set up two cells connected via a gap junction.
+The cells have different leak potentials.
+We will investigate how the equilibrium potentials of the two cells change because of the gap junction connection.
+
+.. Note::
+
+   **Concepts covered in this example:**
+
+   1. Creating a simulation recipe for two cells.
+   2. Placing probes.
+   3. Running the simulation and extracting the results.
+   4. Adding a gap junction connection.
+
+The full code
+*************
+
+You can find the full code of the example at ``python/examples/two_cell_gap_junctions.py``
+
+Executing the script will run the simulation with default parameters.
+
+
+Walk-through
+************
+
+We set up a recipe for the simulation of two cells with some parameters.
+
+.. literalinclude:: ../../python/example/two_cell_gap_junctions.py
+   :language: python
+   :lines: -66
+
+Implementing the ``cell_description`` member function constructs the morphology and sets the properties of the cells as well as the gap junction mechanisms and the discretization policy.
+
+.. literalinclude:: ../../python/example/two_cell_gap_junctions.py
+   :language: python
+   :lines: 61-101
+
+The bidirection gap junction is created in the function ``gap_junctions_on``.
+
+.. literalinclude:: ../../python/example/two_cell_gap_junctions.py
+   :language: python
+   :lines: 103-107
+
+We parse the command line arguments, instantiate the recipe, run the simulation, extract results and plot:
+
+.. literalinclude:: ../../python/example/two_cell_gap_junctions.py
+   :language: python
+   :lines: 110-
+
+The output plot below shows how the potential of the two cells approaches their equilibrium potentials.
+The expected values are denoted by dashed lines.
+
+.. figure:: two_cell_gap_junctions_result.svg
+    :width: 800
+    :align: center
+
+The equivalent circuit used to calculate the equilibrium potentials is depicted below.
+
+.. figure:: gap_junction_circuit.svg
+    :width: 800
+    :align: center

python/example/two_cell_gap_junctions.py

+#!/usr/bin/env python3
+
+import arbor
+import argparse
+import numpy as np
+
+import pandas
+import seaborn  # You may have to pip install these.
+import matplotlib.pyplot as plt
+
+
+class TwoCellsWithGapJunction(arbor.recipe):
+    def __init__(self, probes,
+                 Vms, length, radius, cm, rL, g, gj_g,
+                 cv_policy_max_extent):
+        """
+        probes -- list of probes
+
+        Vms -- membrane leak potentials of the two cells
+        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
+        gj_g -- gap junction conductivity in μS
+
+        cv_policy_max_extent -- maximum extent of control volume in μm
+        """
+
+        # 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.Vms = Vms
+        self.length = length
+        self.radius = radius
+        self.cm = cm
+        self.rL = rL
+        self.g = g
+        self.gj_g = gj_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 2
+
+    def num_sources(self, gid):
+        assert gid in [0, 1]
+        return 0
+
+    def cell_kind(self, gid):
+        assert gid in [0, 1]
+        return arbor.cell_kind.cable
+
+    def probes(self, gid):
+        assert gid in [0, 1]
+        return self.the_probes
+
+    def global_properties(self, kind):
+        return self.the_props
+
+    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 in [0, 1]
+
+        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({'cell' : '(tag 1)', 'gj_site': '(location 0 0.5)'})
+
+        decor = arbor.decor()
+        decor.set_property(Vm=self.Vms[gid])
+        decor.set_property(cm=self.cm)
+        decor.set_property(rL=self.rL)
+
+        # add a gap junction mechanism at the "gj_site" location and label that specific mechanism on that location "gj_label"
+        junction_mech = arbor.junction('gj', {"g" : self.gj_g})
+        decor.place('"gj_site"', junction_mech, 'gj_label')
+        decor.paint('"cell"', arbor.density(f'pas/e={self.Vms[gid]}', {'g': self.g}))
+
+        if self.cv_policy_max_extent is not None:
+            policy = arbor.cv_policy_max_extent(self.cv_policy_max_extent)
+            decor.discretization(policy)
+        else:
+            decor.discretization(arbor.cv_policy_single())
+
+        return arbor.cable_cell(tree, labels, decor)
+
+    def gap_junctions_on(self, gid):
+        assert gid in [0, 1]
+
+        # create a bidirectional gap junction from cell 0 at label "gj_label" to cell 1 at label "gj_label" and back.
+        return [arbor.gap_junction_connection((1 if gid == 0 else 0, 'gj_label'), 'gj_label', 1)]
+
+
+if __name__ == "__main__":
+
+    parser = argparse.ArgumentParser(description='Two cells connected via a gap junction')
+
+    parser.add_argument(
+        '--Vms', help="membrane leak potentials in mV", type=float, default=[-100, -60], nargs=2)
+    parser.add_argument(
+        '--length', help="cell length in μm", type=float, default=100)
+    parser.add_argument(
+        '--radius', help="cell radius in μm", type=float, default=3)
+    parser.add_argument(
+        '--cm', help="membrane capacitance in F/m^2", type=float, default=0.005)
+    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(
+        '--gj_g', help="gap junction conductivity in μS", type=float, default=0.01)
+
+    parser.add_argument('--cv_policy_max_extent',
+                        help="maximum extent of control volume in μm", type=float)
+
+    # parse the command line arguments
+    args = parser.parse_args()
+
+    # set up membrane voltage probes at the position of the gap junction
+    probes = [arbor.cable_probe_membrane_voltage('"gj_site"')]
+    recipe = TwoCellsWithGapJunction(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.01
+    handles = []
+    for gid in [0, 1]:
+        handles += [sim.sample((gid, i), arbor.regular_schedule(dt)) for i in range(len(probes))]
+
+    # run the simulation for 5 ms
+    sim.run(tfinal=5, dt=dt)
+
+    # retrieve the sampled membrane voltages and convert to a pandas DataFrame
+    print("Plotting results ...")
+    df_list = []
+    for probe in range(len(handles)):
+        samples, meta = sim.samples(handles[probe])[0]
+        df_list.append(pandas.DataFrame({'t/ms': samples[:, 0], 'U/mV': samples[:, 1], 'Cell': f"{probe}"}))
+
+    df = pandas.concat(df_list,ignore_index=True)
+
+    fig, ax = plt.subplots()
+
+    # plot the membrane potentials of the two cells as function of time
+    seaborn.lineplot(ax=ax,data=df, x="t/ms", y="U/mV",hue="Cell",ci=None)
+
+    # area of cells
+    area = args.length*1e-6 * 2*np.pi*args.radius*1e-6
+
+    # total conductance and resistance
+    cell_g = args.g/1e-4 * area
+    cell_R = 1/cell_g
+
+    # gap junction conductance and resistance in base units
+    si_gj_g = args.gj_g*1e-6
+    si_gj_R = 1/si_gj_g
+
+    # indicate the expected equilibrium potentials
+    for (i, j) in [[0,1], [1,0]]:
+        weighted_potential = args.Vms[i] + ((args.Vms[j] - args.Vms[i])*(si_gj_R + cell_R))/(2*cell_R + si_gj_R)
+        ax.axhline(weighted_potential, linestyle='dashed', color='black', alpha=0.5)
+
+    # plot the initial/nominal resting potentials
+    for gid, Vm in enumerate(args.Vms):
+        ax.axhline(Vm, linestyle='dashed', color='black', alpha=0.5)
+
+    fig.savefig('two_cell_gap_junctions_result.svg')

scripts/run_python_examples.sh

@@ -22,3 +22,4 @@ $PREFIX python python/example/single_cell_detailed.py python/example/single_cell
 $PREFIX python python/example/single_cell_detailed_recipe.py python/example/single_cell_detailed.swc
 $PREFIX python python/example/single_cell_extracellular_potentials.py python/example/single_cell_detailed.swc
 $PREFIX python python/example/single_cell_cable.py
+$PREFIX python python/example/two_cell_gap_junctions.py