【CARLsim和脉冲神经网络】神经元、突触、组

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