From aeb54ed6bcad6498fd39f350597a1fc086e538db Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Thu, 19 Oct 2023 14:24:59 +0200 Subject: [PATCH] Some more polish. --- .../network_two_cells_gap_junctions.rst | 29 +++++------- .../network_two_cells_gap_junctions.py | 47 +++++++++++-------- 2 files changed, 41 insertions(+), 35 deletions(-) diff --git a/doc/tutorial/network_two_cells_gap_junctions.rst b/doc/tutorial/network_two_cells_gap_junctions.rst index 13516ae331..8e1293ec80 100644 --- a/doc/tutorial/network_two_cells_gap_junctions.rst +++ b/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 diff --git a/python/example/network_two_cells_gap_junctions.py b/python/example/network_two_cells_gap_junctions.py index b5730bc08b..be75d732fe 100755 --- a/python/example/network_two_cells_gap_junctions.py +++ b/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)