# Modeling spiking neural networks with Brian


```python
# If Brian2 is not yet installed, the following will install it
# (otherwise it will print a number of "Requirement already satisfied" lines)
%pip install brian2
```

    Requirement already satisfied: brian2 in /home/marcel/anaconda3/envs/datasci_course/lib/python3.8/site-packages (2.4.1)
    Requirement already satisfied: jinja2>=2.7 in /home/marcel/anaconda3/envs/datasci_course/lib/python3.8/site-packages (from brian2) (2.11.2)
    Requirement already satisfied: sympy>=1.2 in /home/marcel/anaconda3/envs/datasci_course/lib/python3.8/site-packages (from brian2) (1.6.2)
    Requirement already satisfied: setuptools>=24.2 in /home/marcel/anaconda3/envs/datasci_course/lib/python3.8/site-packages (from brian2) (49.6.0.post20201009)
    Requirement already satisfied: numpy>=1.15 in /home/marcel/anaconda3/envs/datasci_course/lib/python3.8/site-packages (from brian2) (1.19.2)
    Requirement already satisfied: pyparsing in /home/marcel/anaconda3/envs/datasci_course/lib/python3.8/site-packages (from brian2) (2.4.7)
    Requirement already satisfied: cython>=0.29 in /home/marcel/anaconda3/envs/datasci_course/lib/python3.8/site-packages (from brian2) (0.29.21)
    Requirement already satisfied: MarkupSafe>=0.23 in /home/marcel/anaconda3/envs/datasci_course/lib/python3.8/site-packages (from jinja2>=2.7->brian2) (1.1.1)
    Requirement already satisfied: mpmath>=0.19 in /home/marcel/anaconda3/envs/datasci_course/lib/python3.8/site-packages (from sympy>=1.2->brian2) (1.1.0)
    Note: you may need to restart the kernel to use updated packages.


Let's first import "everything" from the Brian 2 package.

We also switch off Brian's "code generation" mechanism, which improves the performance for complex models by generating/compiling/executing C++ code "behind the scenes". For the simple models that we are covering in this tutorial, this is not necessary – it even slows down things due to the need for compilation.


```python
from brian2 import *
prefs.codegen.target = 'numpy'
```

We ask the notebook to include plots directly in the notebook (instead of showing them in a separate window). Note that lines starting with % are specific commands for the jupyter notebook, they won't work in a Python script, for example. We also reduce the length of error messages which tend to get very long in Python.


```python
import matplotlib.pyplot as plt
import numpy as np

%matplotlib inline
%xmode minimal
```

    Exception reporting mode: Minimal


## Brian's unit system

Brian has a system for physical units and handles conversions between scales as well as the combination of different units automatically. All units can use a long form (e.g. `mvolt` or `nsiemens`), and for the most common units also a short form is available (e.g. `mV`, `nS`, `ms`).


```python
x = -60*mV
y = 10*nS
print(x, x * y)
```

    -60. mV -0.6 nA


Sometimes we need to remove the units, e.g. to hand them over to matplotlib's plotting functions which do not know how to handle Brian's unit system. The easiest way to do this is to divide by the scale we are interested in.


```python
print(x/mV)  # "value of x in mV"
print(y/nS)  # "value of y in nS"
print(y/uS)  # "value of y in µS"
```

    -60.0
    10.0
    0.01


Brian checks the consistency of units, which can be very useful in particular when writing complex equations. If the units do not match, an error is raised:


```python
print(x + y)  # does not make sense
```

## Neurons

Let's create a group of integrate-and-fire neurons. We use the same equation as in last week's tutorial
$$
\tau\frac{\mathrm{d}V}{\mathrm{d}t} = E_L - V + I_\mathrm{stim}
$$
If $V > V_\mathrm{threshold}$: emit a spike and set $V \leftarrow V_\mathrm{reset}$.


