Brian2学习笔记_Introduction to Brian part1:Neurons

本文根据Brian2官方英文教程进行翻译总结。
Brian2官方安装教程链接：
https://brian2.readthedocs.io/en/stable/introduction/install.html
使用教程Tutorial：
Introduction to Brian part 1 : Neurons 原文
Introduction to Brian part 2 : Synapses 原文
Introduction to Brian part 3 : Simulation 原文

运行的三种环境：Brian提供的基于网页的交互式模拟仿真器“Binder”；本地的Jupyter Notebook；Ipython或标准Python终端。

Introduction to Brian part 1: Neurons

静态不可编辑版本
所有Brian2脚本从以下内容开始：

from brian2 import *
%matplotlib inline

1.1 单位制

可以使用所有基本的SI单位，所有标准前缀，以及一些特殊缩写。

单位	符号	单位	符号	单位	符号
伏特	volt	秒	second(s)	米	meter
安培	amp	赫兹	hertz(Hz)	升	liter
欧姆	ohm	焦耳	joule	克	gram
西门子	siemens(S)	瓦特	watt	千克	kilogram
法拉	farad(F)	帕	pascal	摩尔	mol
库伦	coulomb	开尔文	kelvin	坎德拉	candela

标准前缀及部分缩写

p-	n-	m-	k-	M-	G-	T-
pico	nano	milli	kilo	million	giga	tera
mV	cm	pF

给数据加上单位：用 * 表示；
给数据去除单位：1.用 / 除去单位。2.用numpy的数列函数’asarray()‘或’array()’。
标准格式：

输入	20*volt	1000*amp	1e6*volt	10*nA*5*Mohm
输出	20.0V	1.0kA	1.0MV	50.0mV

rates = [10, 20, 30]*Hz
print(rates/Hz)
print(array(rates))

[10. 20. 30.] #方法1输出
[10. 20. 30.] #方法2输出

1.2 定义神经元

在Brian中，所有模型都是由微分方程系统定义的。下面是一个简单的例子：

start_scope()
tau = 10*ms
eqs = '''
      dv/dt = (1-v)/tau : 1
      '''
G = NeuronGroup(1, eqs，method='exact')

$start\_scope()$用于清空之前的变量。
Python中，$\prime \prime \prime$用于多行字符串的开始和结束。所以方程就是一个字符串，每个方程一行。在一行的末尾，需要加上$:unit$。unit是变量的SI单位。这不是方程两边的单位(1/s)，而是方程定义的变量的单位。在本式中，即$v$。注意：在Brian2中，方程两边的单位必须保持一致。
在Brian中，只需要创建神经元群组。类NeuronGroup用于创建神经元。前两个参数是神经元的数量（本例中为1）和定义的微分方程。
$method=’exact‘$表求解此方程采用数值积分法。
可以通过如下方法观察程序执行前后发生了什么：

print('Before v = %s' %G.v[0])
run(100*ms)
print('After v = %s' %G.v[0])

运行结果如下所示：

Before v = 0.0
After v = 0.9999546000702376

默认情况下，所有变量初始值均为0。在程序运行100ms后，变量$v$会趋近于1。运行结果与预期一致。
还可以绘图来直观地展示$v$如何随着时间而变化。运行以下程序，画出预期行为与程序结果的对比图：

start_scope()
tau = 10*ms
eqs = '''
dv/dt = (1-v)/tau : 1
'''
G = NeuronGroup(1, eqs, method='exact')
M = StateMonitor(G, 'v', record=True)
run(30*ms)
plot(M.t/ms, M.v[0], 'C0', label='Brian')
plot(M.t/ms, 1-exp(M.t/tau), 'C1--', label='Analytic')
xlabel('Time (ms)')
ylabel('v')
legend

运行结果如下：可以发现预期（橙色虚线）与仿真（蓝色实线）重合

在本例中，使用了StateMonitor对象，用于记录仿真运行时神经元变量的值。$StateMonitor(G{,}^{\prime }{v}^{\prime },record=Ture)$中有三个参数。第一个参数G代表记录的神经元组别；第二个参数$v$表示想要记录的变量；第三个参数指定记录方式。也可以指定 $record=0$ ，表示记录神经元的所有值。一般情况下，必须指定想要记录哪一个神经元，因为如果记录所有神经元的变量值将会消耗大量的RAM。
现在改变一些参数：

start_scope()
tau = 10*ms
eqs = '''
dv/dt = (sin(2*pi*100*Hz*t)-v)/tau : 1
'''
# Change to Euler method because exact integrator doesn't work here
G = NeuronGroup(1, eqs, method='euler')
M = StateMonitor(G, 'v', record=0)
G.v = 5 # initial value
run(60*ms)
plot(M.t/ms, M.v[0])
xlabel('Time (ms)')
ylabel('v')

1.3 添加脉冲

接下来使添加的神经元可以发射脉冲，并标注出每次发射的脉冲：

