A ring network

In this example, a small network of cells arranged in a ring will be created, and the simulation will be distributed over multiple threads or GPUs, if available.

Note

Concepts covered in this example:

1. Building a basic arbor.cell with a synapse site and spike generator.

2. Building an arbor.recipe with a network of interconnected cells.

3. Running the simulation and extracting the results.

The cell

Step (1) shows a function that creates a simple cell with a dendrite. We construct the following morphology and label the soma and dendrite:

A 4-segment cell with a soma (pink) and a branched dendrite (light blue).

    # (1) Build a segment tree
    # The dendrite (dend) attaches to the soma and has two simple segments
    # attached.
    #
    #             left
    #            /
    # soma - dend
    #            \
    #             right
    tree = A.segment_tree()
    root = A.mnpos
    # Soma (tag=1) with radius 6 μm, modelled as cylinder of length 2*radius
    soma = tree.append(root, (-12, 0, 0, 6), (0, 0, 0, 6), tag=1)
    # Single dendrite (tag=3) of length 50 μm and radius 2 μm attached to soma.
    dend = tree.append(soma, (0, 0, 0, 2), (50, 0, 0, 2), tag=3)
    # Attach two dendrites (tag=3) of length 50 μm to the end of the first dendrite.
    # Radius tapers from 2 to 0.5 μm over the length of the dendrite.
    l = 50 / sqrt(2)
    _ = tree.append(
        dend,
        (50, 0, 0, 2),
        (50 + l, l, 0, 0.5),
        tag=3,
    )
    # Constant radius of 1 μm over the length of the dendrite.
    _ = tree.append(
        dend,
        (50, 0, 0, 1),
        (50 + l, -l, 0, 1),
        tag=3,
    )

In step (2) we create a label for both the root and the site of the synapse. These locations will form the endpoints of the connections between the cells.

We’ll create labels for the root (red) and a synapse_site (black).

    # Associate labels to tags
    labels = A.label_dict(
        {
            "soma": "(tag 1)",
            "dend": "(tag 3)",
            # (2) Mark location for synapse at the midpoint of branch 1 (the first dendrite).
            "synapse_site": "(location 1 0.5)",
            # Mark the root of the tree.
            "root": "(root)",
        }
    )

After we’ve created a basic arbor.decor, step (3) places a synapse with an exponential decay ('expsyn') on the 'synapse_site'. The synapse is given the label 'syn', which is later used to form arbor.connection objects terminating at the cell.

Note

Mechanisms can be initialized with their name; 'expsyn' is short for arbor.mechanism('expsyn'). Mechanisms typically have some parameters, which can be queried (see arbor.mechanism_info) and set (see arbor.mechanism). In particular, the e parameter of expsyn defaults to 0, which makes it, given the typical resting potential of cell membranes of -70 mV, an excitatory synapse.

Step (4) places a threshold detector at the 'root'. The detector is given the label 'detector', which is later used to form arbor.connection objects originating from the cell.

Note

The number of synapses placed on the cell in this case is 1, because the 'synapse_sites' locset is an explicit location. Had the chosen locset contained multiple locations, an equal number of synapses would have been placed, all given the same label 'syn'.

The same explanation applies to the number of detectors on this cell.

    # (3) Create a decor and a cable_cell
    decor = (
        A.decor()
        # Put hh dynamics on soma, and passive properties on the dendrites.
        .paint('"soma"', A.density("hh")).paint('"dend"', A.density("pas"))
        # (4) Attach a single synapse.
        .place('"synapse_site"', A.synapse("expsyn"), "syn")
        # Attach a detector with threshold of -10 mV.
        .place('"root"', A.threshold_detector(-10 * U.mV), "detector")
    )

The recipe

# (5) Create a recipe that generates a network of connected cells.
class ring_recipe(A.recipe):
    def __init__(self, ncells):
        # Base class constructor must be called first for proper initialization.
        A.recipe.__init__(self)
        self.ncells = ncells
        self.props = A.neuron_cable_properties()

    # (6) Returns the total number of cells in the model; must be implemented.
    def num_cells(self):
        return self.ncells

    # (7) The cell_description method returns a cell
    def cell_description(self, gid):
        return make_cable_cell(gid)

    # Return the type of cell; must be implemented and match cell_description.
    def cell_kind(self, _):
        return A.cell_kind.cable

    # (8) For each gid, provide a list of incoming connections.
    def connections_on(self, gid):
        # This defines the ring by connecting from the last gid. The first src
        # comes from the _last_ gid, closing the ring.
        src = (gid - 1) % self.ncells
        w = 0.01  # 0.01 μS on expsyn
        d = 5 * U.ms
        return [A.connection((src, "detector"), "syn", w, d)]

    # (9) Attach a generator to the first cell in the ring.
    def event_generators(self, gid):
        if gid == 0:
            sched = A.explicit_schedule([1 * U.ms])  # one event at 1 ms
            weight = 0.1  # 0.1 μS on expsyn
            return [A.event_generator("syn", weight, sched)]
        return []

    # (10) Place a probe at the root of each cell.
    def probes(self, gid):
        return [A.cable_probe_membrane_voltage('"root"', "Um")]

    def global_properties(self, _):
        return self.props

To create a model with multiple connected cells, we need to use a recipe. The recipe is where the different cells and the connections between them are defined.