```python
# The "start_scope()" call is necessary if we define several models in the same script or notebook.
# It tells Brian to forget about any model defined previously. Otherwise, when we call "run(...)"
# Brian would not know which model to run.
start_scope()

# Constant that are common to all neurons (or later, synapses), can be defined as Python variables
N = 3
E_L = -50*mV
V_threshold = -30*mV
V_reset = -55*mV
C = 70*pF
g_L = 10*nS
tau = C/g_L
duration = 100*ms

# A NeuronGroup is a group of neurons that are defined by the same equations. The neurons do not need
# to be all identical, they can have different parameters and/or different initial values for their
# variables. In the following model, the variable V is defined as a differential equation, and in
# addition we add a parameter "I_stim" which is different for each neuron.
# Note that we have to add unit definitions like ": volt" at the end of each variable definition, so
# that Brian can verify that all units are consistent.
neurons = NeuronGroup(N, '''dV/dt = (E_L - V + I_stim)/tau : volt
                            I_stim : volt (constant)  # note, this is rescaled with g_L already''',
                      threshold='V>V_threshold', reset='V=V_reset', method='euler')

# All variables of the model (here: V and I_stim) are available as attributes of the NeuronGroup
# object "neuron". We can get and set their values by using neuron.<variable_name> .
# Initialize values
neurons.V = V_reset  # same initial value for all neurons
neurons.I_stim = [19, 21, 30]*mV  # one value for each neuron

# After running a model, we have access to the values of all variables at the end of the simulation.
# In most simulations, however, we are also interested in the values of the variables *during* the
# simulation, i.e. their evolution over time. In order to record these values, we need to create
# so called monitors. A "StateMonitor" records the values of variables over time, and a "SpikeMonitor"
# records action potential (which neuron spiked when)

# Record variable "V" in group "neuron". The argument "record=True" is a shorthand for "record the values
# from all neurons in the group". Alternatively, we could give the indices of the neurons we are interested
# in (e.g. "record=[0, 2, 4, 100]" to record from the neurons with the given indices)
v_mon = StateMonitor(neurons, 'V', record=True)

# Record spikes in the group neuron
spike_mon = SpikeMonitor(neurons)

# run simulation for 100ms
run(duration)
```

We could plot things directly by calling `plt.plot`, but by calling `plt.subplots` first we can easily arrange things in subplots.


```python
fig, ax = plt.subplots(2, 1, sharex=True)
# The recorded values are accessible as attributes of the monitor objects.
# A SpikeMonitor provides the attribute "t" for the time of each spike and the attribute "i" for the
# index of the neuron that spiked. E.g. if spike_mon.t returns [10*ms, 20*ms] and spike_mon.i returns
# [0, 5], this means that neuron 0 spiked at 10ms, and neuron 5 spiked at 20ms. These two vectors are
# what we need to make a raster plot of the spiking activity: we plot the time (in ms) on the x axis
# and the indices of the neurons on the y axis:
ax[0].plot(spike_mon.t/ms, spike_mon.i, 'o')
ax[0].set_ylim(-0.1, 2.1)  # Otherwise matplotlib would not show the empty row

# For a StateMonitor, the "t" attribute contains the times when values were recorded (by default at 
# every time step of the simulation). Since we recorded the "V" variable, the StateMonitor also has
# a V attribute containing the values of V for all neurons for each time steps as a 2D matrix. For
# plotting, this matrix is in the wrong sense: it has one row for each neuron and one column for each
# time step, but we need one row for each time step and one column for each neuron. We therefore
# transpose the matrix by using ".T":
ax[1].plot(v_mon.t/ms, v_mon.V.T/mV);
# Note that we put a ";" after the last line to avoid some output like
# "[<matplotlib.lines.Line2D at 0x7f83c5db33d0>]" in the notebook
```

### f/I curve

We will now plot an f/I curve as in previous week's tutorial.


