Declarative Connectivity in Arbor

Concepts and Requirements

We will assume that you have read the basic recipe and network tutorials.

In addition to Arbor and its requirements matplotlib and networkx need to be installed.

In this tutorial, we are going to demonstrate how to leverage Arbor’s declarative connection description facilities to generate a few common network types. We will gradually build up complexity and generally show the full recipe first before discussing some of the relevant parts. High-level connectivity descriptions can be more intuitive for some types of networks as well as more performant in Python simulations, as the construction is handled entirely in C++.

Prelude: Unconnected Cells

Before we start building actual networks, we will set up the trivial network, which has no connections at all. This will be written as a recipe class, as discussed in other tutorials, and later examples will derive from this class to build upon. If you want, you can skip this part and come back as needed.

Our cells comprise a simple leaky integrate and fire model, not a cable cell, as we want to emphasize building networks. We begin by defining the global settings:

# global parameters
# cell count
N = 5
# total runtime [ms]
T = 1000
# numerical time step [ms]
dt = 0.1

  • N is the cell count of the simulation.

  • T is the total runtime of the simulation in ms.

  • dt is the numerical timestep on which cells evolve.

These parameters are used here:

if __name__ == "__main__":
    rec = unconnected(N)
    sim = A.simulation(rec)
    sim.record(A.spike_recording.all)
    sim.run(T * U.ms, dt * U.ms)
    plot_spikes(sim, T, N, prefix="01-")

where we run the simulation. Before we discuss the relevant details, the recipe reads in full

class unconnected(A.recipe):
    def __init__(self, N) -> None:
        super().__init__()
        self.N = N
        # Cell prototype
        self.cell = A.lif_cell("src", "tgt")
        # random seed [0, 100]
        self.seed = 42
        # event generator parameters
        self.gen_weight = 20
        self.gen_freq = 2.5 * U.kHz

    def num_cells(self) -> int:
        return self.N

    def event_generators(self, gid: int):
        if gid >= 1:
            return []
        seed = self.seed + gid * 100
        return [
            A.event_generator(
                "tgt",
                self.gen_weight,
                A.poisson_schedule(freq=self.gen_freq, seed=seed),
            )
        ]

    def cell_description(self, gid: int):
        return self.cell

    def cell_kind(self, gid: int) -> A.cell_kind:
        return A.cell_kind.lif

In the recipe, we set a prototypical LIF cell:

        self.cell = A.lif_cell("src", "tgt")

and deliver it for all gid:

    def cell_description(self, gid: int):
        return self.cell

With large and complicated cells this can sometimes help with performance, here, it’s just a convenient way to structure our recipe. Also, the first cell has an event generator attached, a Poisson point process seeded with the cell’s gid.

    def event_generators(self, gid: int):
        if gid >= 1:
            return []
        seed = self.seed + gid * 100
        return [
            A.event_generator(
                "tgt",
                self.gen_weight,
                A.poisson_schedule(freq=self.gen_freq, seed=seed),
            )
        ]

All other parameters are set in the constructor:

    def __init__(self, N) -> None:
        super().__init__()
        self.N = N
        # Cell prototype
        self.cell = A.lif_cell("src", "tgt")
        # random seed [0, 100]
        self.seed = 42
        # event generator parameters
        self.gen_weight = 20
        self.gen_freq = 2.5 * U.kHz

We also proceed to add spike recording and generate raster plots using a helper function plot_spikes from util.py, which results in

../_images/01-raster.svg

As only the first cell receives spiking inputs, only it will show up on the plot.

Ring Network

Starting from an unconnected set of cells, we can now start building a simple network. A ring structure is defined by connecting each cell to its predecessor, i.e. the cell with gid = i is connected to the cell with gid = i - 1 and the cell gid = 0 is connected to the last cell gid = N - 1.

We construct such a network by defining a new recpi ring deriving from the unconnected network

#!/usr/bin/env python3

import arbor as A
from arbor import units as U
from util import plot_spikes, plot_network
from unconnected import unconnected


# global parameters
# cell count
N = 5
# total runtime [ms]
T = 1000
# numerical time step [ms]
dt = 0.1


class ring(unconnected):
    def __init__(self, N) -> None:
        super().__init__(N)

    def network_description(self):
        # network structure
        wraps = f"(intersect (source-cell {self.N - 1}) (target-cell 0))"
        cells = f"(gid-range 0 {self.N})"
        chain = f"(chain {cells})"
        ring = f"(join {chain} {wraps})"
        # parameters
        weight = "(scalar 199.99999219)"
        delay = "(scalar 0.5)"
        return A.network_description(ring, weight, delay, {})


if __name__ == "__main__":
    ctx = A.context()
    rec = ring(N)
    sim = A.simulation(rec, ctx)
    sim.record(A.spike_recording.all)
    sim.run(T * U.ms, dt * U.ms)
    plot_spikes(sim, T, N, prefix="02-")
    plot_network(rec, prefix="02-")

The idiomatic way of extending classes with new functionality is to use inheritance. Importantly, the burden of initializing the base class falls on the derived class:

class ring(unconnected):
    def __init__(self, N) -> None:
        super().__init__(N)