Step (5) shows a class definition for a recipe with multiple cells. Instantiating the class requires the desired number of cells as input. Compared to the simple cell recipe, the main differences are connecting the cells (8), returning a configurable number of cells (6) and returning a new cell per gid (7).

Step (8) creates an arbor.connection between consecutive cells. If a cell has gid gid, the previous cell has a gid (gid-1)%self.ncells. The connection has a weight of 0.01 (inducing a conductance of 0.01 μS in the target mechanism expsyn) and a delay of 5 ms. The first two arguments to arbor.connection are the source and target of the connection.

The source is an arbor.cell_global_label object containing a cell index gid, the source label corresponding to a valid detector label on the cell and an optional selection policy (for choosing a single detector out of potentially many detectors grouped under the same label - remember, in this case the number of detectors labeled ‘detector’ is 1). The arbor.cell_global_label can be initialized with a (gid, label) tuple, in which case the selection policy is the default arbor.selection_policy.univalent; or a (gid, (label, policy)) tuple.

The target is a arbor.cell_local_label object containing a cell index gid, the target label corresponding to a valid synapse label on the cell and an optional selection policy (for choosing a single synapse out of potentially many synapses grouped under the same label - remember, in this case the number of synapses labeled ‘syn’ is 1). The arbor.cell_local_label can be initialized with a label string, in which case the selection policy is the default arbor.selection_policy.univalent; or a (label, policy) tuple. The gid of the target cell doesn’t need to be explicitly added to the connection, it is the argument to the arbor.recipe.connections_on() method.

Step (9) attaches an arbor.event_generator on the 0th target (synapse) on the 0th cell; this means it is connected to the "synapse_site" on cell 0. This initiates the signal cascade through the network. The arbor.explicit_schedule in instantiated with a list of times in milliseconds, so here a single event at the 1 ms mark is emitted. Note that this synapse is connected twice, once to the event generator, and once to another cell.

Step (10) places a probe at the "root" of each cell.

# (11) Instantiate recipe
ncells = 4
recipe = ring_recipe(ncells)

Step (11) instantiates the recipe with 4 cells.

The execution

To create a simulation, we need at minimum to supply the recipe, and in addition can supply a arbor.context and arbor.domain_decomposition. The first lets Arbor know what hardware it should use, the second how to destribute the work over that hardware. By default, contexts are configured to use 1 thread and domain decompositons to divide work equally over all threads.

# (12) Create a simulation using the default settings:
# - Use all threads available
# - Use round-robin distribution of cells across groups with one cell per group
# - Use GPU if present
# - No MPI
# Other constructors of simulation can be used to change all of these.
sim = A.simulation(recipe)

# (13) Set spike generators to record
sim.record(A.spike_recording.all)

# (14) Attach a sampler to the voltage probe on cell 0. Sample rate of 10 sample every ms.
handles = [
    sim.sample((gid, "Um"), A.regular_schedule(0.1 * U.ms)) for gid in range(ncells)
]

# (15) Run simulation for 100 ms
sim.run(100 * U.ms)
print("Simulation finished")

Step (12) creates a simulation object from the recipe. Optionally, the simulation constructor takes two more parameters: a arbor.context and a arbor.domain_decomposition. In a followup of this tutorial that will be demonstrated. For now, it is enough to know that for simulations that don’t require customized execution those arguments can be left out. Without further arguments Arbor will use all locally available threads.

Step (13) sets all spike generators to record using the arbor.spike_recording.all policy. This means the timestamps of the generated events will be kept in memory. Be default, these are discarded.

In addition to having the timestamps of spikes, we want to extract the voltage as a function of time.

Step (14) sets the probes (step 10) to measure at a certain schedule. This is sometimes described as attaching a sampler to a probe. arbor.simulation.sample() expects a probeset id and the desired schedule (here: a recording frequency of 10 kHz, or a dt of 0.1 ms). Note that the probeset id is a separate index from those of connection endpoints; probeset ids correspond to the index of the list produced by arbor.recipe.probes() on cell gid.

arbor.simulation.sample() returns a handle to the samples that will be recorded. We store these handles for later use.

Step (15) executes the simulation for a duration of 100 ms.

The results

Step (16) prints the timestamps of the spikes:

# (16) Print spike times
print("spikes:")
for sp in sim.spikes():
    print(" ", sp)

Step (17) generates a plot of the sampling data. arbor.simulation.samples() takes a handle of the probe we wish to examine. It returns a list of (data, meta) terms: data being the time and value series of the probed quantity; and meta being the location of the probe. The size of the returned list depends on the number of discrete locations pointed to by the handle, which in this case is 1, so we can take the first element. (Recall that in step (10) we attached a probe to the "root", which describes one location. It could have described a locset.)

# (17) Plot the recorded voltages over time.
print("Plotting results ...")
dfs = []
for gid in range(ncells):
    samples, meta = sim.samples(handles[gid])[0]
    dfs.append(
        pandas.DataFrame(
            {"t/ms": samples[:, 0], "U/mV": samples[:, 1], "Cell": f"cell {gid}"}
        )
    )
df = pandas.concat(dfs, ignore_index=True)
seaborn.relplot(
    data=df, kind="line", x="t/ms", y="U/mV", hue="Cell", errorbar=None
).savefig("network_ring_result.svg")

Since we have created ncells cells, we have ncells traces. We should be seeing phase shifted traces, as the action potential propagated through the network.

We plot the results using pandas and seaborn:

The full code

You can find the full code of the example at python/example/network_ring.py.