start_scope()
tau = 10*ms
eqs = '''
dv/dt = (1-v)/tau : 1
'''
G = NeuronGroup(1, eqs, threshold='v>0.8', reset='v = 0', method='exact')
statemon = StateMonitor(G, 'v', record=0)
spikemon = SpikeMonitor(G)
run(50*ms)
plot(statemon.t/ms, statemon.v[0])
for t in spikemon.t:
    axvline(t/ms, ls='--', c='C1', lw=3)
xlabel('Time (ms)')
ylabel('v')

运行结果如下：

本例中，我们为NeuronGroup添加了两个新的关键词：$threshold{=}^{\prime }v>0.8^{\prime }$ 和 $reset{=}^{\prime }v=0^{\prime }$ ，这表示当 $v>0.8$ 时会生成一个脉冲，脉冲过后 $v$ 立即复位到0。
上图中，$v$ 的变化开始时与之前一样，直到 $thresholdv>0.8$ ，$v$ 复位到 $0$，生成一个脉冲。
SpikeMonitor对象将我们想要记录的脉冲的那组神经元作为它的参数并将脉冲触发时间存储在变量 $t$ 中。
在这里，我们使用了matplotlib中的axvline命令，在SpikeMonitor记录的每个脉冲发生时绘制了一条橙色的虚线。

1.4 不应期（Refractoriness）

神经元模型是的一个最常见的特征就是不应期。也就是神经元在发放一个脉冲之后，它会在一段时间内保持不响应状态，在这段时间结束之前它都不能再发射另外一个脉冲。执行如下代码：

from brian2 import *
start_scope()
tau = 10*ms
eqs = '''
dv/dt = (1-v)/tau : 1 (unless refractory)
'''
G = NeuronGroup(1, eqs, threshold='v>0.8', reset='v = 0', refractory=5*ms, method='exact')
statemon = StateMonitor(G, 'v', record=0)
spikemon = SpikeMonitor(G)
run(50*ms)
plot(statemon.t/ms, statemon.v[0])
for t in spikemon.t:
    axvline(t/ms, ls='--', c='C1', lw=3)
xlabel('Time (ms)')
ylabel('v')

结果如下：

可以看出，在第一个脉冲之后，$v$ 在恢复到正常行为之前持续 $5ms$ 保持为 $0$。为了做到这一点，我们需要执行两个操作。首先，在NeuronGroup中添加关键词 $refractory=5\ast ms$。这样仅意味着神经元在这段时间内不会发射脉冲，但它不会改变 $v$ 的行为。为了使 $v$ 在不应期内保持恒定，我们必须在微分方程中对 $v$ 的定义后面添加 (unless refractory)。这意味着微分方程决定了 $v$ 的行为，但在不应期时它的作用会失效。
以下是不添加 (unless refractory) 时会发生的情况。在这里，为了更清晰地展示，我们降低了 $tau$ 的值并增加了不应期长度。运行如下代码：

start_scope()
tau = 5*ms
eqs = '''
dv/dt = (1-v)/tau : 1
'''
G = NeuronGroup(1, eqs, threshold='v>0.8', reset='v = 0', refractory=15*ms, method='exact')
statemon = StateMonitor(G, 'v', record=0)
spikemon = SpikeMonitor(G)
run(50*ms)
plot(statemon.t/ms, statemon.v[0])
for t in spikemon.t:
    axvline(t/ms, ls='--', c='C1', lw=3)
axhline(0.8, ls=':', c='C2', lw=3)
xlabel('Time (ms)')
ylabel('v')
print("Spike times: %s" % spikemon.t[:])

结果如下：

可以看出，第一个脉冲的行为与之前是相同的：$v$ 在第 $8ms$ 上升到 $0.8$ 然后立即复位到 $0$ 并发射一个脉冲。由于现在的不应期是 $15ms$，这意味着神经元在 $8+15=23ms$ 后才可以继续放电（发射脉冲）。在第一个脉冲发射后，由于这次我们在 $dv/dt$ 的定义中没有指定 (unless refractory)，所以 $v$ 的值立即开始上升。然而，一旦它的值达到 $0.8$，即使阈值是 $v>0.8$，也不会发射脉冲。这是因为在 $23ms$ 时刻前，神经元一直处在不应期，不会放电。有关不应期工作原理的更多详细信息，请参阅完整文档Refractoriness

1.5 多个神经元

到目前为止，我们只研究了一个神经元，接下使用多个神经元：