Next, we add a new method that is responsible for the network. Note that this — in contrast to most other methods on recipe — does not have an argument of gid, since it is definining the global network.

    def network_description(self):
        # network structure
        wraps = f"(intersect (source-cell {self.N - 1}) (target-cell 0))"
        cells = f"(gid-range 0 {self.N})"
        chain = f"(chain {cells})"
        ring = f"(join {chain} {wraps})"
        # parameters
        weight = "(scalar 199.99999219)"
        delay = "(scalar 0.5)"
        return A.network_description(ring, weight, delay, {})

Similar to the construction of a decor or cv_policy, a light-weight language inspired by LISP or Scheme is used here. For this tutorial, we use Python format strings to compose expressions. Networks comprise a structure and parameters — weight and delay, which can be scalars as shown, or more elaborate expressions, such a drawing from a random distribution.

The structure is defined in terms of combinators reminiscent of relational algebra queries operating on abstract sets of source and target identifiers.

  • source-cell and target-cell construct a set of eligble sources and targets from a single gid.

  • gid-range defines a contiguous range of gids [start, end)

  • intersect constructs the connections between the source and target arguments.

  • chain constructs the set of connections between adjacent neighbours.

  • join takes two sub-structures A and B and returns their union.

Upon close inspection, these combinators directly spell out the prose description of the ring network given above: Connect adjacent cells and close the ring by connecting the beginning and end! Running the network and plotting the spikes we find cells deeper into the ring spiking now

../_images/02-raster.svg

The network structure is rendered via networkx

../_images/02-graph.svg

Excercise: All-to-all Network

Using the unconnected recipe and the network documentation define a fully connected network, i.e. where each cell is connected to every other cell except itself.

Hint

  1. source-cell and target-cell can take a range of ids

  2. Use and intersection with inter-cell to remove self connections

Our solution produces the following output

../_images/03-raster.svg

The network should look like this

../_images/03-graph.svg

For reference, we reproduce it here:

#!/usr/bin/env python3

import arbor as A
from arbor import units as U
from util import plot_spikes, plot_network
from unconnected import unconnected


# global parameters
# cell count
N = 5
# total runtime [ms]
T = 1000
# numerical time step [ms]
dt = 0.1


class all2all(unconnected):
    def __init__(self, N) -> None:
        super().__init__(N)

    def network_description(self):
        # network structure
        full = f"(intersect (inter-cell) (source-cell (gid-range 0 {self.N})) (target-cell (gid-range 0 {self.N})))"
        # parameters
        weight = "(scalar 125)"
        delay = "(scalar 0.5)"
        return A.network_description(full, weight, delay, {})


if __name__ == "__main__":
    rec = all2all(N)
    sim = A.simulation(rec)
    sim.record(A.spike_recording.all)
    sim.run(T * U.ms, dt * U.ms)
    plot_spikes(sim, T, N, prefix="03-")
    plot_network(rec, prefix="03-")

Brunel Network

The Brunel network, or in other words, a inhibition-dominated randomly connected recurrent network, is a common network structure used in computational neuroscience proposed by Nicolas Brunel in 2000. It contains sparsely connected inhibitory and excitatory neurons where a critical balance between inhibition and excitation inputs to each neuron is maintained to ensure a brain-realistic network-wide dynamics. It entails a few typical dynamics of cortical circuits.

Practically, we can describe this network by two populations, called the excitatory and inhibitory populations, such that

  1. Each cell is connected to each other with some probablity \(0 < p < 1\) - There no self-connections

  2. If the pre-synaptic cell is in the excitatory population, the weight is \(w_{exc} > 0\)

  3. If the pre-synaptic cell is in the inhitatory population, the weight is \(w_{inh} < 0\) - \(|w_{inh}| < |w_{exc}|\)

The Brunel network simulation can be implemented like this

class brunel(unconnected):
    def __init__(self, N) -> None:
        super().__init__(N)
        # excitatory population: first 80% of the gids
        self.n_exc = int(0.8 * N)
        # inhibitory population: the remainder
        self.n_inh = N - self.n_exc
        # seed for random number generation
        self.seed = 42
        # excitatory weight
        self.weight = 100
        # relative weight of inhibitory connections
        self.g = 0.8
        # probability of connecting any two neurons
        self.p = 0.1

again using the base class unconnected to define everything except the network. We implement these by writing down the rules above in the recipe

    def network_description(self):
        rand = f"""(intersect (inter-cell)
                                (random {self.seed} {self.p}))"""
        inh = f"(gid-range {self.n_exc} {self.N})"
        weight = f"""(if-else (source-cell {inh})
                              (scalar {self.g * self.weight})
                              (scalar {self.weight}))"""
        delay = "(scalar 0.5)"
        return A.network_description(rand, weight, delay, {})

The rand structure encodes the random connectivity and removes any potential self-connections by intersect with inter-cell, as before. Next, we define the weight according the population of the pre-synaptic neuron. The population is defined by the gid of the neuron; the first 80% of cells is considered excitatory and the remainder inhibitory. The predicate inh reifies this description. The weight function weight then dispatches to one of two values based on the predicate.

Rendering the structure becomes slow and frankly unusable, but showing the adjacency matrix might be helpful

../_images/04-matrix.svg

Note that rendering can be disabled, if things get too slow.

Final Thoughts

Using a few examples we have shown how Arbor’s high-level network description method can be leveraged to generate common structures. The key insight is to build complex layouts from atomic blocks by using set operators like join, difference, and intersect. There are more to explore in the documentation, be especially aware of stochastic distributions. We have also seen how to produce weights and by extension delays using the same, declarative approach. This functionality is quite young and if any useful additions come to mind, do not hesitate to request or implement them!