Some more polish.

doc/tutorial/network_two_cells_gap_junctions.rst

@@ -29,54 +29,51 @@ We set up a recipe for the simulation of two cells
 
 .. literalinclude:: ../../python/example/network_two_cells_gap_junctions.py
    :language: python
-   :lines: 12-46
+   :lines: 13-45
 
 in which we store the relevant parameters for the two cells, all of which are
 shared except the equilibrium potentials in ``Vms``. Next, we define callbacks
 to define our cell population:
+
 - ``num_cells`` returns the number of cells in the network, ie 2
 - ``cell_kind`` specifies that we handle ``cable_cell`` exclusively.
-- ``global_properties`` returns a list of standard parameters based on the
-  defaults of the NEURON simulator.
-- ``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. It returns the finished ``cable_cell``, matching the
-  ``cell_kind`` callback.
+- ``global_properties`` returns a list of standard parameters based on the defaults of the NEURON simulator.
+- ``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. It returns the finished ``cable_cell``, matching the ``cell_kind`` callback.
 
 .. literalinclude:: ../../python/example/network_two_cells_gap_junctions.py
    :language: python
-   :lines: 48-94
+   :lines: 49-93
 
 We build the network conections, here a single, bidirectional gap junction
 
 .. literalinclude:: ../../python/example/network_two_cells_gap_junctions.py
    :language: python
-   :lines: 96-105
+   :lines: 97-105
 
 And, finally, we return a set of probes which are passed in during construction
 
 .. literalinclude:: ../../python/example/network_two_cells_gap_junctions.py
    :language: python
-   :lines: 107-109
+   :lines: 108-109
 
 We parse the command line arguments which are used to set parameters in the recipe
 
 .. literalinclude:: ../../python/example/network_two_cells_gap_junctions.py
    :language: python
-   :lines: 112-146
+   :lines: 113-147
 
-Next, we define a list of probes, construct the recipe and simulation, setting probe
-sampling on a regular grid with width equal to the timestep :math:`dt`
+Next, we define a list of probes and construct the recipe and simulation
 
 .. literalinclude:: ../../python/example/network_two_cells_gap_junctions.py
    :language: python
-   :lines: 148-160
+   :lines: 149-154
 
-Having set up the simulation, we can run it and access the sampling values
+Having set up the simulation, we setting probe sampling on a regular grid with
+width equal to the timestep :math:`dt`. Now, we can run it and access the sampling values
 
 .. literalinclude:: ../../python/example/network_two_cells_gap_junctions.py
    :language: python
-   :lines: 162-176
+   :lines: 156-178
 
 All that is left to do is to put this into a plot. The output plot below shows
 how the potential of the two cells approaches their equilibrium potentials

python/example/network_two_cells_gap_junctions.py

@@ -1,5 +1,6 @@
 #!/usr/bin/env python3
 
+from builtins import enumerate
 import arbor
 import argparse
 import numpy as np
@@ -95,7 +96,6 @@ def cell_description(self, gid):
 
     def gap_junctions_on(self, gid):
         # create a bidirectional gap junction from cell 0 at label "gj_label" to cell 1 at label "gj_label" and back.
-
         if gid == 0:
             tgt = 1
         elif gid == 1:
@@ -108,6 +108,7 @@ def probes(self, gid):
         assert gid in [0, 1]
         return self.the_probes
 
+
 if __name__ == "__main__":
     parser = argparse.ArgumentParser(
         description="Two cells connected via a gap junction"
@@ -152,28 +153,28 @@ def probes(self, gid):
     # configure the simulation and handles for the probes
     sim = arbor.simulation(recipe)
 
+    T = 5
     dt = 0.01
-    handles = []
-    for gid in [0, 1]:
-        handles += [
-            sim.sample((gid, i), arbor.regular_schedule(dt)) for i in range(len(probes))
-        ]
+
+    handles = [
+        sim.sample((gid, i), arbor.regular_schedule(dt))
+        for i, _ in enumerate(probes)
+        for gid in range(recipe.num_cells())
+    ]
 
     # run the simulation for 5 ms
-    T = 5
     sim.run(tfinal=T, 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]
+    for probe, handle in enumerate(handles):
+        samples, meta = sim.samples(handle)[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()
@@ -184,20 +185,28 @@ def probes(self, gid):
     # area of cells
     area = args.length * 1e-6 * 2 * np.pi * args.radius * 1e-6
 
-    # total conductance
-    cell_g = args.g / 1e-4 * area
-
-    # gap junction conductance in base units
+    # total and gap junction conductance in base units
+    cell_g = area * args.g / 1e-4
     si_gj_g = args.gj_g * 1e-6
 
+    # weight
+    w = (si_gj_g + cell_g) / (2 * si_gj_g + cell_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_g + cell_g)
-        ) / (2 * si_gj_g + cell_g)
+        weighted_potential = args.Vms[i] + w * (args.Vms[j] - args.Vms[i])
         ax.axhline(weighted_potential, linestyle="dashed", color="black", alpha=0.5)
-        ax.text(2, weighted_potential, f'$\\tilde U_{j} = U_{j} + w\\cdot(U_{j} - U_{i})$', va='center', ha='center', backgroundcolor='w')
-        ax.text(2, args.Vms[j], f'$U_{j}$', va='center', ha='center', backgroundcolor='w')
+        ax.text(
+            2,
+            weighted_potential,
+            f"$\\tilde U_{j} = U_{j} + w\\cdot(U_{j} - U_{i})$",
+            va="center",
+            ha="center",
+            backgroundcolor="w",
+        )
+        ax.text(
+            2, args.Vms[j], f"$U_{j}$", va="center", ha="center", backgroundcolor="w"
+        )
 
     ax.set_xlim(0, T)