Brunel network

In this tutorial, we will build a classic Brunel network [2] using LIF cells in Arbor, and you can compare it with the simulation using LIF cells in the NEST simulator.

Concepts covered in this example :class: note

1. Build LIF cells by loading certain neuron parameters from the parameter file.

2. Connect neurons in a fixed in-degree manner based on a connection probability.

3. Add Poissonian input to drive the network activity.

4. Record spikes and plot raster plot and peristimulus time histogram (PSTH).

Start a recipe and initiate the parameters

Here, we will follow the protocol for building the recipe, which has been discussed before in Ring Network.

# load packages needed for this project
import arbor as A
from arbor import units as U
import numpy as np
from numpy.random import RandomState
from parameters import *


class brunel_recipe(A.recipe):
    # define recipe for brunel network
    def __init__(
        self,
        NE,
        NI,
        epsilon,
        weight,
        delay,
        g,
        rate,
        seed,
        tau_m,
        V_th,
        V_m,
        E_L,
        V_reset,
        C_m,
        t_ref,
    ):
        A.recipe.__init__(self)

        self.seed_ = seed

        # set up neruon parameters
        self.tau_m = neuron_params["tau_m"] * U.ms
        self.V_th = neuron_params["V_th"] * U.mV
        self.V_m = neuron_params["V_m"] * U.mV
        self.E_L = neuron_params["E_L"] * U.mV
        self.E_R = neuron_params["V_reset"] * U.mV
        self.C_m = neuron_params["C_m"] * U.pF
        self.t_ref = neuron_params["t_ref"] * U.ms

        # set up network parameters
        self.ncells_exc_ = NE
        self.ncells_inh_ = NI
        self.in_degree_exc_ = round(
            epsilon * NE
        )  # fixed indegree based on connection probability
        self.in_degree_inh_ = round(
            epsilon * NI
        )  # fixed indegree based on connection probability
        self.delay_ = delay * U.ms

        # set and scale synaptic weights based on C_m due to the conductance based property of LIF neuron in Arbor
        self.weight_exc_ = weight * C_m
        self.weight_inh_ = -1 * g * weight * C_m
        self.weight_ext_ = weight * C_m

        # set the firing rate of the Poissonian input
        self.lambda_ = rate * U.Hz

We can find all the parameters listed in a separate parameter file:

# define network parameters
order = 2500
NE = order * 4
NI = order
g = 8.0
eta = 1.5
epsilon = 0.1

J = 0.1
weight = J
delay = 1.5
theta = 20.0

CE = round(epsilon * NE)
CI = round(epsilon * NI)

# define neuron parameter
tau_m = 20
V_th = theta
C_m = 250.0
E_L = 0.0
V_reset = 10.0
V_m = 0.0
t_ref = 2.0


# for nest
neuron_model = "iaf_psc_delta"

neuron_params = {
    "C_m": C_m,
    "tau_m": tau_m,
    "t_ref": t_ref,
    "E_L": E_L,
    "V_reset": V_reset,
    "V_m": V_m,
    "V_th": theta,
}


# define external input
nu_th = theta / (J * CE * tau_m)
nu_ex = eta * nu_th
rate = 1000.0 * nu_ex * CE

# simulation parameters ####
tfinal = 1000
dt = 1
seed = 42

We define the network size with num_cells and cell type with cell_kind. Then load all the neuron parameters to the LIF cells with cell_description function.

def num_cells(self):
    return self.ncells_exc_ + self.ncells_inh_

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

def cell_description(self, gid):
    return A.lif_cell(
        "src",
        "tgt",
        tau_m=self.tau_m,
        V_th=self.V_th,
        C_m=self.C_m,
        E_L=self.E_L,
        E_R=self.E_R,
        V_m=self.V_m,
        t_ref=self.t_ref,
    )

The Brunel network is randomly sparsely connected with a fixed in-degree regulated by a connection probability (\(\epsilon\)). We, therefore, define a function to enable random connectivity. This funciton draws random connections from the pre-selected pool of source neruons defined by the \(gid\) defined within \(start\) and \(end\). The total number of random connections are regulated by the fix-indegree value, which is \(m\) here in the function and \(CE\) and \(CI\) in the parameters file.

def sample_subset(self, gen, gid, start, end, m):
    # this function is used to generate random connection based on fixed indegree
    idx = np.arange(start, end)
    if start <= gid < end:
        idx = np.delete(idx, gid - start)
    gen.shuffle(idx)
    return idx[:m]