```python
start_scope()
N = 100  # 100 neurons to have more values of I_stim
E_L = -50*mV
V_threshold = -30*mV
V_reset = -55*mV
C = 70*pF
g_L = 10*nS
tau = C/g_L
duration = 100*ms
neurons = NeuronGroup(N, '''dV/dt = (E_L - V + I_stim)/tau : volt
                            I_stim : volt (constant)  # note, this is rescaled with g_L already''',
                      threshold='V>V_threshold', reset='V=V_reset', method='euler')
# Initialize values
neurons.V = V_reset
neurons.I_stim = np.linspace(0, 50, N)*mV  # linspace(0, 50, N) gives us N equally spaced values between 0 and 50

# record membrane potential and spikes
v_mon = StateMonitor(neurons, 'V', record=True)
spike_mon = SpikeMonitor(neurons)

# run simulation
run(duration)
```


```python
# A SpikeMonitor has another attribute "count", which stores the total number of spikes for each neuron
spike_mon.count
```




    <spikemonitor_1.count: 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,  0,  0,  0,  0,  2,  4,  4,  5,  5,  6,  6,  7,  7,  7,  8,
            8,  8,  9,  9, 10, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13,
           14, 14, 14, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 18, 18, 18, 18,
           19, 19, 20, 20, 20, 20, 21, 21, 21, 22, 22, 22, 22, 23, 23],
          dtype=int32)>




```python
# Plotting the spike rate (number of spikes per second) against the input current for each neuron gives
# us an f/I curve:
plt.plot(neurons.I_stim/mV, spike_mon.count/duration);
```

Why do we get a "staircase" instead of a smooth curve? This could be because we are not using enough different values of input currents, e.g. imagine what would the curve look like if we only tested 10 different values. But it
turns out that this is not the case here. Let us simulate things again with 1000 neurons (i.e. 1000 different values of I_stim):


```python
start_scope()
N = 1000  # 1000 neurons to have more values of I_stim
E_L = -50*mV
V_threshold = -30*mV
V_reset = -55*mV
C = 70*pF
g_L = 10*nS
tau = C/g_L
duration = 100*ms
neurons = NeuronGroup(N, '''dV/dt = (E_L - V + I_stim)/tau : volt
                            I_stim : volt (constant)  # note, this is rescaled with g_L already''',
                      threshold='V>V_threshold', reset='V=V_reset', method='euler')
# Initialize values
neurons.V = V_reset
neurons.I_stim = np.linspace(0, 50, N)*mV  # linspace(0, 50, N) gives us N equally spaced values between 0 and 50

# record membrane potential and spikes
v_mon = StateMonitor(neurons, 'V', record=True)
spike_mon = SpikeMonitor(neurons)

# run simulation
run(duration)
```


```python
plt.plot(neurons.I_stim/mV, spike_mon.count/duration);
```

This looks even more like a staircase! If we look closely we see that the values jump up in steps of 10, i.e. there is a line at 100Hz, one at 110Hz, etc. This is because we simulate for only 100ms. To calculate the firing rate "per second", we multiply the number of spikes we get during 100ms with 10. Since there cannot be "half spikes", this means we will never see a firing rate like, say, 105Hz – either a neuron spikes 10 times and therefore the firing rate is 100Hz, or a neuron spikes 11 times and therefore the firing rate is 110Hz. The line gets smoother if we use a longer simulation time:


```python
start_scope()
N = 100  # back to 100 neurons
E_L = -50*mV
V_threshold = -30*mV
V_reset = -55*mV
C = 70*pF
g_L = 10*nS
tau = C/g_L
duration = 1000*ms  # 1s instead of 0.1s
neurons = NeuronGroup(N, '''dV/dt = (E_L - V + I_stim)/tau : volt
                            I_stim : volt (constant)  # note, this is rescaled with g_L already''',
                      threshold='V>V_threshold', reset='V=V_reset', method='euler')
# Initialize values
neurons.V = V_reset
neurons.I_stim = np.linspace(0, 50, N)*mV  # linspace(0, 50, N) gives us N equally spaced values between 0 and 50

# record membrane potential and spikes
v_mon = StateMonitor(neurons, 'V', record=True)
spike_mon = SpikeMonitor(neurons)

# run simulation
run(duration)
```


```python
plt.plot(neurons.I_stim/mV, spike_mon.count/duration);
```

