CARLsim4:http://uci-carl.github.io/CARLsim4/index.html
CARLsim4代码:https://github.com/UCI-CARL/CARLsim4
前言
CARLsim4是一个高效、易用的、基于GPU加速的大规模脉冲神经网络仿真系统。由“[Cognitive Anteater Robotics Laboratory]
这一节主要介绍脉冲神经网络的几个核心概念——神经元、连接、突触。脉冲神经网络的一大特点就是生物学概念还原的很好,而传统的神经网络已经是高度抽象的结果。
Cognitive Anteater Robotics Laboratory
开源,C++编写,支持windows和linux平台。
即将入坑SNN,新开一个关于CARLsim的主题,主要的内容都是基于CARLsim4文档和源码,记录自己的理解、总结、翻译和搬运。
神经元
CARLsim目前支持基于电流和基于电导的Izhikevich脉冲神经元。
使用函数createSpikeGeneratorGroup
和createGroup
可以创建从一维到三维的神经元组。
微分方程可以由forward-Euler法或Runge-Kutta法实现。
1.1 Izhikevich神经元(4参数模型)
Izhikevich是一个动力学系统模型,可以由以下二阶常微分方程描述:
d v d t = 0.04 v 2 + 5 v + 140 − u + I \begin{array}{l} \frac{{dv}}{{dt}} = 0.04{v^2} + 5v + 140 - u + I\\ \end{array} dtdv=0.04v2+5v+140−u+I
d u d t = a ( b v − u ) \begin{array}{l} \frac{{du}}{{dt}} = a(bv - u) \end{array} dtdu=a(bv−u)
第一个公式描述了给定电流 I I I的膜电位 v v v,其中 I I I是所有突触和外部电流之和( I = i s y n + I e x t I = {i_{syn}} + {I_{ext}} I=isyn+Iext);
第二个式子表述了恢复变量 u u u; a a a是恢复变量的速率常数; b b b描述了恢复变量对于膜电位阈值下波动的敏感性;两个公式的变量都是无量纲的,但是第一个公式的 v v v是 m V 级 mV级 mV级的, t t t是 m s ms ms级的。
和其他简单的模型(如LIF)不同,Izhikevich神经元可以自行生成上行脉冲,因此电压在峰值(
v
c
u
t
o
f
f
=
+
30
{v_{cutoff}} = + 30
vcutoff=+30)复位,而不是在阈值。
当
v
v
v达到峰值截断时,动作电位下行通过瞬时膜电位复位建模,并加上回复变量的步进(这句话很拗口,但是看了下面的公式,就能秒懂),具体由下面的公式描述:
v ( v > 30 ) = c \begin{array}{l} v(v > 30) = c \end{array} v(v>30)=c
u ( v > 30 ) = u + d \begin{array}{l} u(v > 30) = u + d \end{array} u(v>30)=u+d
不同的神经元abcd取值不同,如
规则脉冲(RS)取值
a
=
0.02
a=0.02
a=0.02,
b
=
0.2
b=0.2
b=0.2,
c
=
−
65
c=-65
c=−65,
d
=
8
d=8
d=8;
快速脉冲(FS)取值
a
=
0.1
a=0.1
a=0.1,
b
=
0.2
b=0.2
b=0.2,
c
=
−
65
c=-65
c=−65,
d
=
2
d=2
d=2。
1.2 Izhikevich神经元(9参数模型)
对于9参数Izhikevich模型定电流 I I I的膜电位 v v v,如下式:
d v d t = ( k ( v − v r ) ( v − v t ) − u + I ) / C \frac{{dv}}{{dt}} = (k(v - vr)(v - vt) - u + I)/C dtdv=(k(v−vr)(v−vt)−u+I)/C
其中, K K K是膜电位的速度常量,由神经元的输入电阻和基强度(rheobase,足引起刺激的最小强度), v r vr vr是静止膜电位, v t vt vt是瞬时阈电位(threshold potential), C C C是膜容。
对于恢复变量
u
u
u,如下式:
d
v
d
t
=
a
(
b
(
v
−
v
r
)
−
u
)
\frac{{dv}}{{dt}} = a(b(v - vr) - u)
dtdv=a(b(v−vr)−u)
对于动作电位下行建模如下:
v
(
v
>
v
p
e
a
k
)
=
c
\begin{array}{l} v(v > vpeak) = c \end{array}
v(v>vpeak)=c
u
(
v
>
v
p
e
a
k
)
=
u
+
d
\begin{array}{l} u(v > vpeak) = u + d \end{array}
u(v>vpeak)=u+d
所谓的九变量是指:
C
C
C,
k
k
k,
v
r
vr
vr,
v
t
vt
vt,
a
a
a,
b
b
b,
v
p
e
a
k
vpeak
vpeak,
c
c
c,
d
d
d
RS神经元:
C
=
100
C=100
C=100,
k
=
0.7
k=0.7
k=0.7,
v
r
=
−
60
vr=-60
vr=−60,
v
t
=
−
40
vt=-40
vt=−40,
a
=
0.03
a=0.03
a=0.03,
b
=
−
2.0
b=-2.0
b=−2.0,
v
p
e
a
k
=
35
vpeak=35
vpeak=35,
c
=
−
50
c=-50
c=−50,
d
=
100
d=100
d=100
FS神经元:
C
=
20
C=20
C=20,
k
=
1
k=1
k=1,
v
r
=
−
55
vr=-55
vr=−55,
v
=
−
40
v=-40
v=−40,
a
=
0.15
a=0.15
a=0.15,
b
=
8
b=8
b=8,
v
p
e
a
k
=
25
vpeak=25
vpeak=25,
c
=
−
50
c=-50
c=−50,
d
=
200
d=200
d=200
一般在九参数的Izhikevich模型中,使用Runge-Kutta而不是forward-Euler。
1.3 Multi-Compartment神经元
使用CARLsim::setCompartmentParameters
函数将一组神经元拓展到多个单元,使用 CARLsim::connectCompartments
连接Compartment。
这时,总电流 I I I被进一步扩展,包含了树突电流,即:
I = i s y n + i d e n d r + i e x t I=i_{syn}+i_{dendr}+i_{ext} I=isyn+idendr+iext,
每个单元的树突电流来自(down,mother)下单元和(up,daughter)上单元:
i d e n d r = G d o w n ( v − v d o w n ) + ∑ u p G u p ( v − v u p ) {i_{dendr}} = {G_{down}}(v - {v_{down}}) + \sum\limits_{up} {{G_{up}}} (v - {v_{up}}) idendr=Gdown(v−vdown)+up∑Gup(v−vup)
电导率
G
d
o
w
n
G_{down}
Gdown和
G
u
p
G_{up}
Gup由CARLsim::setCompartmentParameters
确定。
同样的,这里推荐使用Runge-Kutta而不是forward-Euler。
这一段直接翻译了CARlsim的文档, 我根本就没有就看懂在讲什么,Multi-Compartment是什么鬼,下图附一个轴突、树突图来大胆猜想
突触
脉冲神经网络使用突触传递信息。生物学中,突触通过化学信号在神经元之间传递信息,(电突触也是存在的,但是Carlsim没有实现)。
通常来说,突触传递信息是单项的。
突触分为两中,一种在神经元发送信息(前突触神经元),另一种在神经元接受信息(突触后神经元)。
这两种神经元被突出间隔隔开。突触前神经元释放神经递质和突触后神经元受体结合,产生电流。当多个突触的产生的电流之和大于阈值电压,则产生脉冲。
下图为高中生物中介绍的突触结构:
CARLsim支持两种突触模型,基于电流的模型(CUBA)只考虑突触电流;基于电导的模型使用更复杂的电导方程计算每个突触受体电流。
CUBA和COBA都受突触权重的影响。
接下来将详细讨论这两个模型。
2.1 CUBA
使用CUBA模型时,不考虑电导。产生的电流强度和突触权重成正比。
突触后神经元产生的电流 I j s y n I_j^{syn} Ijsyn由突触前神经元产生脉冲时给出:
I j s y n = ∑ i = 1 N s i j w i j I_j^{syn} = \sum\limits_{i = 1}^N {{s_{ij}}} {w_{ij}} Ijsyn=i=1∑Nsijwij
上式中, S i j = 1 S_{ij}=1 Sij=1意味着神经元产生了脉冲, S i j = 0 S_{ij}=0 Sij=0则没有; w i j w_{ij} wij表示突触后神经元 j j j和突触前神经元 i i i之间的权值; N N N是突触前神经元和突触后神经元 j j j产生的连接数量。
基于电流的突触模型是CARLsim的默认模型,目前CARLsim不支持COBA模型和其他模型混用。
基于电流的神经元模型电流是一次性产生的,而基于电导的神经元模型的突触电流随时间衰减。
由于COBA模型产生的电流要大得多,当从CUBA模型转换到COBA时,需要用户去减少突触权重。
2.2 COBA
当使用COBA模型时,存在指数衰减的电导(电导,表示某一种导体传输电流能力强弱程度)。如果突触连接是兴奋的,则使用AMPA和NMDA衰减。如果突出连接是抑制的,则使用GABA A和GABA B连接。
总电流既可以来自兴奋区( i e {i_e} ie),也来自抑制区( i e {i_e} ie)。
i
e
=
i
N
M
D
A
+
i
A
M
P
A
{i_e} = {i_{NMDA}} + {i_{AMPA}}
ie=iNMDA+iAMPA
i
i
=
i
G
A
B
A
A
+
i
G
A
B
A
B
{i_i} = {i_{GAB{A_A}}} + {i_{GAB{A_B}}}
ii=iGABAA+iGABAB
兴奋电流由NMDA电流( i N M D A i_{NMDA} iNMDA)和AMPA电流( i A M P A {i_{AMPA}} iAMPA )组成:
i A M P A = g A M P A ( v − v A M P A r e v ) {i_{AMPA}} = {g_{AMPA}}(v - v_{AMPA}^{rev}) iAMPA=gAMPA(v−vAMPArev)
i N M D A = g N M D A [ v + 80 60 ] 2 1 + [ v + 80 60 ] 2 ( v − v N M D A r e v ) {i_{NMDA}} = {g_{NMDA}}\frac{{{{[\frac{{v + 80}}{{60}}]}^2}}}{{1 + {{[\frac{{v + 80}}{{60}}]}^2}}}(v - v_{NMDA}^{rev}) iNMDA=gNMDA1+[60v+80]2[60v+80]2(v−vNMDArev)
上式中 g g g是电导, v v v是突触后神经元电压, v r e v vrev vrev是逆转电位, g g g和 v r e v vrev vrev都是针对特定的离子通道和受体;
抑制神经元通道由GABA A和GABA B电流组成:
i G A B A A = g G A B A A ( v − v G A B A A r e v ) \begin{array}{l} {i_{GAB{A_A}}} = {g_{GAB{A_A}}}(v - v_{GAB{A_A}}^{rev})\end{array} iGABAA=gGABAA(v−vGABAArev)
i G A B A B = g G A B A B ( v − v G A B A B r e v ) \begin{array}{l} {i_{GAB{A_B}}} = {g_{GAB{A_B}}}(v - v_{GAB{A_B}}^{rev}) \end{array} iGABAB=gGABAB(v−vGABABrev)
每一个粒子通道或受体 k k k的电导 g g g,由下式给出:
g k = ∑ f e t − t f / τ Θ ( t − t f ) {g_k} = \sum\limits_f {{e^{t - {t_f}/\tau }}} \Theta (t - {t_f}) gk=f∑et−tf/τΘ(t−tf)
上式中 f f f指一个脉冲事件, t f tf tf是脉冲时间的时间, T h e t a Theta Theta指Heaviside函数。
函数sim.setConductances(true)
表示使用COBA模式,该模式默认参数:
t
d
A
M
P
A
=
5
m
s
tdAMPA=5ms
tdAMPA=5ms,
t
d
N
M
D
A
=
150
m
s
tdNMDA=150ms
tdNMDA=150ms,
t
d
G
A
B
A
a
=
6
m
s
tdGABAa=6ms
tdGABAa=6ms,
t
d
G
A
B
A
b
=
150
m
s
tdGABAb=150ms
tdGABAb=150ms(瞬时上升时间)
也可以自定义参数:
sim.setConductances(true,tdAMPA,tdNMDA,tdGABAa,tdGABAb);
还可以自定义下降时间常数(我不知道是啥),上升时间常数(默认为0):
sim.setConductances(true,tdAMPA,tdNMDA,trNMDA,tdGABAa,tdGABAb,trGABAb)
也可以修改默认电导时间常数:
sim.setDefaultConductanceTimeConstants(tdAMPA,tdNMDA,trNMDA,tdGABAa,tdGABAb,trGABAb);
组(group)
要创建一个Izhikevich神经元组,需要参数:名字(如output),神经元数量(如100):
int gOut = sim.createGroup("output", 100, EXCITATORY_NEURON);
EXCITATORY_NEURON
:组中神经元是谷氨酸能的
INHIBITORY_NEURON
:GABAergic突触神经元
返回组号gOut
给Izhikevich设置参数,举一个列子,给第一类兴奋神经元(规则神经元)设置参数:
sim.setNeuronParameters(gOut, 0.02f, 0.2f, -65.0f, 8.0f);
分别对应前面介绍的四参数Izhikevich神经元的 a a a, b b b, c c c, d d d参数。
创建一组脉冲生成器(名称、大小、类型)
int gIn = sim.createSpikeGeneratorGroup("input", 10, EXCITATORY_NEURON);
3.1 Topography
Grid3D结构体可以将一组神经元排列三维网格(以三维结构的中心为原点),并根据通过CARLsim::connect
连接神经元。这允许创建具有复杂空间结构的网络。
具体的排列规则如下:
Grid3D(Nx,Ny,Nz)
表示坐标范围
[
−
(
N
x
−
1
)
/
2
,
(
N
x
−
1
)
/
2
]
[-(Nx-1)/2, (Nx-1)/2]
[−(Nx−1)/2,(Nx−1)/2],
[
−
(
N
y
−
1
)
/
2
,
(
N
y
−
1
)
/
2
]
[-(Ny-1)/2, (Ny-1)/2]
[−(Ny−1)/2,(Ny−1)/2],和
[
−
(
N
z
−
1
)
/
2
,
(
N
z
−
1
)
/
2
]
[-(Nz-1)/2, (Nz-1)/2]
[−(Nz−1)/2,(Nz−1)/2]。
Grid3D(1,1,1)
在位置$ (0,0,0)$创建一个神经元
Grid3D(2,1,1)
创建两个神经元, 第一个 (ID 0) 在$ (-0.5,0,0)$, 第二个 (ID 1) 在
(
0.5
,
0
,
0
)
(0.5,0,0)
(0.5,0,0).
Grid3D(1,1,2)
创建两个神经元, 第一个 (ID 0) 在
(
0
,
0
,
−
0.5
)
(0,0,-0.5)
(0,0,−0.5), 第二个 (ID 1) 在
(
0
,
0
,
0.5
)
(0,0,0.5)
(0,0,0.5).
Grid3D(2,2,2)
创建八个神经元, 第一个 (ID 0) 在
(
−
0.5
,
−
0.5
,
−
0.5
)
(-0.5,-0.5,-0.5)
(−0.5,−0.5,−0.5), 第二个
(
0.5
,
−
0.5
,
−
0.5
)
(0.5,-0.5,-0.5)
(0.5,−0.5,−0.5), 第三个
(
−
0.5
,
0.5
,
−
0.5
)
(-0.5,0.5,-0.5)
(−0.5,0.5,−0.5), 以此类推,如上图
Grid3D(3,3,3)
创建
3
x
3
x
3
=
27
3x3x3=27
3x3x3=27 个神经元,$ (-1,-1,-1)
,
,
,,
,,(0,-1,-1)
,
,
,(1,-1,-1)$,
(
−
1
,
0
,
−
1
)
(-1,0,-1)
(−1,0,−1),…,最后一个是
(
1
,
1
,
1
)
(1,1,1)
(1,1,1).
神经元的3D位置可以用CARLsim::getNeuronLocation3D
函数查询。
创建一个有500个神经元的3D神经元组
int gOut = sim.createGroup("output", Grid3D(10,10,5), EXCITATORY_NEURON);
以下两行代码等价
int gOut = sim.createGroup("output", N, EXCITATORY_NEURON);
int gOut = sim.createGroup("output", Grid3D(N,1,1), EXCITATORY_NEURON);