def connections_on(self, gid):
    gen = RandomState(gid + self.seed_)
    # Add incoming excitatory connections.
    connections = [
        A.connection((i, "src"), "tgt", self.weight_exc_, self.delay_)
        for i in self.sample_subset(
            gen, gid, 0, self.ncells_exc_, self.in_degree_exc_
        )
    ]
    # Add incoming inhibitory connections.
    connections += [
        A.connection((i, "src"), "tgt", self.weight_inh_, self.delay_)
        for i in self.sample_subset(
            gen,
            gid,
            self.ncells_exc_,
            self.ncells_exc_ + self.ncells_inh_,
            self.in_degree_inh_,
        )
    ]

    return connections

To enable the network activity, we apply Poissonian input via event_generators to to each neuron in the network. It aims to achieve a similar effect as the Poisson_generator in the NEST simulator.

def event_generators(self, gid):
    # add poissonian input to each neuron
    return [
        A.event_generator(
            "tgt",
            self.weight_ext_,
            A.poisson_schedule(
                freq=self.lambda_,
                seed=gid + (self.ncells_exc_ + self.ncells_inh_) * self.seed_,
            ),
        )
    ]

In the end, we build the network, run the simulation, and record the spikes.

if __name__ == "__main__":
    recipe = brunel_recipe(
        NE,
        NI,
        epsilon,
        weight,
        delay,
        g,
        rate,
        seed,
        tau_m,
        V_th,
        V_m,
        E_L,
        V_reset,
        C_m,
        t_ref,
    )

    context = A.context()
    if A.config()["profiling"]:
        A.profiler_initialize(context)
    print(context)
    meters = A.meter_manager()
    meters.start(context)

    hint = A.partition_hint()
    hint.cpu_group_size = 5000
    hints = {A.cell_kind.lif: hint}
    decomp = A.partition_load_balance(recipe, context, hints)
    print(decomp)

    sim = A.simulation(recipe, context, decomp)
    sim.record(A.spike_recording.all)
    sim.run(tfinal * U.ms, dt * U.ms)

    # get data
    spikes = sim.spikes()
    sources = spikes["source"]["gid"]
    times = spikes["time"]

    # save data
    np.save("times.npy", times)
    np.save("sources.npy", sources)

One can also use the code below to visualize the raster plot of the entire nework and a few selected cells, and the peristimulus time histogram (PSTH) of the entire network. The parameters used here are supposed to achieve asynchronous irregular dynamics.

import numpy as np
import matplotlib.pyplot as plt
from parameters import *


def plot_raster_and_PSTH(times, sources, bin=30):
    plt.figure(figsize=(60, 80))

    plt.subplot(311)
    plt.plot(times, sources, "|", color="k")
    plt.ylim([100, 200])
    plt.xlim(0, tfinal)
    plt.ylabel("neuron ID")

    plt.subplot(312)
    plt.plot(times, sources, "|", color="k")
    plt.xlim(0, tfinal)
    plt.ylabel("neuron ID")

    plt.subplot(313)
    counts = np.histogram(times, range(0, tfinal))[0]
    lefts = np.histogram(times, range(0, tfinal))[1][0:-1]
    plt.bar(lefts, counts, color="k")
    plt.xlabel("time | ms")
    plt.ylabel("counts")
    plt.xlim(0, tfinal)

    plt.suptitle("raster plot and PSTH")
    plt.savefig("dynamics.svg")


def analysis_rate(sources, times, tfinal):
    rates = []

    for neuron in np.arange(0, NE + NI, 1):
        idx = np.where(sources == neuron + 1)
        time = times[idx]
        rate = len(time) / (tfinal / 1000.0)
        rates.append(rate)
    return rates


def plot_rate_histogram(sources, times, tfinal):
    rates = analysis_rate(sources, times, tfinal)
    plt.figure(2)
    plt.hist(rates)
    plt.xlabel("firing rate | Hz")
    plt.ylabel("counts")
    plt.savefig("rates.svg")


times = np.load("times.npy")
sources = np.load("sources.npy")
plot_raster_and_PSTH(times, sources)
plot_rate_histogram(sources, times, tfinal)
plt.show()

Raster plot and PSTH of the Brunel network simulated in Arbor.

Firing rate distribution of all neurons in the Brunel network simulated in Arbor.

The full code ————- You can find the same network architecture simulated in the NEST simulator in the same repo python/examples/brunel/nest_brunel.py. The average firing rate of neurons and network dynamics look similar in both cases.

Raster plot and PSTH of the Brunel network simulated in NEST.

Firing rate distribution of all neurons in the Brunel network simulated in NEST.

References

[2]

Brunel, Journal of Computational Neuroscience 8: 183-208 (2000); https://link.springer.com/article/10.1023/A:1008925309027.