This looks much better 👍

## Synapses

As mentioned in the lecture, synapses can be modeled on various models of detail. Here, we will use an "exponential current-based" synapse. Let us first try this out with some artficial stimulation:


```python
start_scope()
E_L = -50*mV
V_threshold = -30*mV
V_reset = -55*mV
C = 70*pF
g_L = 10*nS
tau = C/g_L
tau_syn = 2*ms
duration = 100*ms

# SpikeGeneratorGroup is an "artificial" neuron group. Instead of writing an equation for the membrane
# potential, etc. we just say: "emit spikes at certain times". This can be very useful for testing
# things. Here we define two groups with a single neuron (size 1). For each group we give a list of the
# neurons that should spike (here this is always "neuron 0", since we only have a single neuron in the
# group), and a list of the times when the neurons should spike. Note that this is the same convention
# as for the attributes "i" and "t" of a SpikeMonitor.
exc_stim = SpikeGeneratorGroup(1, [0, 0], [10, 85]*ms)
inh_stim = SpikeGeneratorGroup(1, [0, 0], [30, 80]*ms)

# Our neuron model is almost the same as previously, but now I_stim is replaced by I_syn. It is no
# longer a constant value, but instead defined by a differential equation: in between two spikes, the
# synaptic current I_syn decays to zero with time constant tau_syn
neurons = NeuronGroup(3, '''dV/dt = (E_L - V + I_syn)/tau : volt
                            dI_syn/dt = -I_syn/tau_syn : volt # decay between spikes''',
                      threshold='V>V_threshold', reset='V=V_reset', method='euler')
# Initialize values
neurons.V = E_L

# Synapses
# For each type of connection, we first describe the synaptic model. This consists of two parts: first
# we say which group of neurons are connected with each another (but not yet which specific neurons in
# these groups), and then what should happen if a spike arrives at the synapses. Here, excitatory
# synapses will increase the post-synaptic current I_syn by 1mV for each spike, and inhibitory neurons
# will decrease it by the same amount.
exc_syn = Synapses(exc_stim, neurons, on_pre='I_syn_post += 1*mV')  # "on_pre" = "for each pre-synaptic spike"
inh_syn = Synapses(inh_stim, neurons, on_pre='I_syn_post -= 1*mV')
# Now we need to define which neurons in the "source group" are connected to which neurons in the
# "target group". For very small networks like this, we can define each connection individually. For
# bigger networks, we would rather give some rule (connect all neurons to all neurons, connect neurons
# with a given probability, etc.). Here, we connect neuron 0→0, and 0→2 for the excitatory synapses, and
# 0→1 and 0→2 for the inhibitory synapses. This means that the first target neuron (neuron 0) gets only
# excitatory input, the second target neuron (neuron 1) gets only inhibitory input, and the third neuron
# (neuron 2) gets both.
exc_syn.connect(i=[0, 0],
                j=[0, 2])  # i = "source index", j = "target index"
inh_syn.connect(i=[0, 0],
                j=[1, 2])

# record membrane potential and spikes
v_mon = StateMonitor(neurons, 'V', record=True)
spike_mon = SpikeMonitor(neurons)

# run simulation
run(duration)
```


```python
# We plot the membrane potential of each target neuron in a separate line
fig, axs = plt.subplots(3, 1, sharex=True)
axs[0].plot(v_mon.t/ms, v_mon.V[0]/mV)
axs[1].plot(v_mon.t/ms, v_mon.V[1]/mV)
axs[2].plot(v_mon.t/ms, v_mon.V[2]/mV);
```

### Network

We now create a network of excitatory and inhibitory neurons, with connections among and across the groups. These kind of networks have been used to study network dynamics in the cortex, for example.


