Brian2学习笔记

前言

运行环境

系统：Ubuntu 16.04
架构：x86-64
核数：8核
内存：16G

写点有用的

原文链接：http://brian2.readthedocs.io/en/stable/introduction/index.html

从2018年初至现在使用Brian进行脉冲神经网络(spiking neural network，SNN)模拟一年时间，主要研究内容是将SNN运用到图像识别中，曾将Dielh的论文¹开源代码从Brian转移到Brian2，然后基于Brian2进行了一些SNN学习算法的探索性实验，主要的使用体验如下：

• 神经元和突触的动态特性定义比较简单，因为公式和参数都是可查阅到的，基于Brian的SNN模拟实验设置主要内容在于所设计的连接的建立，所以对于突触连接的几种定义方式要格外熟悉，灵活运用。
• 第二点想说的是C++ Code Generation模式对于拥有复杂模型的网络加速明显，但我也遇到很多奇怪的情况，实验经常跑崩掉，网络规模越大跑崩的几率越大（runtime generation模式倒是没有出现跑崩的情况），而且发生地没有规律，主要问题应该出在STDP的加速模块上。
• Brian这个模拟器用起来倒不难，而且如果细心着灵活着用可以发现很多可以优化的细节，收益也是可观的。
• 最后，我觉得Brian主要有两个不足：第一是缺乏运用多核甚至计算集群的能力，据说并行模拟做的较好的是NEST，它是面向计算神经科学的，模型复杂度低，规模可以很大，重视对网络信息处理、网络动态特性以及网络的学习和可塑性的探索。第二，每次模拟编译器都需要对模型编译，尤其对于大规模网络的模拟是不利的。虽然可以运用standalone code generation mode将python脚本转换为C++代码，但是在C++基础上修改和使用是不容易的，尤其对C++不熟悉的同志。

没用的

初稿写于2018年6月，最近一看，写的什么东西啊，怎么就发出去了！对不起了！我重写！
本篇博客主要内容是翻译自官方文档，英语好的还是推荐看源文。我这个大致浏览一下，可以了解个大概。
另外，关于Brian和Brian2区别的翻译文章在另一篇博客Brian到Brian2的转换

简介

Brian是一款面向SNN的开源模拟器。编程语言为Python，具有简单易学、灵活和可扩展的功能，旨在降低研究人员的编程时间和计算机的模拟时间。2018年十月发布最新版本为 Brian 2.2。
下面一段代码定义了一个随机连接的神经元网络，神经元模型为IF（integrate and fire），神经元数量为4000，激活和抑制的电流为指数型，模拟整个网络的运行1秒钟，并绘制神经元的fire/spike记录。

from brian2 import *
eqs = '''
dv/dt  = (ge+gi-(v+49*mV))/(20*ms) : volt
dge/dt = -ge/(5*ms)                : volt
dgi/dt = -gi/(10*ms)               : volt
'''
P = NeuronGroup(4000, eqs, threshold='v>-50*mV', reset='v=-60*mV')
P.v = -60*mV
Pe = P[:3200]
Pi = P[3200:]
Ce = Synapses(Pe, P, on_pre='ge+=1.62*mV')
Ce.connect(p=0.02)
Ci = Synapses(Pi, P, on_pre='gi-=9*mV')
Ci.connect(p=0.02)
M = SpikeMonitor(P)
run(1*second) 
plot(M.t/ms, M.i, '.')
show()

除了在本地安装使用，Brian也提供了基于网页的联网交互模式（Interactive demo），挂载在mybinder.org。

作者推荐的使用和学习流程是这样：安装；学习tutorial；运行例程；学习User guide，里面有每个模块和方法的详细参考链接以及底端的代码；尝试自行编写或者修改模拟文件；学习完整的Documents。如果遇到问题可以通过email support list联系作者，然后Brian也放在了Github上管理了brian2，可以在Issues上提问和提建议。

Brian团队还开发了几个辅助功能的独立工具：

• brian2tools：用于Brian结果的绘图和分析
• Brian Hears：建模听觉系统
• Model fitting toobox：自动根据电生理数据调整模型

引用

[1] Goodman DF and Brette R (2009). The Brian simulator. Front Neurosci doi:10.3389/neuro.01.026.2009
[2] Stimberg M, Goodman DFM, Benichoux V, Brette R (2014). Equation-oriented specification of neural models for simulations. Frontiers Neuroinf, doi: 10.3389/fninf.2014.00006.

安装

三种安装方法：作者推荐的Anaconda（Anacoda是一种流行的Python数据科学的开发平台）， Anaconda发行的版本自带必须的依赖包，不需要使用源文件在本地编译，且自动的连续性测试的工具traivs和appveyor都是基于Anaconda，作者对于Anaconda的发行版本更有信心。也有其他一些Python开发的平台自带了Brian，例如Enthought Canopy, Python(x,y) for Windows, …，当然最直接的就是在纯净的Python中编译或者Python包管理工具pip安装。
在此仅介绍我使用的方法：

• 从Github存档的源编译安装
• pip安装
  需要注意的是这两种安装方法最好不要在系统的Python中进行，因为Brian安装的包可能和系统所需的冲突，我是付出惨痛代价的！最好的方法是使用Pipenv，Virtualenv，pyenv等工具为Brian建立独立的虚拟环境，这个可以参考我的博客我的pipenv工作流程和方法以及另一篇virtualenv的文章

python编译安装

安装步骤

