Contents

2024-07-01 · AdExNet-Brian

2024-07-01 · AdExNet-Brian

Using Brian to simulate a network of spiking neurons (specifically: conductance-based (COBA) adaptive exponential (AdEx) integrate and fire neurons).

from brian2 import *

Neuron model

Copying from https://github.com/tfiers/phd/blob/main/nb/pylib/neuron.py :

Parameters for a cortical regular spiking (RS) neuron, from Naud 2008.

C   = 104  * pF
gL  = 4.3  * nS
EL  = -65  * mV
VT  = -52  * mV
DT  = 0.8  * mV
Vs  =  40  * mV
Vr  = -53  * mV
a   = -0.8 * nS
b   =  65  * pA
tau_w = 88 * ms

Ee  =   0  * mV
Ei  = -80  * mV
tau_g = 7  * ms;

Strength of noisy external input current:

Iext = 8 * pA * sqrt(second)
\[8e-12\,s^{0.5}\,A\]
eqs = """
dV/dt = ( -gL*(V - EL) + gL * DT * exp((V-VT)/DT) -I -w -(Iext * xi)) / C : volt
dw/dt = (a*(V - EL) - w) / tau_w : amp

I = ge * (V - Ee) + gi * (V - Ei) : amp

dge/dt = -ge / tau_g : siemens
dgi/dt = -gi / tau_g : siemens
""";

def COBA_AdEx_neuron(N = 1):
    G = NeuronGroup(
        N, eqs, threshold="V > Vs", reset="V = Vr; w += b", method='euler',
    )
    G.V = EL
    # Rest of vars are auto set to 0
    return G

COBA_AdEx_neuron()
NeuronGroup 'neurongroup' with 1 neurons.
Model: \begin{align*}I &= ge \left(- Ee + V\right) + gi \left(- Ei + V\right) && \text{(unit of $I$: $\mathrm{A}$)}\\ \frac{\mathrm{d}V}{\mathrm{d}t} &= \frac{DT gL e^{\frac{V - VT}{DT}} - I - Iext \xi - gL \left(- EL + V\right) - w}{C} && \text{(unit of $V$: $\mathrm{V}$)}\\ \frac{\mathrm{d}ge}{\mathrm{d}t} &= - \frac{ge}{\tau_{g}} && \text{(unit of $ge$: $\mathrm{S}$)}\\ \frac{\mathrm{d}gi}{\mathrm{d}t} &= - \frac{gi}{\tau_{g}} && \text{(unit of $gi$: $\mathrm{S}$)}\\ \frac{\mathrm{d}w}{\mathrm{d}t} &= \frac{a \left(- EL + V\right) - w}{\tau_{w}} && \text{(unit of $w$: $\mathrm{A}$)}\end{align*} Spiking behaviour:
  • Threshold condition: V > Vs
  • Reset statement(s): V = Vr; w += b

Num neurons.

N = 100;

Number of excitatory neurons:

Ne = N * 4//5

80

G = COBA_AdEx_neuron(N)  # NeuronGroup

Ge = G[:Ne]
Gi = G[(Ne+1):];

Connectivity

delay = 0 * ms

p_conn = 0.05

we = 10 * nS / (N*p_conn)  # Based on pkgs Nto1AdEx, and VoltoMapSim

wi = 4*we;

seed(1234)

Se = Synapses(Ge, G, on_pre="ge += we", delay=delay)
Si = Synapses(Gi, G, on_pre="gi += wi", delay=delay);

prefs.codegen.target = "numpy"  # Until we have enough internet to download C++ build tools ('visual studio').

Se.connect(condition='i != j', p=p_conn)
Si.connect(condition='(i+Ne+1) != j', p=p_conn)

Convert to connection matrix.

W = np.zeros((N,N), dtype=int)

for (i,j) in zip(Se.i, Se.j):
    W[i,j] = 1

for (i,j) in zip(Si.i, Si.j):
    i += (Ne+1)  # Offset for the inhibitory 'from' neurons
    W[i,j] = -4

W

array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ...,
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]])

A sanity check:

np.count_nonzero(W) / np.size(W)

0.0488

Save to disk.

data_dir = "../data/2024-07-01__AdExNet-Brian/";

import os
os.makedirs(data_dir, exist_ok=True)

np.save(data_dir + "connectivity_matrix.npy", W)

Simulate net

M = StateMonitor(G, ["V"], record=[0])
S = SpikeMonitor(G)
net = Network(G, Se, Si, M, S);

T = 60 * second;

net.run(T, report='stdout', report_period=10*second)

Starting simulation at t=0. s for a duration of 60. s
7.8048 s (13%) simulated in 10s, estimated 1m 7s remaining.
16.2596 s (27%) simulated in 20s, estimated 54s remaining.
24.8611 s (41%) simulated in 30s, estimated 42s remaining.
33.5618 s (55%) simulated in 40s, estimated 32s remaining.
42.116 s (70%) simulated in 50s, estimated 21s remaining.
50.6588 s (84%) simulated in 1m 0s, estimated 11s remaining.
59.2927 s (98%) simulated in 1m 10s, estimated 1s remaining.
60. s (100%) simulated in 1m 10s

Save spike trains to disk.

trains = {}
for (i,train) in S.spike_trains().items():
    trains[f"neuron_{i}"] = train / second

trains["neuron_0"]