```python
start_scope()
N_E = 4000  # excitatory neurons
N_I = 1000  # inhibitory neurons
E_L = -50*mV
V_threshold = -30*mV
V_reset = -55*mV
C = 70*pF
g_L = 10*nS
tau = C/g_L
tau_syn = 2*ms
w_e = 1*mV
w_i = 8*mV
duration = 1000*ms

neurons = NeuronGroup(N_E + N_I,
                      '''dV/dt = (E_L - V + I_syn)/tau : volt
                         dI_syn/dt = -I_syn/tau_syn : volt # decay between spikes''',
                      threshold='V>V_threshold', reset='V=V_reset', method='euler')
# Initialize values
neurons.V = E_L

exc_neurons = neurons[:N_E]  # Uses "slicing" to get subpopulations of neurons
inh_neurons = neurons[N_E:]

# Connect neurons to each other
exc_syn = Synapses(exc_neurons, neurons, on_pre='I_syn_post += w_e')  # "on_pre" = "for each pre-synaptic spike"
inh_syn = Synapses(inh_neurons, neurons, on_pre='I_syn_post -= w_i')
# We now have 5000 neurons in our networks, so we can no longer type out all the connnections. Instead, we
# ask to make a randomly connected network.
exc_syn.connect(p=0.02)  # 2% connection probability between each pair of neurons  
inh_syn.connect(p=0.02)

# record spikes of all neurons, but the membrane potential only of 1 neuron
v_mon = StateMonitor(neurons, 'V', record=0)
spike_mon = SpikeMonitor(neurons)

# run simulation
# for simulations that take longer to run, it is useful to add "report='text'", which will give us some
# output during the simulation telling us how long the simulation has been run and how long it will
# still take.
run(duration, report='text')
```

    Starting simulation at t=0. s for a duration of 1. s
    1. s (100%) simulated in 1s



```python
fig, ax = plt.subplots(2, 1, sharex=True, figsize=(12, 6))
ax[0].plot(spike_mon.t/ms, spike_mon.i, ',')
ax[1].plot(v_mon.t/ms, v_mon.V.T/mV);
```

As we can see, the network isn't doing anything yet. This is because while the neurons are interconnected between each other (and could therefore spike due to the input of the other neurons), none of the neurons spike spontaneously and therefore no neuron propagates any activity to the other neurons. To get the activity going, we need to add some external stimulation, e.g. a current as in the previous models:


```python
start_scope()
N_E = 4000  # excitatory neurons
N_I = 1000  # inhibitory neurons
E_L = -50*mV
V_threshold = -30*mV
V_reset = -55*mV
C = 70*pF
g_L = 10*nS
tau = C/g_L
tau_syn = 2*ms
w_e = 1*mV
w_i = 8*mV
duration = 1000*ms

# The equation for V now includes I_ext, and I_ext is given as a parameter of the model
neurons = NeuronGroup(N_E + N_I,
                      '''dV/dt = (E_L - V + I_syn + I_ext)/tau : volt
                         I_ext : volt (constant)
                         dI_syn/dt = -I_syn/tau_syn : volt # decay between spikes''',
                      threshold='V>V_threshold', reset='V=V_reset', method='euler')
# Initialize values
neurons.V = E_L
neurons.I_ext = 21*mV  # stimulate the neurons enough so that they spike

exc_neurons = neurons[:N_E]  # Uses "slicing" to get subpopulations of neurons
inh_neurons = neurons[N_E:]

# Connect neurons to each other
exc_syn = Synapses(exc_neurons, neurons, on_pre='I_syn_post += w_e')  # "on_pre" = "for each pre-synaptic spike"
inh_syn = Synapses(inh_neurons, neurons, on_pre='I_syn_post -= w_i')
# We now have 5000 neurons in our networks, so we can no longer type out all the connnections. Instead, we
# ask to make a randomly connected network.
exc_syn.connect(p=0.02)  # 2% connection probability between each pair of neurons  
inh_syn.connect(p=0.02)

# record spikes of all neurons, but the membrane potential only of 1 neuron
v_mon = StateMonitor(neurons, 'V', record=0)
spike_mon = SpikeMonitor(neurons)

# run simulation
# for simulations that take longer to run, it is useful to add "report='text'", which will give us some
# output during the simulation telling us how long the simulation has been run and how long it will
# still take.
run(duration, report='text')
```

    Starting simulation at t=0. s for a duration of 1. s
    1. s (100%) simulated in 2s



