【Brian2学习教程之一】Introduction to Brian part 1: Neurons

前言：

目前正在学习脉冲神经网络建模，由于研究该领域的人员相对较少，所以很难在网上找到很多系统全面的中文学习教程，只能自己找英文原版资料进行学习。为了能够吸引更多脉冲神经网络建模领域的同学相互交流进步，于是决定开始开设Brian2学习教程系列专栏，总结自己所学，也分享知识给大家。由于水平有限，所以教程中可能会存在一些错误，欢迎大家批评指正。

本文内容根据Brian2官网提供的技术文档进行整理总结。

链接：https://brian2.readthedocs.io/en/stable/resources/tutorials/1-intro-to-brian-neurons.html

目录：

•  单位系统
• 一个简单模型
• 添加脉冲
• 不应期
• 多个神经元
• 参数
• 随机的神经元
• 结束教程

正文：

• 单位系统

在Brian中，存在一个专门的管理物理单位的系统。所有标准的国际单位制（A,K,s,m,kg,cd,mol）都可以在Brian中表示，同时这些单位的导出量及标准的单位前缀（V,mV,pF等）也可以在Brian中表示。下面列出几个例子，如表1所示：

表1 Brian单位系统实例

在Brian单位系统中，支持单位量之间进行数学运算，如表1中最后一个例子所示。但是，如果是不匹配的两个单位量之间进行无物理意义的运算，则会出现维度不匹配的错误。例如，将5A+10V在Brian中表示成5*amp+10*volt，则程序会弹出如下错误：

• 一个简单模型

接下来我们开始定义一个简单的神经元模型。在Brian中，所有的模型都是通过微分方程系统定义的。下面我们通过一个例子来详细说明。

tau = 10*ms
eqs = '''
dv/dt = (1-v)/tau : 1
'''

在Python中，'''符号用来开始和结束多行字符。所以上面这个例子中eqs其实只有一行等式dv/dt = (1-v)/tau : 1，它是用标准的数学符号表示的。在等式的最后我们可以看到: 1，它表示前面等式最终运算出来的单位，在Brian中，必须保证等式两边的单位保持一致。

现在让我们使用前面定义的微分方程创建一个神经元：G=NeuronGroup(1,eqs)。在Brian中，你只能使用类NeuronGroup创建一组神经元。在G=NeuronGroup(1,eqs)中存在两个参数，其中第一个参数1代表了这组神经元中包含的神经元的个数（在本例中，神经元的个数为1），第二个参数eqs代表了定义的微分方程。

接下来我们看一下，如果去掉等式中的变量tau会发生什么。执行如下代码。

from brian2 import *
eqs = '''
dv/dt = 1-v : 1
'''
G = NeuronGroup(1, eqs)
run(100*ms)

运行程序，可以看到发生如下错误：微分方程等式左右单位不一致。所以在Brian中，为了避免这种错误的发生，必须保证等式两边单位的一致性。

现在，我们将上面的等式修改正常，并在Spyder（Python3.6）中执行如下代码。

from brian2 import *
start_scope()

tau = 10*ms
eqs = '''
dv/dt = (1-v)/tau : 1
'''

G = NeuronGroup(1, eqs)
run(100*ms)

程序运行结果如下图所示，我们可以发现出现了一个“INFO”信息：没有指定数值积分方法。这对程序的执行没有什么影响，仅仅是提醒我们选择了哪种数值积分方法，我们可以通过显示指定的方法来避免出现这样一个“INFO”信息。

仅仅通过上面这段程序我们没法直观地看出程序执行完后发生了什么了？接下来我们将变量值v在仿真前后的结果打印出来。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallback
start_scope()

tau = 10*ms
eqs = '''
dv/dt = (1-v)/tau : 1
'''

G = NeuronGroup(1, eqs, method='exact')
print('Before v = %s' % G.v[0])
run(100*ms)
print('After v = %s' % G.v[0])

程序运行结果如下图所示，默认情况下，所有的变量初值都为0。但由于等式dv/dt=(1-v)/tau的作用，仿真执行一段时间后，变量v会趋近于1。这与程序运行结果一致。

现在我们绘制一个图来直观地看变量v是如何随着时间变化的。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallback
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])
xlabel('Time (ms)')
ylabel('v');

程序运行结果如下图所示。在这个例子中，我们使用了StateMonitor对象，它是用来记录仿真运行时神经元变量的值的。StateMonitor(G,'v',record=True)中有三个参数，第一个参数G表示记录哪一组神经元，第二个参数’v’表示想要记录哪个变量，第三个参数指定了记录方式。我们也可以指定record=0，这意味着我们将记录神经元0的所有值。一般情况下，我们必须指定我们想要记录哪一个神经元，因为在大规模仿真中存在很多神经元，如果记录所有神经元的变量值将会消耗大量的RAM。

接下来我们修改一下等式及一些参数看会发生什么变化。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallback
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');

程序运行结果如下图所示。

• 添加脉冲

前面一部分教程仅仅使用微分方程定义了一些简单的神经元模型，接下来我们让它们可以发放脉冲。执行如下代码：

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallback
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')

M = StateMonitor(G, 'v', record=0)
run(50*ms)
plot(M.t/ms, M.v[0])
xlabel('Time (ms)')
ylabel('v');

程序运行结果如下图所示。在这个例子中我们为NeuronGroup添加了两个新的关键词：threshold='v>0.8'和reset='v = 0'，这意味着当v>0.8时将会生成一个脉冲，脉冲过后v立即复位到0。

在上图中，不能直接看到脉冲的生成，但是在Brian内部已经以脉冲的形式记住了这样一个事件。SpikeMonitor对象将你想要记录的脉冲的那组神经元作为它的参数并将脉冲发放时间存储在变量t中。接下来我们将在图中直观地将这些脉冲展示出来。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallback
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');

程序运行结果如下图所示。在这个例子中，我们使用了matplotlib中的axvline命令绘制了一个橘黄色的垂直虚线来表示通过SpikeMonitor记录的不同时间点发放的脉冲。

• 不应期

对神经元模型来说最常见的一个特征是不应期。这意味着神经元在发放一个脉冲之后，它会在一段时间内保持不应期，在这段时间结束之前它都不能再发放另外一个脉冲。接下来我们将展示一下在Brian中如何实现它，执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallback
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*ms。这样仅仅意味着神经元在这段时间内不能放电，但它不会改变v的行为。为了使v在不应期内保持常量，在微分方程中v定义的末位添加了（unless refractory）。这意味着微分方程将决定v的行为除非在不应期这种作用将会失效。

接下来，我们将看一下如果没有添加（unless refractory）将会发生什么。要注意的是，为了使结果更清楚，我们减小了tau的值同时增加了不应期的长度。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallback
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');

程序运行结果如下图所示。从图中可以看出，第一个脉冲的行为和之前是一样的：v上升到0.8然后在8ms时刻神经元在立即复位到0之前发放一个脉冲。由于现在的不应期是15ms，这意味着神经元将在第8+15=23ms的时刻才会放电。在第一次放电之后，因为我们这次没有在dv/dt的定义中指定（unless refractory），所以v的值马上就会开始上升。然而，一旦它的值在大约持续8ms之后达到0.8（绿色虚线）甚至超过阈值v>0.8，它都不会发放一个脉冲。这是因为在23ms时刻之前，神经元一直处在不应期。

• 多个神经元

到目前为止我们只研究了一个神经元，接下使用多个神经元做一些有趣的事情。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallback
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=’rand()’将每个神经元的v指初始化为0到1之间均匀分布的随机值，用变量spikemon.i代表了对应神经元每个脉冲的发放时间。在图中用黑点表示，图中x轴表示时间，y轴表示神经元索引。这是神经科学领域中使用的标准的“raster plot”。

• 参数

为了使这多个神经元做一些更有趣的事情，我们为每个神经元引入一些与微分方程无关的参数。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallback
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为每个神经元申明了一个单位为1的参数v0。G.v0 = 'i*v0_max/(N-1)'为每个神经元初始化了一个从0增加到v0_max的v0，式中i代表了神经元索引。在这个例子中，对于v0<1的神经元，永远无法发放脉冲，但随着v0的逐渐增大，它将以一个更高的比率发放脉冲。右图展示了发放率作为v0的函数的结果，这是该神经元模型的I-f曲线。

• 随机的神经元

通常在制作神经元模型的时候，我们会包含一个随机项来模拟各种形式的神经元噪声。在Brian中，我们可以在微分方程中使用符号xi来做这件事情。严格地说，这个符号是一个“随机微分”但是你可以把它当作是一个均值为0标准差为1的高斯随机变量。考虑到单位一致性的问题，我们给它成了一个tau**-0.5。还有一点需要注意，我们这次使用‘euler’改变了method的关键字参数；我们之前使用的'exact'方法不能适用于随机微分方程。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallback
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)');

程序运行结果如下图所示。从图中可以看出，这个结果和上一节的结果类似，但是添加了一些噪声。请注意这个曲线是如何改变形状的：I-f曲线不再是从0一下子变成一个正值，而是以类似S形的趋势逐渐增加。这是由于无论多么小的驱动力，随机性都可能会导致神经元发放脉冲。

• 教程结尾

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

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallback
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=ones(len(spikemon))/(N*defaultclock.dt))
xlabel('Time (ms)')
ylabel('Instantaneous firing rate (sp/s)');

程序运行结果如下图所示。