start_scope()
N = 100
tau = 10*ms
eqs = '''
dv/dt = (2-v)/tau : 1
'''
G = NeuronGroup(N, eqs, threshold='v>1', reset='v=0', method='exact')
G.v = 'rand()'
spikemon = SpikeMonitor(G)
run(50*ms)
plot(spikemon.t/ms, spikemon.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index')

运行结果：

在上述代码中，相较之前有几个变化。首先，我们使用 $N$ 代表神经元的个数。其次，用 $G.v{=}^{\prime }rand({)}^{\prime }$ 将每个神经元的 $v$ 值初始化为 $0$ 到 $1$ 之间均匀分布的随机值，这样是为了让每一个神经元有一定的差异。
当然，最大的差异还是我们的数据展现形式。用变量 spikemon.t 代表了对应神经元 $i$ 的每个脉冲的发放时间，spikemon.i 代表对应的第 $i$ 个神经元，在图中用黑点表示。其中x轴表示时间，y轴表示神经元索引。这是神经科学领域中经常使用的标准“栅格图”(raster plot)。

1.6 参数

为了让这些神经元做更多的事情，可以为每个神经元引入一些参数。这些参数与微分方程无关。运行如下代码：

start_scope()
N = 100
tau = 10*ms
v0_max = 3.
duration = 1000*ms
eqs = '''
dv/dt = (v0-v)/tau : 1 (unless refractory)
v0 : 1
'''
G = NeuronGroup(N, eqs, threshold='v>1', reset='v=0', refractory=5*ms, method='exact')
M = SpikeMonitor(G)
G.v0 = 'i*v0_max/(N-1)'
run(duration)
figure(figsize=(12,4))
subplot(121)
plot(M.t/ms, M.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index')
subplot(122)
plot(G.v0, M.count/duration)
xlabel('v0')
ylabel('Firing rate (sp/s)')

运行结果如下：

$v0:1$ 为每个神经元声明了一个新参数 $v0$ 。单位为 $1$，即无量纲。
$G.v0{=}^{\prime }i\ast v0\_max/(N-1{)}^{\prime }$ 将每个神经元的 $v0$ 初始化为一个从 $0$ 增加到 $v0\_max$ 的序列。
在这个例子中，我们以指数的方式驱使神经元电压值向 $v0$ 靠近，但是当 $v$ 越过 $1$ 时，它会触发神经元发射一个脉冲并复位。其效果是，它触发峰值的速率将与 $v0$ 的值相关。当 $v0<1$ 时，它将永远不会发射脉冲，当 $v0$ 变大时，它将以更高的速率发射脉冲。右边的图显示了发射速率(firing rate)作为 $v0$ 的函数的结果。这是神经元模型的 $I-f$ 曲线。
在本例中，我们使用了SpikeMonitor的count变量。这是一组神经元中每个神经元发射的脉冲数目的数组。用它除以运行时间duration就得到了发射速率(firing rate)。

1.7 随机神经元

通常在神经元模型中会包含一个随机参数来模拟各种形式的神经元噪声。在Brian2中，我们可以在微分方程中引入符号 $xi$ 来做这件事情。严格地说，这个符号是一个“随机微分”，但是可以把它当作是一个均值为 $0$ ，标准差为 $1$ 的高斯随机变量。我们还应该考虑到单位一致性的问题，所有我们在下面的代码中引入了一个 $tau\ast \ast -0.5$ 。还有一点需要注意，我们这次使用欧拉法（‘euler’）来改变method的关键字参数，而在前面使用的’exact’方法不适用于随机微分方程。

start_scope()
N = 100
tau = 10*ms
v0_max = 3.
duration = 1000*ms
sigma = 0.2
eqs = '''
dv/dt = (v0-v)/tau+sigma*xi*tau**-0.5 : 1 (unless refractory)
v0 : 1
'''
G = NeuronGroup(N, eqs, threshold='v>1', reset='v=0', refractory=5*ms, method='euler')
M = SpikeMonitor(G)
G.v0 = 'i*v0_max/(N-1)'
run(duration)
figure(figsize=(12,4))
subplot(121)
plot(M.t/ms, M.i, '.k')
xlabel('Time (ms)')
ylabel('Neuron index')
subplot(122)
plot(G.v0, M.count/duration)
xlabel('v0')
ylabel('Firing rate (sp/s)');

这与前一节中的结果类似，但增加了一些噪声干扰。注意曲线是如何改变形状的：它现在以s型曲线的方式增加，而不是从0到正值的跳跃。这是因为无论驱动力有多小，随机性都可能导致神经元发射一个脉冲。

1.8 End of tutorial 结束教程

最后，给出一个例子，看你能否明白它是在做什么。
尝试添加一个 StateMonitor 记录其中一个神经元的变量值来帮助理解这个例子。
也可以在这个例子中尝试一下在本教程中学到的其他内容。如果你完成了，就可以进入下一个关于突触的教程了。

start_scope()

N = 1000
tau = 10*ms
vr = -70*mV
vt0 = -50*mV
delta_vt0 = 5*mV
tau_t = 100*ms
sigma = 0.5*(vt0-vr)
v_drive = 2*(vt0-vr)
duration = 100*ms

eqs = '''
dv/dt = (v_drive+vr-v)/tau + sigma*xi*tau**-0.5 : volt
dvt/dt = (vt0-vt)/tau_t : volt
'''

reset = '''
v = vr
vt += delta_vt0
'''

G = NeuronGroup(N, eqs, threshold='v>vt', reset=reset, refractory=5*ms, method='euler')
spikemon = SpikeMonitor(G)

G.v = 'rand()*(vt0-vr)+vr'
G.vt = vt0

run(duration)

_ = hist(spikemon.t/ms, 100, histtype='stepfilled', facecolor='k', weights=list(ones(len(spikemon))/(N*defaultclock.dt)))
xlabel('Time (ms)')
ylabel('Instantaneous firing rate (sp/s)');

下一节：Synapses突触