• 从Github下载源码到本地（Github上还存档了Brian的最终版Brian1.4.4以及你可能感兴趣的工具包，可以在brian-team中浏览下载安装）
• 解压
• 激活所建立的Python虚拟环境
• 用pip安装Brian2依赖的包：numpy（科学计算），matplotlib（绘图），cython， scipy（科学计算），nose（单元测试）。
• 进入解压后的含有setup.py的主目录，终端运行“python setup.py install”，安装就完成了。
• 终端中运行“pip list”查看所安装的Brian和其他包的版本，有的包的版本如果太新而不兼容可以通过“pip install packagename==x.x”来指定版本安装。

pip安装

pip安装比从源文件要简单，安装步骤如下：

• 激活所建立的Python虚拟环境
• 用pip安装Brian2依赖的包：numpy（科学计算），matplotlib（绘图），cython， scipy（科学计算），norse（单元测试）。
• 终端运行“pip install brian2”，等待安装完成。此种方式安装的版本会略低于源码安装。
• 终端中运行“pip list”查看所安装的Brian和其他包的版本，有的包的版本如果太新而不兼容可以通过“pip install packagename==x.x”来指定版本安装。

C++ code generation的安装要求

作者建议在模拟中使用C++ code generation模式，会极大地提升模拟速度（但也会出错或者不稳定），如果要使用C++ code generation模式必须要有C++的编译器（例如linux的g++）Cython或者weave（weave只可以在Python2.x中使用），使用Anaconda方式安装会自动安装Cython或者weave，如果Python的模拟环境中没有可以通过pip安装。Windows在C++编译器的选择和安装上会比Unix系统复杂，而且也经常碰到问题，所以建议在Unix系统上进行实验。

测试

brian2的测试是基于nose而不是基于python自带的unittest，在测试之前需要安装nose包：

• 激活所需的虚拟环境
• 在终端中运行“python”
• 在python中运行“import brian2”和“brian2.test()”来测试，测试成功了输出的底部是“OK”而且只有success和skip。

使用教程 Tutorial

练习Tutorial有三种环境：Brian提供的基于网页的交互式模拟仿真器“Binder”； 本地Jupyter notebook（inline plotting需要加%matplotlib inline）；Ipython或者标准的python终端。

part 1: 神经元 Neurons

物理单位

单位	符号	单位	符号	单位	符号
伏	volt	秒	second	米	meter
安培	amp	凯尔文	kelvin	坎德拉	candela
欧姆	ohm	瓦特	watt	西蒙子	siemens
千克	kilogram	帕斯卡	pascal	升	liter
焦耳	joule	赫兹	hertz	克	gram
摩尔 6.022 e 23 6.022e^{23} 6.022e23	mole/mol	库伦	coulomb	法拉	farad

数量级除了用科学计数法表示也可以可以用前缀标明p, n, u, m, k, M, G, T，例如毫伏mV、纳安namp、兆欧Mohm；两个特例kilogram不使用任何前缀；metre and meter多一种“centi”前缀表示厘米cmetre/cmeter。还其他便利的缩写：cm ( cmetre/cmeter), nS (nsiemens), ms (msecond), Hz (hertz), mM (mmolar)，pF(picofarad)。

给一个数据加上单位的方法是‘乘以’一个单位，例如‘3*ms’。去除单位变为纯数字的方法有三个，最常用的是‘除以’，例如‘t/ms’得到的结果就是纯数字；第二种方法是用numpy的数列函数‘asarray()’或者‘array()’，例如‘asarray(speed)’；第三种方法是在变量后加‘_’，例如：

print(M.v_[0])

需要特别说明的一点是Brian2对单位的要求更加严格规范，而且也支持数组的单位。更多关于单位。

定义神经元

神经元的特性包括：胞膜电压动态特性用微分方程组来表示，阈值threshold（threshold的条件一般是胞膜电压‘v’的不等式，也可以是包括时间‘t’在内的其他自定义变量的不等式），复位reset（复位操作可以对神经元的变量进行操作），不应期refractory。以下为一个简易的定义神经元（1个神经元）的例子：

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

神经元的动态特性是由一元微分方程来定义的，方程最后要标明变量的SI单位，注意不是方程运算结果的单位，方程左右两边单位的运算结果是相同的。
下面这个例子定义了一个神经元，并且定义了StateMonitor来记录神经元的状态变量，可以在模拟完成后进行绘图可视化，观察神经元的模型运行是否是准确。对于StateMonitor例子中定义了“record=0”表示记录所有神经元的变量‘v’，在大规模模拟中可以定制需要记录的变量，否则消耗大量内存。

start_scope()
G = NeuronGroup(1, eqs, method='exact')
M = StateMonitor(G, 'v', record=0)
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();

输出的绘图如下所示：

下面这段代码是指定几个神经元的某些变量进行记录：

G = NeuronGroup(...)
M = StateMonitor(G, ('v', 'u'), record=[0, 10, 100])
run(...)
plot(M.t/ms, M.v[0]/mV, label='v')
plot(M.t/ms, M.u[0]/mV, label='u')

start_scope()作用是用在jupyter notebook中，保证在这个函数之前定义的对象不纳入下一次的模拟中。
更多关于神经元、方程定义、方程运算以及变量记录（Brian2的变量监测记录的内容包括神经元和突触的内变量，spike的时间和索引序列、计数，spike的频率）

生成 spike