```python
fig, ax = plt.subplots(2, 1, sharex=True, figsize=(12, 6))
ax[0].plot(spike_mon.t/ms, spike_mon.i, ',')
ax[1].plot(v_mon.t/ms, v_mon.V.T/mV);
```

Since all neurons get the same input current "I_ext", and all neurons start with the same value of "V", all neurons will emit a spike at the same time. This creates the vertical line at the start of the raster plot. This activity
then gets propagated to the other neurons. Since the connections are probabilistic, some neurons will now get a bit more excitatory or a bit more inhibitory input compared to the other neurons, and therefore the activity of the network starts to look "chaotic".

To get a closer look at the activity, we can only record from a subset of the neurons and plot their activity:


```python
start_scope()
N_E = 4000  # excitatory neurons
N_I = 1000  # inhibitory neurons
E_L = -50*mV
V_threshold = -30*mV
V_reset = -55*mV
C = 70*pF
g_L = 10*nS
tau = C/g_L
tau_syn = 2*ms
w_e = 1*mV
w_i = 8*mV
duration = 1000*ms

# The equation for V now includes I_ext, and I_ext is given as a parameter of the model
neurons = NeuronGroup(N_E + N_I,
                      '''dV/dt = (E_L - V + I_syn + I_ext)/tau : volt
                         I_ext : volt (constant)
                         dI_syn/dt = -I_syn/tau_syn : volt # decay between spikes''',
                      threshold='V>V_threshold', reset='V=V_reset', method='euler')
# Initialize values
neurons.V = E_L
neurons.I_ext = 21*mV  # stimulate the neurons enough so that they spike

exc_neurons = neurons[:N_E]  # Uses "slicing" to get subpopulations of neurons
inh_neurons = neurons[N_E:]

# Connect neurons to each other
exc_syn = Synapses(exc_neurons, neurons, on_pre='I_syn_post += w_e')  # "on_pre" = "for each pre-synaptic spike"
inh_syn = Synapses(inh_neurons, neurons, on_pre='I_syn_post -= w_i')
# We now have 5000 neurons in our networks, so we can no longer type out all the connnections. Instead, we
# ask to make a randomly connected network.
exc_syn.connect(p=0.02)  # 2% connection probability between each pair of neurons  
inh_syn.connect(p=0.02)

# record spikes of the first 100 neurons, but the membrane potential only of 1 neuron
v_mon = StateMonitor(neurons, 'V', record=0)
spike_mon = SpikeMonitor(neurons[:100])  # We can consider "neurons" like an array of neurons

# run simulation
# for simulations that take longer to run, it is useful to add "report='text'", which will give us some
# output during the simulation telling us how long the simulation has been run and how long it will
# still take.
run(duration, report='text')
```

    Starting simulation at t=0. s for a duration of 1. s
    1. s (100%) simulated in 2s



```python
fig, ax = plt.subplots(2, 1, sharex=True, figsize=(12, 6))
ax[0].plot(spike_mon.t/ms, spike_mon.i, '|')  # we often use a vertical bar for spikes
ax[1].plot(v_mon.t/ms, v_mon.V.T/mV);
```

89eb83dcacb50d20881fa28bde6cfaecac5f66c9
ipynb
Jupyter Notebook
tutorials/T08-Spiking-Neural-Networks.ipynb
mgraupe/DataSciPy2021
94454b2651e864531f0e2ce4a1691b92508fa3a2
CC-BY-4.0
2021-11-17T20:40:33.000Z
2021-12-31T11:21:49.000Z
tutorials/T08-Spiking-Neural-Networks.ipynb
mgraupe/DataSciPy2021
94454b2651e864531f0e2ce4a1691b92508fa3a2
CC-BY-4.0
tutorials/T08-Spiking-Neural-Networks.ipynb
mgraupe/DataSciPy2021
94454b2651e864531f0e2ce4a1691b92508fa3a2
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES