[PYTHON] Simulate neural activity with Brian2

Research on simulating neural activity on a computer has been conducted for some time. NEURON of Yale University is famous as a simulator, but there are various other software.

This time, I will use a software called Brian2, which was published at eLife last year (2019).

• It's easy to make a mistake, but the name is Brian **, not ** Brain.

Home page: https://briansimulator.org/ GitHub：https://github.com/brian-team/brian2 Paper: https://elifesciences.org/articles/47314 The developer is Sorbonne University in France.

Installation

The environment is Ubuntu 16.04 Anaconda Python 3.6 is.

It seems that installation can be done with either conda or pip.

pip install brian2

You have now installed. The version was brian 2-2.3.0.2.

From the next section, I will try the basics while referring to Document Tutorial.

Try to move

First, let's do a simple membrane potential simulation.

Import. (In this section, to make it easier to understand which module is, I dare to import it normally instead of from brian2 import * as per the manual.)

import brian2

** Brian2 has a module as a physical unit, and you can define Hz, ms, etc. as they are. ** ** You can use it to describe parameters and constants concisely. In addition, it is possible to prevent setting mistakes because an error will occur if the physical units are not aligned during calculation.

v0 =  -10 * brian2.mV
tau = 10 * brian2.ms

To start the simulation, first insert start_scope ().

brian2.start_scope()

The differential equation (which defines the change in membrane potential) seems to be defined as a String type.

eqs = 'dv/dt = (v0-v)/tau : volt (unless refractory)'

The above volt is the physical unit of the variable. Expressed in SI units. Here, the membrane potential is a variable. (unless refractory) is required when using the refractory period.

Nerve cells are defined in the Neuron Group module.

G = brian2.NeuronGroup(1, eqs, threshold='v > -40*mV', reset='v = -60*mV', refractory=5*brian2.ms, method='exact')

The first argument is the number of nerve cells. This time I set it with one. The second argument is the differential equation of the membrane potential. The other arguments are threshold is the threshold potential that becomes a spike reset is the reset potential after spike refractory is refractory method is a numerical calculation method Is set.

Set the initial value of the membrane potential.

G.v = -70 * brian2.mV

You can monitor with StateMonitor.

M = brian2.StateMonitor(G, 'v', record=True)

The simulation is started by run (<time>).

brian2.run(50*brian2.ms)

You can now simulate.

Plot the results.

brian2.plot(M.t/brian2.ms, M.v[0]/brian2.mV)
brian2.xlabel('Time (ms)')
brian2.ylabel('Ptential (mV)');

This result is a simulation of neuronal spikes in the LIF model. (Although there is no peak potential.)

Now you know the flow of the simulation.

As an evolution, it seems that it is also possible to handle stochastic differential equations containing noise by using xi (random variable of standard normal distribution).

Spike simulation of multiple nerve cells

Next, we simulate the spikes of multiple nerve cells. This time, I will give each cell a different input potential (I).

I will start from import.

from brian2 import *

start_scope()

N = 100 #Number of cells
tau = 10*ms #Membrane time constant
v0 = -50*mV #Resting membrane potential

eqs = '''
dv/dt = (v0-v+I)/tau : volt (unless refractory)
I : volt
'''

G = NeuronGroup(N, eqs, threshold='v>-40*mV', reset='v=-60*mV', refractory=1*ms, method='exact')
G.v = -70*mV

Use SpikeMonitor if you want to see only spikes, not changes in membrane potential.

M = SpikeMonitor(G)

Define a different input for each cell.

I_max = 40*mV
G.I = 'i*I_max/(N-1)'

i is the cell index. You have now defined different inputs.

run

duration = 1000*ms
run(duration)

Create a raster plot.

plot(M.t/ms, M.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index');

As the cell index goes up, the input potential goes up, so there are more spikes.

I will also make an F-I curve that you often see.

plot(G.I/mV, M.count/duration)
xlabel('Input (mV)')
ylabel('Firing rate (sp/s)');

it is complete. This is said to be similar to the ReLU function.

To simulate a neural circuit

I will omit the details, but it seems that you can easily simulate neural circuits.

Synapses can be defined in the Synapses module. The connect method connects the cells. (Which cell to connect to which cell)

It is possible to define a differential equation for the strength of synapses, and it seems that an STDP model can also be created.

Besides, it seems that you can define the network with the Network module.

Impressions

I thought it was easy to use because it had a lot of functions for simulation.

There are many other examples in the documentation. https://brian2.readthedocs.io/en/stable/examples/index.html There was also an example of a compartment model. I would like to try these as well.

• Addition In January of this year, Brian2GeNN was announced. https://www.nature.com/articles/s41598-019-54957-7 It seems that GPU can now be used as well.