现在考虑让神经元产生Spike并记录产生Spike的时刻（spike在模拟器中是一个事件），要产生spike就需要加入threshold、reset的定义，用SpikeMonitor来检测和记录（SpikeMonitor还可以记录spike时刻神经元的内在变量）。

start_scope()
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');

输出的绘图如下所示：

不应期 refractory

神经元有一个普遍的特性是不应期，即在神经元spike之后会在一段时间进入一个静息状态，在这段时间内不接受外界刺激（状态变量维持在静息的水平）也不能再次spike。refractory不是用公式定义，而直接用数字表示。神经元定义中的refractory是让神经元在5ms之内不再spike，神经元模型中的 “(unless refractory)”则使得神经元的包膜电压在那段时间不受微分方程的控制，也就不变了。

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')

更多关于refractory

参数 parameters

实际的模拟中你需要为神经元定义很多参数来做一些有意思的模拟，这些参数可以是per-neuron的。这些变量通常在神经元的模型中定义，而后可以为每个神经元的参数赋予不同的值。如下：

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)

'v0’就是每个神经元的一个定制参数，这样每个神经元就有不同的动态参数了，这个参数的赋值通常是随机的或者根据索引（内部定义的变量‘i’）有规律地赋值，从而在模拟中探索规律。再举一个有意思的例子就是给神经元加入坐标变量，那样就可以根据神经元坐标关系定制模拟实验。

随机神经元 Stochastic neuron

Brian2中内置了一个随机变量‘xi’(符合高斯随机分布平均值为0标准差为1)，可以构造随机差分方程（此时差分方程求解方法是‘euler’）。

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)'

part 2: 突触 Synapses

神经元通过突触连接，前级神经元通过突触将spike传递给后级神经元，并且可以定义操作，改变神经元和突触的状态量，这是神经元网络功能的基础。
如下为一个简单的例子，“on_pre=‘v_post += 0.2’”的意义为接收到前级(_pre)的spike后后级神经元的包膜电压(_post)增加0.2V，更加复杂的计算均可以通过on_pre的公式定制。需要注意的是synapse的定义和连接是分两步进行的。

eqs = '''
dv/dt = (I-v)/tau : 1
I : 1
tau : second
'''
G = NeuronGroup(2, eqs, threshold='v>1', reset='v = 0', method='exact')
G.I = [2, 0]
G.tau = [10, 100]*ms

# Comment these two lines out to see what happens without Synapses
S = Synapses(G, G, on_pre='v_post += 0.2')
S.connect(i=0, j=1)

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

run(100*ms)

plot(M.t/ms, M.v[0], label='Neuron 0')
plot(M.t/ms, M.v[1], label='Neuron 1')
xlabel('Time (ms)')
ylabel('v')
legend();

输出结果如下所示：

更加实际的情况是发生spike后改变突触的权值，Hebbian学习法则、shor term plasticity以及STDP都是以此为基础的。如下例子：

eqs = '''
dv/dt = (I-v)/tau : 1
I : 1
tau : second
'''
G = NeuronGroup(3, eqs, threshold='v>1', reset='v = 0', method='exact')
G.I = [2, 0, 0]
G.tau = [10, 100, 100]*ms
# Comment these two lines out to see what happens without Synapses
S = Synapses(G, G, 'w : 1', on_pre='v_post += w')
S.connect(i=0, j=[1, 2])
S.w = 'j*0.2' # add a weight respectively
S.delay = 'j*2*ms' # introducing a delay
M = StateMonitor(G, 'v', record=True)
run(50*ms)
plot(M.t/ms, M.v[0], label='Neuron 0')
plot(M.t/ms, M.v[1], label='Neuron 1')
plot(M.t/ms, M.v[2], label='Neuron 2')
xlabel('Time (ms)')
ylabel('v')
legend();

注意这里的权值是电压型的(而且是瞬时的，instantaneous)，也就是直接增加胞膜电压，突触模型也有电流型和电导型的突触，最广泛的是电导型的。connect函数中’i’默认是前级神经元的索引， 'j’是后级神经元的索引。定义突触的连接方式是SNN模拟中最关键的部分，其他更加灵活的连接定义方式可以查阅用户手册的突触

延时 introducing a delay

在connect后加如下语句即可，也可以定义随机的延时或者其他基于索引或突触变量描述延时。

S.delay = 'j*2*ms'

更加复杂的连接

对于大规模的网络我们不可能一个一个’i’和’j’显式地声明连接，通常是利用神经元索引或其他突触的变量。另外，还可以用条件（condition）、概率§、字符串表达式以及生成器建立连接关系。
条件加概率的，适用几乎全部连接的情况 ：

S.connect(condition='i!=j', p=0.2)
S.connect(condition='abs(i-j)<4 and i!=j')
S.connect(j='k for k in range(i-3, i+4) if i!=k', skip_if_invalid=True) # 生成器语法
S.connect(j='i') # 生成器语法

‘skip_if_invalid=True’可以避免边界错误，例如连接神经元1至神经元-2。值得注意的是相同的连接形式，不同的定义方法是有效率高低之分的，个人感觉你的描述方法越简单所需的计算越少则越有效率，这是一件值得思考的事情。

更加复杂的权值赋值

也可以像定义连接那样用字符串来给突触定义权值，如下例子为根据突触的空间位置来决定连接权值。

N = 30
neuron_spacing = 50*umetre
width = N/4.0*neuron_spacing

# Neuron has one variable x, its position
G = NeuronGroup(N, 'x : metre')
G.x = 'i*neuron_spacing'

# All synapses are connected (excluding self-connections)
S = Synapses(G, G, 'w : 1')
S.connect(condition='i!=j')
# Weight varies with distance
S.w = 'exp(-(x_pre-x_post)**2/(2*width**2))'

scatter(S.x_pre/um, S.x_post/um, S.w*20)
xlabel('Source neuron position (um)')
ylabel('Target neuron position (um)');

输出结果如下：

连接可视化

可以利用这个函数以两种形式可视化单层连接，直接将突触对象传给它就可以。

def visualise_connectivity(S):
    Ns = len(S.source)
    Nt = len(S.target)
    figure(figsize=(10, 4))
    subplot(121)
    plot(zeros(Ns), arange(Ns), 'ok', ms=10)
    plot(ones(Nt), arange(Nt), 'ok', ms=10)
    for i, j in zip(S.i, S.j):
        plot([0, 1], [i, j], '-k')
    xticks([0, 1], ['Source', 'Target'])
    ylabel('Neuron index')
    xlim(-0.1, 1.1)
    ylim(-1, max(Ns, Nt))
    subplot(122)
    plot(S.i, S.j, 'ok')
    xlim(-1, Ns)
    ylim(-1, Nt)
    xlabel('Source neuron index')
    ylabel('Target neuron index')
visualise_connectivity(S)

说明书 User guide

导入 import

from brian2 import *

普通的导入brain2不仅可以导入brian的module还包含了pylab，其中包括了matplotlib中关于绘图的模块和numpy/scipy所有的包。但要注意从Brain2中导入的matplotlib和scipy是单位敏感的。所以在brain2的模块导入后不能再重复导入pyplot包中的模块并构成冲突。
但也可以仅仅导入brian的模块，如下方法：

from brian2.only import *

按上面的方法导入就没有包含matplotlib和numpy/scipy。如果再需要导入pylab和numpy，brian2也提供了一些可以在后面导入的模块（经过brian修改过成为unit-aware），例如：

import brian2.numpy_ as np
import brian2.only as br2

单位 units

基本单位

见上文

温度

使用基本单位制中的kelvin，在brian2中已经去掉了celsius的单位。

常数

The brian2.units.constants 包提供了很多物理常数，供模型定义中使用。

数值in-place操作的单位

带单位的inplace操作会影响underlying array，非in-palce只影响操作变量本身。

a = [1,2]*mV
b = a 
a += 1mV
print(a)
>>> [2., 3.] mV
print(b)
>>> [2., 3.] mV

a = [1,2]*mV
b = a 
a += 1
print(a)
>>> [2., 3.] mV
print(b)
>>> [1., 2.] mV

神经元组 NeuronGroup

模型 model

神经元的定义是整个模拟的核心，神经元模型就是一组（一个或多个）方程（要求为差分方程形式）。模型方程描述的语法如下：

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

带有子表达式的描述如下：

eqs = '''
dv/dt = I_leak/Cm : volt
I_leak = g_l*(E_l -v) : amp
'''

表达式是所有NeuronGroup中定义的神经元的模型，如果有因神经元而异的参数就要单独把这个参数定义一下。

eqs = '''
dv/dt = (v_i-v)/tau : 1
v_i : 1
'''

不应期 Refractory

在brian1中，refractory有两个一并的效果：薄膜电压被钳住；神经元不进行spike。在brian2中这两个效果是可以分开施加的，包膜电压的管制是通过‘（unless refractory）’实现的，下面是brian2中对refractory的运用。

group = NeuronGroup(N, 'dv/dt = (I - v)/tau : volt (unless refractory)',
                    threshold='v > -50*mV',
                    reset='v = -70*mV',
                    refractory=3*ms)

阈值和复位操作 threshold & reset

阈值的条件一旦被满足就会执行’reset’。threshold和reset和模型的公式一样可以有外部变量、参数和单位。

v_r = -70*mV  # reset potential
G = NeuronGroup(10, '''dv/dt = -v/tau : volt
                       v_th : volt  # neuron-specific threshold''',
                threshold='v > v_th', reset='v = v_r')

阈值的参考也可以是神经元的其他状态变量甚至是时间。

访问状态变量

方法：get_states(); set_states()

 >>> group = NeuronGroup(5, '''dv/dt = -v/tau : 1
...                           tau : second''')
>>> initial_values = {'v': [0, 1, 2, 3, 4],
...                   'tau': [10, 20, 10, 20, 10]*ms}
>>> group.set_states(initial_values)
>>> group.v[:]
array([ 0.,  1.,  2.,  3.,  4.])
>>> group.tau[:]
array([ 10.,  20.,  10.,  20.,  10.]) * msecond
>>> states = group.get_states()
>>> states['v']

还可以通过Pandas数据格式（需要安装有Pandas）对状态变量(不带单位)进行访问。

>>> df = group.get_states(units=False, format='pandas')
>>> df
   N      dt  i    t   tau    v
0  5  0.0001  0  0.0  0.01  0.0
1  5  0.0001  1  0.0  0.02  1.0
2  5  0.0001  2  0.0  0.01  2.0
3  5  0.0001  3  0.0  0.02  3.0
4  5  0.0001  4  0.0  0.01  4.0
>>> df['tau']
0    0.01
1    0.02
2    0.01
3    0.02
4    0.01
Name: tau, dtype: float64
>>> df['tau'] *= 2
>>> group.set_states(df[['tau']], units=False, format='pandas')
>>> group.tau
<neurongroup.tau: array([ 20.,  40.,  20.,  40.,  20.]) * msecond>