array([ 0.0981,  0.2154,  0.4087,  0.8497,  0.8514,  1.2013,  1.5335,
        1.5479,  1.7766,  1.9157,  2.2367,  2.3956,  2.6073,  2.6266,
        2.8691,  2.8704,  3.2107,  3.6489,  3.7416,  4.0888,  4.7359,
        4.7523,  5.2979,  5.5197,  5.6059,  5.6559,  6.0879,  6.3774,
        6.6031,  6.6055,  6.6094,  7.0454,  7.4738,  7.733 ,  8.1622,
        8.1809,  8.9604,  9.231 ,  9.787 , 10.0557, 10.3248, 10.5081,
       11.0537, 11.3813, 11.3848, 11.5985, 11.9418, 12.0614, 12.5731,
       13.0706, 13.1415, 13.4369, 13.7976, 14.2316, 14.6229, 14.8934,
       15.121 , 15.1301, 15.1887, 15.1938, 15.9013, 15.9265, 16.2067,
       16.2113, 16.3349, 16.5827, 16.8747, 17.319 , 17.3577, 17.5984,
       17.9209, 17.9249, 17.9301, 18.4096, 18.5597, 18.6846, 19.0331,
       19.2043, 19.5641, 19.8067, 19.8588, 19.8636, 20.2356, 20.2536,
       20.7705, 21.2847, 21.5383, 21.9267, 21.9302, 22.081 , 22.4182,
       22.8306, 22.9873, 23.0414, 23.3584, 23.6576, 23.7582, 23.948 ,
       23.9501, 24.169 , 24.5981, 24.6032, 24.9426, 25.0947, 25.6163,
       25.6229, 26.1376, 26.1411, 26.3475, 26.3511, 26.7188, 27.3444,
       27.4142, 27.7603, 27.9331, 28.3518, 28.3528, 28.7086, 28.9969,
       29.0029, 29.2595, 29.7441, 29.8995, 30.0619, 30.6914, 30.9326,
       31.1707, 31.1772, 31.5057, 31.6431, 31.8948, 32.0425, 32.239 ,
       32.6661, 33.3309, 33.4437, 33.7612, 33.9635, 34.2051, 34.3509,
       34.356 , 34.3588, 34.8407, 35.0223, 35.1304, 35.5372, 35.6959,
       35.9651, 35.9821, 36.4458, 36.4536, 36.4549, 36.8519, 37.2377,
       37.4299, 37.5938, 37.7423, 37.8918, 38.1242, 38.1255, 38.4648,
       38.5299, 38.8984, 39.1581, 39.1602, 39.3277, 39.8567, 40.1646,
       40.1669, 40.5115, 40.5466, 40.7938, 40.9032, 41.1605, 41.6738,
       41.8306, 41.8332, 42.1683, 42.5166, 42.5207, 42.5279, 42.7583,
       42.7632, 43.6027, 43.9832, 44.0637, 44.308 , 44.5589, 44.5612,
       45.0076, 45.3233, 45.7163, 45.7292, 46.0697, 46.0956, 46.4798,
       47.0739, 47.0793, 47.2125, 47.4582, 47.7903, 47.8063, 48.3805,
       48.6393, 48.8619, 48.8832, 48.8922, 49.4908, 49.4942, 49.4957,
       49.4983, 49.7505, 50.206 , 50.4721, 50.4893, 50.8539, 50.926 ,
       51.1346, 51.3462, 51.3499, 51.835 , 51.9299, 51.9425, 52.2705,
       52.2719, 52.6566, 53.0916, 53.6611, 53.9234, 54.0125, 54.2602,
       54.3598, 54.4548, 54.9532, 54.9552, 55.9542, 55.9814, 56.2207,
       56.3134, 56.5939, 56.5981, 56.9588, 57.6817, 57.965 , 58.1998,
       58.3994, 58.6284, 58.937 , 59.4005, 59.5231, 59.5373])

np.savez(data_dir + "spiketimes.npz", **trains)

%run pylib/plot.py

importing mpl … ✔
importing brian … ✔

def ceil_spikes_brian(M: StateMonitor, S: SpikeMonitor, var='V', n=0, V_ceil=Vs):
    "For nice plots, set the voltage trace to some constant at spike times"
    V = getattr(M, var)[n]
    spikes = S.t[S.i == n]
    return ceil_spikes(V, M.t, spikes, V_ceil)

def ceil_spikes(V, t, spiketimes, V_ceil=Vs):
    "For nice plots, set the voltage trace to some constant at spike times"
    i = searchsorted(t, spiketimes)
    V[i] = V_ceil
    return V

Voltage of one (excitatory) neuron:

plotsig(ceil_spikes_brian(M, S), tlim=[0, 2]*second);
../_images/2024-07-01__AdExNet-Brian_39_0.png

Rasterplot:

tmax = 10*second
sel = S.t < tmax
ax = plot(
    S.t[sel], S.i[sel], "k.", xlim=[0, tmax], ylim=[0, N], ms=3,
    ylabel="Neuron number", xlabel="Simulation time",
)
ax.set_title(f"Spiketimes of {Ne} excitatory and {N - Ne} inhibitory neurons", fontsize="small");
../_images/2024-07-01__AdExNet-Brian_41_0.png

previous

Voltage-to-wiring, the simulation

next

2024-05-29 · Residuals linefit