共享变量或子式

这里的共享是指一个参数对神经元组中所有神经元相同，共享变量与外部参数不同的是它可以在模拟过程中变化，通过使用’run_rugularly()’。比如给神经元添加泊松阵列输入要保证每个神经元接收相同的输入就必须定义为共享变量。注意，对于共享变量的访问只能访问到一个变量。

G = NeuronGroup(10, '''shared_input : volt (shared)
                       dv/dt = (-v + shared_input)/tau : volt
                       tau : second''')
>>> G.shared_input = '(4.0/N)*mV'
>>> G.shared_input
<neurongroup.shared_input: 0.4 * mvolt>

注：共享变量的使用是有限制的，第一共享变量不能添加在只可能作用在一部分神经元的模型中，例如复位‘reset’操作；第二，如果还有向量变量，那么共享变量要写在前面。

链接变量

链接变量是允许一个神经元组引用存储在另外一个神经元组中的状态变量。在Brian2中链接变量必须现实地标注出来。

group1 = NeuronGroup(N,
                     'dv/dt = -v / tau : volt')
group2 = NeuronGroup(N,
                     '''dv/dt = (-v + w) / tau : volt
                        w : volt (linked)''')
group2.w = linked_var(group1, 'v')

另外一个例子：

inp = NeuronGroup(1, 'dnoise/dt = -noise/tau + tau**-0.5*xi : 1')

neurons = NeuronGroup(100, '''noise : 1 (linked)
                              dv/dt = (-v + noise_strength*noise)/tau : volt''')
neurons.noise = linked_var(inp, 'noise')

如果两个组的神经元数量相同那这种链接是一对一的，如果，源组只有一个神经元这个变量将连接所有的目标阻的神经元。如果是其他情况就必须显示地指明对应关系。

# two inputs with different phases
inp = NeuronGroup(2, '''phase : 1
                        dx/dt = 1*mV/ms*sin(2*pi*100*Hz*t-phase) : volt''')
inp.phase = [0, pi/2]

neurons = NeuronGroup(100, '''inp : volt (linked)
                              dv/dt = (-v + inp) / tau : volt''')
# Half of the cells get the first input, other half gets the second
neurons.inp = linked_var(inp, 'x', index=repeat([0, 1], 50))

子组 subgroup

Brian2中除去了’subgroup’这一方法，保留了‘切片’的方法。

group1 = G.[ :3200]

注意这里的group1中神经元的数量是3200，最大索引是3199。

突触 Synapses

Brian2中的突触与生物中实际的突触略有不同，模拟中的突触可以接收前后神经元的spike，并启动更新，生物中大部分情况突触是单向的接收信息的。更新的顺序是’on_pre’或‘on_post’代码，然后更新’model’ 中的变量。

定义突触模型

一个最简单的突触定义如下：

w = 1*mV
S = Synapses(P, G, on_pre = 'v += w')

P和G是NeuroGroup定义的神经元组，分别为源和目标，如果仅有一个神经元组，那么这个神经元组自连。关键词‘on_pre’表示的是每当有前神经元的脉冲到达突触，突触进行的动作。‘v’不是突触中定义的变量那就默认为后神经元中的变量（postsynaptic variable）。
突触模型也可以用公式来定义，最简单的形式如下。

S = Synapses(P, G, model='w:1', on_pre = 'v+=w')

上式定义了一个参数’w’，‘w’是一个synapse-specific的权值。突触也可以定义在后神经元发射脉冲后的动作，用关键词’on_post’来标明。另外还可以为每个突触接收’pre-synaptic’的脉冲定义一个固定的延迟。
**注：**突触模型中出现的变量优先指向突触中的变量，其次是神经元中的（默认是后神经元的），再其次是外部变量。如果要显示指明是来自哪个神经元，那么应该在变量后加‘_post’，如‘v_post’。

事件触发更新

通常情况下，差分方程的计算是时钟驱动的，但是这样计算是很费时间的，如果确认是需要连续记录状态，那么可以用’(clock-driven)‘来标明。另外一种是事件（spike）驱动，用’(event-driven)'标明，最典型的例子就是脉冲序列。

model = '''
w : 1
dApre/dt = -Apre/tau_pre : 1 (event-driven)
dApost/dt = -Apost/tau_post : 1 (event-driven)

事件触发的情况下，突触在接收到无论是前向还是后向神经元发出的脉冲后都会对公式进行计算，也就是说是在’on_pre’和’on_post’的代码调用之前进行计算的。另外，事件触发的计算形式只适合一小部分公式，尤其适合一维线性方程，事件触发的公式也不能依赖其它公式的结果。

pre&post code

pre&post代码分别在前神经元脉冲和后神经元脉冲的触发下执行，代码的内容可以很灵活。

S=Synapses(input,neurons,model="""w : 1
                              p : 1""",
                         on_pre="v+=w*(rand()<p)")

建立突触连接

定义了模型后需要create突触才能实现真正的连接。连接方法有如下几种方式：

S = Synapses(P, G, model='w: 1', on_pre='v+=w')
S.connect() #
S.connect(i=[2, 3], j=[5,6,7])
S.connect(i=numpy.arange(10), j=1)
S.connect(condition='abs(i-j)<=5') # reference to i, j
S.connect(condition='sqrt((x_pre-x_post)**2 + 
                   (y_pre-y_post)**2) < 250*umeter') # reference to neuron variables
S.connect(p=0.1) # probability
S.connect(condition='i != j',
          p='p_max*exp(-(x_pre-x_post)**2+(y_pre-y_post)**2) / (2*(125*umeter)**2)')
S.connect(j='i') # one to one
S.connect(j='int(i/2) if i % 2 == 0')
# generator: j='EXPR for VAR in RANGE if COND'
S.connect(j='k for k in range(0, i+1)') # range(start, end , step)
S.connect(j='k for k in sample(10, p=0.1)') # sample(N , p=x) or sample(low, high, step, p=x)
S.connect(j='i+(-1)**k for k in range(2)', skip_if_invalid=True) #skip invalid index

访问突触变量

访问突触变量参考的索引可以有三种：神经元序号的二维索引，神经元组，神经元变量的条件公式。

S.w[2, 5] = 1*nS
S.w[1, :] = 2*nS
S.w = 1*nS # all synapses assigned
S.w[2, 3] = (1*nS, 2*nS)
S.w[group1, group2] = "(1+cos(i-j))*2*nS"
S.w[:, :] = 'rand()*nS'
S.w['abs(x_pre-x_post) < 250*umetre'] = 1*nS

延时

‘delay’是Brian中预定义的突触变量，在只有’on_pre’代码的情况下是指从前神经元到突触的延迟，可以用’pre.delay’来访问（也可以用Synapse.delay访问），那当有’on_post’的情况下也可以定义’pre.post’。需要注意的是这个delay是不能’synapse-specific’的，而是统一的。要是想单独为每个突触定义，那么要将delay定义为普通变量（有个疑问就是如果定义为普通变量还能起到延迟的作用吗）。

synapses = Synapses(sources, targets, '...', on_pre='...', delay=1*ms)

突触的监测

对突触监测的基本语法如下：

M = StateMonitor(S, 'w', record=[0,1])

需要注意的是这里的’record=[0,1]‘是突触的自身序号的索引，并非神经元的索引，如果是通过神经元的索引需要通过’S’:

M = StateMonitor(S, 'w', record= S[0,1]) # neuron 0 & neuron 1
M = StateMonitor(S, 'w', record= S[0,:]) # from neuron 0
M = StateMonitor(S, 'w', record=S['i!=j'])  # all synapses excluding autapses
M = StateMonitor(S, 'w', record=S['w>0'])  # all synapses with non-zero weights (at this time)

另外对监测变量的访问方法如下：

plot(M.t / ms, M[0].w / nS)  # first synapse
plot(M.t / ms, M[S[0, :]].w / nS)  # all synapses originating from neuron 0
plot(M.t / ms, M[S['w>0*nS']].w / nS)  # all synapses with non-zero weights (at this time)

注意不能在standalone模式中使用Synapse对象索引，也不能用‘record= True’，因为此时Synapse还没生成，Brian没法索引突触。

突触的总和型变量

后神经元可以定义变量的总和，例如来自所有突触的变量，称作“summed variable”。

neurons = NeuronGroup(1, model='''dv/dt=(gtot-v)/(10*ms) : 1
                                  gtot : 1''')
S = Synapses(input, neurons,
             model='''dg/dt=-a*g+b*x*(1-g) : 1
                      gtot_post = g : 1  (summed)
                      dx/dt=-c*x : 1
                      w : 1 # synaptic weight''', on_pre='x+=w')

Another example is gap junctions:

neurons = NeuronGroup(N, model=''' dv/dt=(v0-v+Igap)/tau : 1
                                   Igap : 1''')
synapses = Synapses(neurons, model=''' w :1
                                       Igap_post=w*(v_pre-v_post): 1 (summed)''')

如果一个后神经元的突触来自多个对象，则必须分开说明，并在后神经元中相加。

neurons = NeuronGroup(1, model='''dv/dt=(gtot-v)/(10*ms) : 1
                                  gtot = gtot1 + gtot2: 1
                                  gtot1 : 1
                                  gtot2 : 1''')
S1 = Synapses(input, neurons,
              model='''dg/dt=-a1*g+b1*x*(1-g) : 1
                       gtot1_post = g : 1  (summed)
                       dx/dt=-c1*x : 1
                       w : 1 # synaptic weight
                    ''', on_pre='x+=w')
S2 = Synapses(input, neurons,
              model='''dg/dt=-a2*g+b2*x*(1-g) : 1
                       gtot2_post = g : 1  (summed)
                       dx/dt=-c2*x : 1
                       w : 1 # synaptic weight
                    ''', on_pre='x+=w')

多个突触

这里的多个突触是在一对神经元中添加多个突触，基本语法如下：

S.connect(i=numpy.arange(10), j=1, n=3)

可以定义多个突触的作用显然是探索两个神经元间不同参数连接的作用，例如有不同延时的突触，那么此时最好为这多个突触定义索引以便后续访问。

neurons = NeuronGroup(N, model='dv/dt=(v0-v)/dt', threshold='v>vt', method='exact')
synapse = Synapses(inputs, neurons, 'w: 1', on_pre='v+=w', multisynaptic_index='synapse_num')
synapse.connect(i=numpy.arange(10), j=1, n=10)
synapse.delay = '1*ms + synapse_num*2*ms'
S.w['synapse_number<5'] = 0.5
S.w['synapse_number>=5'] = 1
S.w[:, :, 5:] = 1 : 3D index

多条更新路径

来自同样一组神经元足可能存在多条更新代码路径，这对同一个前突触脉冲需要执行不同动作是相当有意义的。为了达到这个目的，用字典来定义路径的名字和需要执行的代码。

on_pre={'pre_transmission':'v+=w',
        'pre_plasticity': '''w=clip(w+Apost,0,wmax)
                             Apre+=dApre'''}

pre pathways通常情况下是在post pathways之前执行的，至于pre-pathway和post-pathway内部则是按照字母表排序执行的(i.e. pre_plasticity will be executed before pre_transmission)。如果需要显示定义，则可以如下定义：

S.pre_transmission.order = -2

这就将保证在每个时间怖中pre_transmission代码可以在pre_plasticity代码前面执行。

数字积分 numerical integration

如果要对变量进行积分，可以采用与NeuronGroup中相同的关键词‘dt’。如果运用了积分函数那么变量将在每个时间步后计算，那么也将带来复杂性的提高。

显示的事件触发更新

可以为突触变量定义显示的事件触发更新。为此，Brian内部提供了两个内部变量， ‘t’是代码执行时的时间‘lastupdate’是突触上次更新的时间 (either through on_pre or on_post code)。一个例子是short-term plasticity：

S=Synapses(input,neuron,
           model='''x : 1
                    u : 1
                    w : 1''',
           on_pre='''u=U+(u-U)*exp(-(t-lastupdate)/tauf)
                  x=1+(x-1)*exp(-(t-lastupdate)/taud)
                  i+=w*u*x
                  x*=(1-u)
                  u+=U*(1-u)''')

效率的考虑

如果只有有限的神经元需要连接，那么为连接指名前后神经元的序号是建议的，再其次用生成器；如果是一对一的形式就用‘S.connect(j=‘i’)’，这会比条件描述‘S.connect(condition=‘i=j’)’高效；当只有很少的神经元不连接时，则可以采用条件方法，条件方法的复杂度是O(N*N)。

关于模拟设置 simulation

开启模拟的方法有两种，一种是定义’Network’，然后使用对应的’Network.run()’；或者直接无脑使用’run()’，那么此时Brian将收集当前名字空间中的所有对象进行模拟。如果直接用’run()’，那么Brian将首先自动收集(magic system)显示定义的对象，让后再进行模拟。如果对象不是被显示定义的那需要进行手动添加。

G = NeuronGroup(10, 'dv/dt = -v / tau : volt')
S = Synapses(G, G, model='w:1', on_pre='v+=w')
S.connect('i!=j')
mon = SpikeMonitor(G)
run(10*ms)  # will include G, S, mon
#####################################
G = NeuronGroup(10, 'dv/dt = -v / tau : volt')
S = Synapses(G, G, model='w:1', on_pre='v+=w')
S.connect('i!=j')
monitors = [SpikeMonitor(G), StateMonitor(G, 'v', record=True)]
# a simple run would not include the monitors
net = Network(collect())  # automatically include G and S
net.add(monitors)  # manually add the monitors
net.run(10*ms)

设置模拟时间精度

时间精度必须能被模拟时间整除
全局的时间精度

defaultclock.dt = 0.05*ms

个别的对象时间精度的定义

s_mon = StateMonitor(group, 'v', record=True, dt=1*ms)

修改时间精度

全局修改可以在对象定义完成后，或者一段模拟结束后。

defaultclock.dt = 0.1*ms
# Set the network
# ...
run(initial_time)
defaultclock.dt = 0.01*ms
run(full_time - initial_time)

局部修改不能直接更改对象的’dt’属性，而是修改’clock’属性的’dt’值。To change the time step between runs for objects that do not use the defaultclock, you cannot directly change their dt attribute (which is read-only) but instead you have to change the dt of the clock attribute. If you want to change the dt value of several objects at the same time (but not for all of them, i.e. when you cannot use defaultclock.dt) then you might consider creating a Clock object explicitly and then passing this clock to each object with the clock keyword argument (instead of dt). This way, you can later change the dt for several objects at once by assigning a new value to Clock.dt.

继续或重复模拟

store()和restore()默认保存和返回的状态名字是’default’，也可以在执行时自己定义名称防止多个状态之间的混乱和覆盖，例如’store(‘initialized’); restore(‘initialized’);store(‘training’);restore(‘training’)’。

# set up the network
G = NeuronGroup(...)
...
spike_monitor = SpikeMonitor(G)

# Snapshot the state
store()

# Run the trials
spike_counts = []
for trial in range(3):
    restore()  # Restore the initial state
    run(...)
    # store the results
    spike_counts.append(spike_monitor.count)

store()中可以加入’filename’关键词变量，将状态保存在文件中（这样可以保存模拟的初始部分，在今后可以不必重复模拟），但是仅仅保存的状态，网络还需要重新定义
在Brian2中只能通过’store(), restore(), run()'影响模拟时间，‘Network.t’是只读变量，想要回到之前的模拟时间，比如要加入新的噪声实例，那么必须使用’store(), restore(), run()’。

模拟规划

在Brian2中每个模拟对象都有三个属性’dt, when, order’。在一个时间步中，模拟对象首先根据它们在’when’变量中的位置进行更新。这个规划的顺序决定于’Network.schedule’（是一个字符串的列表决定执行计算的位置和顺序’execution slots’），默认为 ‘start’, ‘groups’, ‘thresholds’, ‘synapses’, ‘resets’, ‘end’，除了这些，名字还可以为类似于’before_thresholds’ 或者’after_synapses’对应于具体的位置。'when’默认属性的对大部分对象来说是敏感值(resets will happen in the reset slot, etc.)，但有的时间改变它是有意义的， 例如如果一个人想StateMonitor（which by default records in the end slot）在复位前记录胞膜电压，如果不这样设置将不会看到胞膜电压超越阈值。在一个’execution slots’中仿真顺序由’order’决定，'order’相同由词典顺序决定。
要查看一个网络的对象是如何规划的，可以通过scheduling_summary() 函数查看。

>>> group = NeuronGroup(10, 'dv/dt = -v/(10*ms) : 1', threshold='v > 1',
...                     reset='v = 0')
>>> mon = StateMonitor(group, 'v', record=True, dt=1*ms)
>>> scheduling_summary()
                object                  |           part of           |        Clock dt        |    when    | order | active
----------------------------------------+-----------------------------+------------------------+------------+-------+-------
statemonitor (StateMonitor)             | statemonitor (StateMonitor) | 1. ms (every 10 steps) | start      |     0 |  yes
neurongroup_stateupdater (StateUpdater) | neurongroup (NeuronGroup)   | 100. us (every step)   | groups     |     0 |  yes
neurongroup_thresholder (Thresholder)   | neurongroup (NeuronGroup)   | 100. us (every step)   | thresholds |     0 |  yes
neurongroup_resetter (Resetter)         | neurongroup (NeuronGroup)   | 100. us (every step)   | resets     |     0 |  yes

一次模拟中多个run()

如果模拟的对象均是新定义的，那么从时间0开始一段新的模拟；如果全是老对象，那么就从结束的时间继续模拟；如果需要混合新旧对象，那么必须显示地定义’Network’。在C++ code generation模式下需要有另外的要求和设置。

模拟过程的输出

主要是对’report’和’report_period’的定义，‘report’可以是’stdout’、‘text’和’stderr’。'report_period’可以定义除了开始和结束，输出模拟进程的间隔。需要时也可以进行定制的进程输出。

模拟耗时详情输出

如果要得到模拟的详情，就必须在’run()'时加入’profile=True’的关键词变量。模拟完成后可以通过’Network.profiling_info’取回，通过’profiling_summary(show=5)'打印出来。

>>> profiling_summary(show=5)  # show the 5 objects that took the longest
Profiling summary
=================
neurongroup_stateupdater    5.54 s    61.32 %
synapses_pre                1.39 s    15.39 %
synapses_1_pre              1.03 s    11.37 %
spikemonitor                0.59 s     6.55 %
neurongroup_thresholder     0.33 s     3.66 %

数值积分的方法

默认是自动选择计算方法，但因为默认的不能提前知晓参数的情况，所以最好是自己说明。如果选择的计算方法不合适将有可能导致计算错误（not a number or large oscillation）。
方法选择
模拟开始后你将得到若干条 INFO信息来告诉你模拟器最终使用的哪种积分方法以及其他几种积分方法所用的时间。可以手动设置的方法：

• ‘exact’: exact integration for linear equations (alternative name: ‘linear’)
• ‘exponential_euler’: exponential Euler integration for conditionally linear equations
• ‘euler’: forward Euler integration (for additive stochastic differential equations using the Euler-Maruyama method)
• ‘rk2’: second order Runge-Kutta method (midpoint method)
• ‘rk4’: classical Runge-Kutta method (RK4)
• ‘heun’: stochastic Heun method for solving Stratonovich stochastic differential equations with non-diagonal multiplicative noise.
• ‘milstein’: derivative-free Milstein method for solving stochastic differential equations with diagonal multiplicative noise

Debug

编译出现错误要从最后一行看起，最后一行通常会指名错误的性质和原因。再向上看就是和错误有关的函数和语句（包括具体的文件以及用箭头标示的行），这些与错误相关的语句是层次追溯的（traceback），一般情况错误在于最顶层的函数或语句，也就是自己的脚本。如果排除了自己的错误，那么可以联系官方修正BUG。

总结

需要关注的重点，对模拟的准确性和高效性有影响的点

• 编译的规划：standalone or runtime generation;
• 执行实际的计算： Python(numpy code);Cython or C; standalone code using OpenMP(in development with limited function) 。调代码和探索阶段建议使用Python code generation，需要进一步大规模或多次模拟时使用Cython code generation；有条件可以尝试OpenMP。
• 模拟规划：Network.schedule，[‘start’, ‘groups’, ‘thresholds’, ‘synapses’, ‘resets’, ‘end’]
• 可以查看哪个模块模拟占用的时间最多：profile=True，profiling_summary()
• 数值积分方法的选择：默认是自动选择计算方法，但因为默认的不能提前知晓参数的情况，所以最好是自己说明。
• 模拟对象的定义：Network
• 模拟的拆分：store(); restore()
• 共享变量：（shared_variable）

容易导致隐性错误的点
inplace operations：在python中inplace operation会连同underlying的数列一起改变，在brian中也是如此，但是当inplace operation没有带单位时，只改变参与运算的数列本身。

写在最后

我写博客的原则：

• 问题的系统性阐述：简洁、美观、全面
• 关键问题的阐述：准确、清晰
  写完后对照看一下是否是粗制滥造的一篇

---

1. Diehl P U, Cook M. Unsupervised learning of digit recognition using spike-timing-dependent plasticity.[J]. Frontiers in Computational Neuroscience, 2015: 99-99. ↩︎