帕金森状态的脉冲神经网络建模及最优控制
上QQ阅读APP看本书,新人免费读10天
设备和账号都新为新人

2.1 单神经元脉冲模型

我们首先从神经元生物模型开始,认识神经元的特性和建模机理。

2.1.1 神经元生物模型

研究发现,人的大脑是由数以十亿计的神经元组成的,神经元是大脑处理信息的基本单元,主要由细胞核(胞体)、树突和轴突组成(见图2.1)。神经元输出的有电信号和化学信号,其中最主要的是沿着轴突细胞膜表面传播的电脉冲。树突和轴突都有大量分支,轴突末端通常连接到其他神经元细胞的树突上,连接点被称为突触,一个神经元的输出通过突触传递给成千上万个下游神经元。有的突触促进下游细胞的兴奋,被称为兴奋性突触;有的突触抑制下游细胞的兴奋,被称为抑制性突触。同样,一个神经元也会接收来自上游成千上万个神经元的突触输入,神经元会对这些突触输入进行整合,当这些输入超过神经元的放电阈值时,神经元细胞就会产生一个动作电位(Action Potential),这个动作电位也被称为脉冲(Spike)。

img

图2.1 神经元细胞的结构示意图

事实上,神经元之间是不能直接传递脉冲序列的,神经元的输入脉冲改变神经元的突触后电压,当这个电压超过阈值时,才产生一个动作电位。因此,神经元之间的脉冲传递是一个“脉冲—电压—脉冲”的转化过程。那么神经元是如何用脉冲序列表达信息的呢?即怎样将激励信号转化为脉冲序列,在多数的建模研究中,都将激励信号等价为神经元接收到的输入电流[57]。本章中,我们将重点介绍几种神经元的数学模型。

2.1.2 单间室电生理模型

神经元电生理模型是指从神经元膜的生物物理机制出发,建模单个离子通道对神经元兴奋性及动力学特性的影响效果,不同离子通道在激活电压、失活电压、电导等电生理特性上有显著差异,这些差异使不同的离子在神经信息的传递中发挥不同的作用。通过对神经元进行电生理建模,能研究不同种类离子通道在神经元兴奋性的产生及传导中的作用效果,按照电生理特性的划分,神经元电生理模型又分为单间室电生理模型和双间室电生理模型,在单间室电生理模型中,最经典的刻画离子特性的电生理模型就是Hodgkin和Huxley建立的Hodgkin-Huxley(HH)模型。

1.Hodgkin-Huxley模型

1952年,英国生物学家Hodgkin和Huxley利用电压钳技术将膜电流区分为Na+电流、K+电流,结合非线性动力学特性建立了能够描述长枪乌贼巨大轴突的微分方程,称为Hodgkin-Huxley(HH)方程,并因此获得了诺贝尔生理学或医学奖。HH模型与电缆理论方程相结合时,可以准确模拟乌贼神经的许多电特性,如动作电位的形态、阻抗的变化、传导速度及离子交换等现象。因此HH模型被认为是最接近实际生物神经元的数学模型之一,它的各项参数都具有明确的生理意义,能准确描述神经元脉冲序列的产生和传递机理[58]。HH模型以电容和电阻的模式构建表示神经元信息传递的电路模型,如图2.2所示。图中C表示神经元的膜电容,每个并行通路表示一种离子电流的动态,分别由电阻和该离子的平衡电势组成,I(t)表示来自其他神经元的突触输入电流。

img

图2.2 HH模型的电路图

HH模型的动态由一个四阶多变量系统组成,微分方程如下:

img

(2.1)

其中

img

(2.2)

式中,V表示神经元的膜电压;C表示膜电容;I(t)表示外部施加电流;gNagKgL分别表示Na+电导、K+电导及漏电导,反映了离子通道在细胞膜上的分布密度;VNaVKVL分别对应离子的平衡电势;变量mhn分别表示电压相关的Na+通道激活门变量、Na+通道失活变量、K+通道激活变量;ambmahbhanbn表示相应离子与V相关的门变量速率函数。

图2.3描述了HH模型中动作电位的产生机制。发现对于HH模型而言,当神经元处于静息态时,该模型的膜电压为0mV,这也就是说流向细胞内的电流与流向细胞外的电流处于平衡状态,因此净电流为0。静息态是稳定的,因为当一个很小的脉冲I(t)注入神经元时,会对细胞膜电压产生一个较小的正扰动,称为去极化,相应地导致一个很小的净电流,此电流驱使细胞膜电压恢复到静息态电压。然而当较大的电流脉冲注入神经元时会产生一个显著放大的扰动,因为膜电导系数依赖于膜电压V。这样非线性的放大导致膜电压从静息态迅速脱离,此现象称为动作电位或者脉冲。同时,图2.3也展示了HH模型系统中动作电位的一个典型时间过程。较强的膜电压去极化使激活变量mn增大而使失活变量h减小,因此激活变量m的时间常数τm(V)非常小,激活变量m相对比较大,从而使Na+通道迅速激活,激活的Na+通道使Na+大量流入细胞,因此细胞膜电压趋近于Na+的反电势(ENa),这就导致了加深的去极化并进一步激活了Na+电导。这个过程形成一个正反馈,从而使膜电压上升也就是动作电位的上升支。当膜电压达到Na+反电势时,门变量打开速度变慢。变量h趋近于0,导致Na+电流的失活,以及变量n趋向于1,也就是说流向细胞外的K+电流缓慢激活。随后漏电流使膜电压复极化趋向静息电位。

当膜电压趋向于静息电位时,电压敏感时间常数τn(V)、τh(V)比较大,因此变量nh恢复缓慢。特别的是,外流K+电流持续被激活,直到动作电位下降,从而驱使膜电压低于静息态,接近K+电流反电势,这个现象称为后超极化。除此之外,在变量h比较小时Na+电流持续失活,因此神经元不具备持续放电的功能。对于HH模型而言它不能在绝对恢复期产生另一个动作电位。如果刺激电流比较大,电流去失活,系统恢复产生动作电位的能力。

因为HH模型能复现很多复杂的电生理现象,因此被认为是建模电生理神经元的经典模型,并且具有深远的影响,主要分为3个方面。一是在分子水平,HH模型构建了描述离子通道结构、功能特性的框架,如离子渗透、选择等机制;二是在细胞水平,HH模型不仅可以预测动作电位,而且可以预测控制动作电位起始时刻的条件,如阈值、恢复期等;三是在电路水平,HH模型的成功预测使它在实验方面成为利用基于实验数据建模的范例,为现在蓬勃发展的计算神经科学奠定了基础。因此HH模型在模拟某些生物结构方面意义非凡且应用广泛,如Pospischil及合作者利用最小限度的HH模型建立了皮层和丘脑神经元[59],或者用来分析神经元的动力学来阐释神经元的编码机制[60]。但是,它具有维数高、变量多的高复杂性特点,因此,不适合大规模网络仿真的建模。尽管如此,HH模型依然产生了比较深远的影响。

img

图2.3 HH模型中的动作电位[66]

2.Morris-Lecar模型

Morris-Lecar模型是C.Morris等人在1981年利用电流钳位对北极鹅的肌肉纤维进行实验时发现的[61],他们认为在北极鹅的纤维中存在一种简单的电导特性,而这种特性与肌肉纤维细胞膜中的Ca2+和K+相关。基于这种特性,C.Morris等人建立了Morris-Lecar模型,初始的电导模型用一个三阶的微分方程表示,虽然这个模型能很好地解释复杂的电压振荡行为,但对它进行数值仿真需要许多参数的调节,而且要从三阶系统的数值研究中获得电压振荡的必备条件是非常困难的。因此C.Morris等人将三阶模型简化成二阶Morris-Lecar模型,本节先介绍其三阶模型,然后再介绍其二阶模型。

Morris-Lecar模型具有与HH模型相似的电路结构,其等效电路图如图2.4所示。根据电路图,可得到描述肌肉纤维电导特性的数学表达式:

img

(2.3)

img

图2.4 Morris-Lecar模型的等效电路图

式中,MN分别表示Ca2+和K+浓度的门变量;τM(V)、τN(V)分别表示Ca2+和K+浓度变化的时间变量;I为外部输入电流,作为模型的分岔参数,当输入电流变化时,模型表现出不同的放电特性。其中,M(V)、N(V)及时间变量τM(V)、τN(V)又可以表示为:

img

(2.4)

如前所述的三阶Morris-Lecar模型能较好地解释复杂的电压振荡行为,但从中获得电压振荡的必要条件非常困难。C.Morris等人[61]通过研究发现,Ca2+和K+在变化过程中具有不同的时间尺度,即时间变量呈10倍以上差异。因此,利用奇异摄动理论,可以对具有不同时间尺度的变量进行化简。C.Morris等人认为Ca2+系统比K+系统变化快很多,以至于在任意时刻的稳定状态下满足Ca2+电导是一个稳态值,即满足M=M(V)。在此基础上,式(2.3)可以简化为一个二阶模型,即后来被大量用于模型研究的Morris-Lecar等式:

img

(2.5)

其中的稳定状态变量定义为:

img

(2.6)

式(2.5)和式(2.6)中,gL=0.5、gK=2、gCa=1、V1=-1mV、V2=15mV、V3=10mV、VCa=100mV、VL=-50mV、V4=14.5mV、VK=−70mV。这里我们主要考虑模型呈单峰放电状态时的动态行为,因此取外部输入电流I=9mA。

2.1.3 双间室电生理模型

为了探索电刺激对神经元动态的影响,需要在神经元建模时考虑电场的影响,而电场刺激被认为是具有空间分布的特征,以上所述HH电导模型属于单间室电生理模型(简称“单间室模型”),不能体现出神经元的空间分布特征,因此不能有效建模电场效应。双间室电生理模型(简称“双间室模型”)将空间上连续的单个神经元的胞体和树突建成双间室的形式,能描述神经元的时空变化,尤其可以量化外部电场对神经元的场效应。

1.Pinsky-Rinzel模型

大量实验表明,当抑制减少时,海马体CA3区域会产生同步癫痫放电。目前已有数学模型可以建模CA3神经元的本质特性及突触交互等行为,如Traub等人建立了19间室的CA3区域锥体神经元的计算模型[62]。在此基础上,Pinsky与Rinzel[63]建立了锥体细胞Pinsky-Rinzel(PR)双间室模型,此简化模型将Na+放电分为近端胞体间室的快速电流、树突间室的慢速Ca2+电流及Ca2+调节电流。在每个模型的周期簇放电被替代为重复胞体放电,当胞体刺激电流增加时,稳定的树突刺激可以产生更加显著的周期簇放电,其频率范围为8~20Hz;而胞体的簇放电频率小于8Hz。PR模型中的簇放电现象出现在电子耦合电导处于中间位置时,且与离子通道类型及胞体与树突间室之间的耦合电流有关。当胞体和树突的耦合强度较大时,此模型简化为单间室模型,则不会产生簇放电现象。用兴奋α-氨基-3-羟基-5-甲基-4-异恶唑丙酸受体(AMPA受体)和N-甲基-D-天冬氨酸受体(NMDA受体)突触连接的PR模型结果和Traub等人建立的19间室复杂电缆模型相似。处于静息态的网络单个神经元的简单刺激会产生多种同步集群簇放电,快速AMPA受体突触对同步机制贡献巨大,且簇放电的数目随着最大NMDA受体电导水平的增加而增加。

此模型可以很好地阐释CA3锥体神经元的簇放电、适应性等动力学特性。更重要的是,此模型能够很好地描述锥体神经元树突上分布的慢尺度Ca2+电流对簇放电产生的重要影响。双间室PR模型包含树突间室和胞体间室,其机制如图2.5所示。

img

图2.5 Pinsky-Rinzel(PR)模型机制[63]

图2.5中模型的动力学方程式如下:

img

(2.7)

式中,VSVD分别表示胞体间室、树突间室的膜电压;C表示膜电容;ISID分别表示应用到胞体间室、树突间室的电流;p表示胞体面积与整个细胞面积的比例。胞体间室电流包括压控内流Na+电流INa和外流延迟修正K+电流IK-DR。Na+电流瞬时激活,mm(VS),Na+通道的失活变量h和延迟修正K+通道的激活变量n在几毫秒的时间常数内迅速上升。胞体间室拥有3种压控电流,分别是内流Ca2+电流ICa、Ca2+激活K+电流IK-CIK-AHP、突触电流IAMPAINMDAICa的激活变量s是快速激活变量,IK-C正比于c倍快速激活变量s的饱和函数X(Ca);IK-AHP拥有慢速激活变量q,此变量依赖Ca2+电流而不是压控变量。q的时间常数在低Ca2+的1000ms至高Ca2+的100ms内变化。快速激活变量的有效时间常数约为6ms。AMPA受体的电容上升很快,延迟时间较短;NMDA受体绑定(Receptor Binding)上升迅速,但是延迟时间较长,时长约为150ms。所有电流性质总结如下:

img

(2.8)

式中,变量hnscq满足下式:

img

(2.9)

当变量yhn时,UVS;当变量y表示sc时,UVD。对于Ca2+y代表q。Ca2+动力学方程及NMDA受体突触电导权重Si(t)和AMPA受体突触电导权重Wi(t)的动力学方程为:

img

(2.10)

式中,img表示所有突触前神经元的代数和;H(x)表示Heaviside阶跃函数,当变量x不小于0时,H(x)恒为1,当x小于0时,H(x)恒为0。

Park等人在2005年通过实验数据对模型中的参数进行匹配,其中门变量率函数及其他参数取值分别用表2.1和表2.2表示。

表2.1 离子通道门变量率函数[100]

img

表2.2 PR模型参数

img

续表

img

2.双间室模型的场效应模型

由于电场刺激具有空间特性,为了建模神经元的电场效应,科学家构建了如图2.6所示的神经元电阻阵列模型,即双间室模型的场效应模型,该模型由两个并联的双间室神经元电导模型通过K+电流耦合嵌入一个电阻阵列中。Vapp表示外部电场的电压,Z表示各种等价电阻,图中TD、DS、SG、SS、DD分别表示顶部-树突连接、树突-胞体连接、胞体-接地连接、胞体-胞体连接、树突-树突连接,具体参数参见文献[64]。

img

图2.6 神经元电阻阵列模型

3.简化双间室模型

双间室模型的场效应模型是一种简单的神经元模型,它可以刻画神经元的动力学特性,却不能表现神经元形态对神经元行为的影响。目前,伊国胜等人在PR模型的基础上,提出了简化双间室模型,它能很好地将两者结合起来,既可以表示神经元形态,也可以刻画其动力学特性[65],如图2.7所示。

img

图2.7 简化双间室模型结构

简化神经元中的两个间室分别表示神经元的树突和胞体,相应间室的膜电压用VDVS表示。胞体和树突之间用内连电导gc连接。其中,胞体间室含有3种离子电流,分别是流向胞外的K+主动电流IK、流向胞内的主动Na+电流INa和流向胞外的漏电流IL。为了体现树突形态特性和树突间室膜电压VD,树突间室只有一种流向胞外的漏电流IDLISID分别表示外部注入胞体和树突的刺激电流。描述树突间室膜电压VD和胞体间室膜电压VS的动力学方程如下:

img

(2.11)

img

(2.12)

式中,C表示膜电容,大小为2μF/cm2p表示形态参数,标准值为0.5,变化范围为0.01~0.99。式(2.11)等号右边的5项分别表示胞体中的突触输入电流IS、树突流向胞体的内部电流IDS、Na+电流INa、K+电流IK和漏电流ILENaEKESL分别表示Na+通道、K+通道及漏电流的平衡电势,数值分别为50mV、-100mV、-70mV。gNa=20mS/cm2gK=20mS/cm2gSL=2mS/cm2分别表示Na+电流、K+电流及漏电流的电导。m(VS)表示Na+通道开通概率的稳态值,而n表示压控K+通道的激活变量,n的动力学方程为

img

(2.13)

式中,n(VS)与τn(VS)分别表示变量n的稳态值和弛豫时间;m(VS)、n(VS)和τn(VS)都是胞体膜电压VS的函数,表达式如下:

img

(2.14)

式(2.12)等号右边的3项分别表示树突间室的突触输入电流ID、树突流向胞体的内部电流IDS和漏电流IDL。树突间室的漏电流平衡电势及电导分别为ESL=-70mV、gDL=2mS/cm2;而树突流向胞体的电流IDS如下:

img

(2.15)

式中,内连电导gc的标准值为1mS/cm2,取值范围为0.05~10mS/cm2

尽管此模型描述的树突形态与胞体结构比较简单,却可以刻画出锥体神经元在电场刺激下的空间极化效应,并且用形态参数p刻画了神经元形态。除此之外,此模型保留了PR模型中产生动作电位必需的离子种类——Na+和K+,更有利于研究电场调节神经电活动的动力学生物物理机制。

2.1.4 电缆模型

真实神经元的主要组成结构有树突、胞体和轴突。树突可以接收来自空间不同位置的突触信息,并将此信号进行整合后传入胞体,然后经过胞体的处理由轴突传递至下一神经元。在神经元信号传导方面,上述神经元模型就体现出它们的缺陷,因为它们仅刻画了神经元动力学机制及神经元形态的特性。因此电缆模型经常用来分析电信号在神经元中的传导机制,最早提出的电缆理论用于分析水下电缆的电报传输过程,随后Hodgkin和Rushton首次提出用该理论探究神经元轴突动作电位传导机制,后来Rall等人利用电缆理论分析了树突的电信号传导机制。更重要的是电缆模型可以较完美地展示真实神经元的形态,而神经元形态在电信号传导过程中又扮演着重要角色。电缆模型将神经元视为可以导电的电缆,本节结合真实神经元模型及等效电缆电路图对电缆模型进行描述。

图2.8为电缆模型,其中图2.8(b)为图2.8(a)中等效的电缆结构,图2.8(c)详细展示了图2.8(b)中片段圆柱体电缆的等效电路。在图2.8(c)中,每个主体的细胞膜电容可表示为cmΔx;图中方块表示细胞膜上的离子通道,总的电流为Im(xx,且此电流随着圆柱体轴线方向变化;Vi(x)和Ve(x)表示所在位置的细胞膜内外电势,并假设这些电势沿着圆柱体轴线方向变化,因此细胞膜的膜电位为V(x)=Vi(x)-Ve(x)。由于膜电压在不同空间位置是不同的,因此细胞膜内外电流也是浮动的,分别用Ii(x)和Ie(x)表示。图2.8中的riΔx表示第n个和第n+1个子圆柱体中心之间的细胞外等效电阻。

img

图2.8 电缆模型

在生理实验中,通常考虑用电极施加外部刺激,如细胞内电极刺激、细胞膜外电刺激,这些带电电极会影响神经元的电活动。因此将细胞内外电极刺激视为电流源输入,用符号IAiIAe表示,则电极电流在子圆柱体细胞内外空间引入的电流为IAiΔxIAeΔx,根据欧姆定律得到:

img

(2.16)

img

(2.17)

当Δx→0时,对式(2.16)和式(2.17)取极限得到:

img

(2.18)

img

(2.19)

根据电路节点定律有:

img

(2.20)

img

(2.21)

相应地,式(2.20)和式(2.21)的极限形式表达式为:

img

(2.22)

img

(2.23)

因此,联合式(2.18)、式(2.19)、式(2.22)和式(2.23)可得:

img

(2.24)

img

(2.25)

若用Iion(x, V(x), t)表示细胞膜离子通道的一般表达式,则图2.8(c)中方框#2表示的细胞膜离子通道的电流可以描述为:

img

(2.26)

联合式(2.24)和式(2.25)可以得到:

img

(2.27)

将式(2.26)代入式(2.27)得到非线性的电缆方程为:

img

(2.28)

式中,Istim(x)表示外部电极施加电刺激的等效激励项;其他两项表示空间某一位置的膜电流。Istim(x)的表达式为:

img

(2.29)

对于神经元而言,无论是细胞膜上还是细胞质中都包含多种多样的离子通道或离子,胞体主要包含Na+、K+、Cl,而树突还包含Ca2+、超极化激活阳离子等。在NEURON软件中,系统默认的离子电流为HH类离子电流,包括Na+电流、K+电流和漏电流。这类离子电流的等效电路图如图2.9(a)所示,图2.9(b)描述了图2.9(a)等效的单个反电势-电导电路。

img

图2.9 离子电流的等效电路图

由于电缆模型可以较为真实地体现神经元的形态等特征,在研究神经元形态对神经元响应的影响方面应用广泛。例如,Berzhanskaya等人[66]利用更加符合生理意义的海马体CA1区模型研究了阈下电场对θ调制γ节律的调节作用,并发现当树突直径不同时,电场对其极化作用也是不同的;伊国胜等人[67]利用海马体CA1区域神经元研究其形态对匀强电场刺激下神经元反应的影响,发现神经元的直径、长度、弯曲角度和树突分支及轴突终端通常都会不同程度地影响胞体的去极化程度,因此电缆模型在揭示神经元的生理意义方面具有非常重要的作用。

2.1.5 简化模型

在前面的章节中,我们介绍了几种具有明确生理意义的电导模型、双间室模型和电缆模型,这类电生理模型能够准确描述神经元的电生理现象,各参数具有明确的生理意义,但是模型都比较复杂,不适合大规模网络建模研究。从计算效率和数学分析的角度,能描述神经元动态的简化模型也发挥了非常重要的作用。

1.Izhikevich模型

为了理解大脑工作原理,需要将实验研究和大规模脑模型的数值研究结合起来。用于建模大规模脑网络的脉冲模型必须满足两个特征:①计算简单;②能够复现真实神经元丰富的放电模式。上述HH模型能复现各种放电现象,但结构复杂,不适合大规模网络的研究。形式简洁的Integrate-and-Fire模型计算效率很高,适合大规模网络研究,但是该模型不能复现各种丰富的脉冲放电现象及簇放电现象。基于此问题,Izhikevich结合分叉方法论,于2003年提出了一种结构简单、能复现非常多生物神经元放电现象,且适合大规模网络计算的脉冲模型——Izhikevich模型[68]。该模型结构简单,由一个二阶连续微分方程和一个离散重置部分组成,公式如下:

img

(2.30)

重置部分满足,当img时,有:

img

(2.31)

式(2.30)和式(2.31)中,Vuabcd均为无量纲变量;变量V表示神经元的膜电压;u表示膜的恢复变量,解释K+电流的激活特性和Na+电流的失活特性,并作为一个反馈变量给膜电压。当神经元的膜电压上升到30mV时,根据式(2.31)的定义,膜电压和恢复变量会被重置为一个新的值。变量I表示突触输入电流或外部直流输入电流。各参数的取值和定义按照IzhiKevich的规定如下:

(1)参数a描述恢复变量的时间尺度,其值越小,恢复过程越慢,一个典型的取值为a=0.02。

(2)参数b用来描述恢复变量对膜电压阈下振荡的敏感性,取值过大可能导致阈下振荡和低阈值脉冲放电动态,一个典型的取值为b=0.2。当b<a时,模型的静息态会发生鞍节点分岔;当b>a时,则对应Hopf分岔。

(3)参数c刻画膜电压超过阈值被重置后的状态,一个典型的取值为c=-65。

(4)参数d描述恢复变量被重置后的状态,一个典型的取值为d=2。

实际上,除上述典型取值以外,这些参数还可以取不同的值,不同的参数取值对应不同的放电模式,Liu等人在最近的论文中分析了在不同放电模式下各种参数的取值范围[69]。根据不同取值的组合,Izhikevich模型能复现几类主要的神经元放电模式,如图2.10所示。其中,RS、IB、CH表示皮层兴奋性神经元;FS、LTS表示皮层抑制性中间神经元,对每类神经元来说,接收的外部输入为I=10mA的直流阶跃输入。

我们以图2.10中CH神经元展现的簇放电为例,给出实现该放电模式的Matlab代码:

img
img

图2.10 Izhikevich模型复现的各种放电模式[68]

2.跨临界混合模型

经典的HH模型及所有基于HH模型的简化模型主要描述Na+和K+在动作电位产生中的作用机制,Na+电流是一个快速去极化电流,K+电流是一个慢速超极化电流。实验研究发现Ca2+电流在多巴胺能神经元的放电机制中发挥了重要的作用,在与帕金森病相关的STN核团同步中也起到了调制作用。因此,Drion等人[70]采用之前的模型简化方法,将Ca2+动态添加到模型中,通过化简得到一个结构类似Izhikevich模型的跨临界混合模型。这个模型具有新的相平面特性,丰富了Izhikevich模型的建模能力,最重要的是,该模型中有一个关键参数能建模生物神经元的Ca2+动态。在低Ca2+动态下,该模型与早期的简化模型动态相似,但一般的简化模型不能描述高Ca2+的动态,而这个跨临界混合模型能准确描述实验中观察到的高Ca2+时的电生理现象[71],被认为是建模基底核运动回路的理想模型。

描述跨临界混合模型动态的微分方程如下:

img

(2.32)

与Izhikevich模型相似,该跨临界混合模型也具有一个重置函数,当img时,有:

img

(2.33)

式(2.32)和式(2.33)中,v表示膜电压的快速变量;w表示慢速恢复变量;z表示慢速适应变量,用于描述细胞内Ca2+变化的影响及相关的Ca2+激活泵电流变化情况。模型中参数的数值为:a=0.25、b=-1.5、c=-10、d=20、ε=1、vth=40、εz=0.1、dz=150。参数w0表示Ca2+电导,w0>0表示低Ca2+电导状态,w0<0表示高Ca2+电导状态。仿真研究表明,在不同的Ca2+浓度下,神经元具有不同的放电模式。

图2.11和图2.12分别描述了跨临界混合模型在高Ca2+和低Ca2+模式下的放电现象,可以明显地看出,高Ca2+状态会出现簇放电现象,且放电模式具有脉冲延迟(Spike Latency)、去极化后电势(ADP)及Plateau振荡等多巴胺能神经元在高Ca2+状态才具有的典型电生理现象;而在低Ca2+状态下,不会出现上述3类放电现象。为了分析高、低Ca2+状态时放电模式的不同,我们绘制了跨临界混合模型的膜电压相平面图,如图2.13所示。

img

图2.11 高Ca2+电导情况下(w0=-3),神经元在不同参数下的放电模式

img

图2.12 低Ca2+电导情况下(w0=4),神经元在不同参数下的放电模式

img

图2.13 跨临界混合模型在低Ca2+和高Ca2+状态时的膜电压及相平面图

选取参数使得神经元在高、低Ca2+模式下均呈现脉冲放电行为,根据其相平面图我们可以发现,在低Ca2+模式下,w零线与V零线的上支相交,当神经元发放一个动作电位后,由于相交点的吸引力作用,会很快开始下一个脉冲放电。而在高Ca2+模式下,w零线与V零线的下支相交,当神经元发放一个动作电位后,由于受到下方相交点的吸引力作用,会沿着w零线运行一段,然后越过相交点进入下一个脉冲放电,因此产生了去极化后电势的现象。不同Ca2+电导模式下的放电差异也表现在模型的频率-电流曲线上,如图2.14所示。图2.14(a)为高Ca2+电导下的频率-电流图,图2.14(b)为低Ca2+电导下的频率-电流图。可以发现,在高Ca2+模式下,神经元可以以任意低的频率开始放电,因此属于Ⅰ型兴奋性神经元;在低Ca2+模式下,当电流较小时,神经元处于静息态,只有当电流足够大时,神经元才能以某一频率开始放电,即属于Ⅱ型兴奋性神经元。

img

图2.14 跨临界混合模型频率−电流曲线

3.Fitzhugh-Nagumo模型

Hodgkin和Huxley的首创工作及后续的相关研究在神经脉冲传导方面建立了很好的数学模型。这些模型大都是耦合的微分方程,为了更方便地定性描述神经元行为,FitzHugh和Nagumo将四维HH模型简化,从而提出二维Fitzhugh-Nagumo模型[72,73]。该模型忽略了Na+和K+流动的电化学特性传播及激活造成的本质数学性质,其表达式如下:

img

(2.34)

式中,J表示外部输入电流;y表示电压变量,该变量具有立方非线性动力学特性,可以通过正反馈进行自我激活再生;x表示恢复变量,此变量具有线性动力学特性,允许缓慢的负反馈;abc为常数,其约束条件为:

img

(2.35)

二维Fitzhugh-Nagumo模型的等价电路图如图2.15所示。

img

图2.15 二维Fitzhugh-Nagumo模型等价电路图[72]

尽管此模型比较简单,但是仍然可以用来对重要的与兴奋性和放电机制相关的生物现象进行数学解释。二维Fitzhugh-Nagumo模型的相平面图如图2.16所示。

img

图2.16 二维Fitzhugh-Nagumo模型的相平面图[80]

其中,直线区间的x零线是式(2.34)中恢复变量x的微分方程为0时得到的,并且参数取值为a=0.7、b=0.8、c=3、z=0。可以看出,电压变化比恢复变量变化缓慢,除了接近电压零线的位置。如果电压为常数,也就是电压一阶微分方程为零,即电压不是变化的,则相应相平面的水平线可以视为仅有变量x的简化系统的相线。此水平线穿过静息态P,如图2.16中标记的一样,此直线有3个与x零线相交的奇点。其中一个奇点是不稳定的,表示阈值;在另外两个奇点中,左侧的奇点为稳定的兴奋点,右侧的奇点为稳定的静止点。从P点的相点移动到不稳定阈值点的左侧时会在简化的系统中产生兴奋,并且相点接近兴奋奇点。当恢复变量x向左侧移动时,电压变量遵循式(2.34)中的电压一阶微分方程缓慢变化,导致相线向上移动,直到与兴奋奇点和阈值奇点相遇、消失。在单一变量x的简化系统中,静止点在右支。最终电压变量缓慢减小,相平面中的相点接近静止点P

在图2.16中穿过P点的水平点线也是由z中瞬时脉冲替换的静息相点轨迹。这样的脉冲没有直接改变膜电压。在x−y相平面中,与x轴相比,阈值现象不会在单一奇点处产生,而是准类(QTP)。此类允许结余全或无之间的任何响应,只在有准确刺激的幅值时出现。计算机求解中出现了从右侧传递到左侧的邻近路径。用近似的计算,遵循分离线不会偏离分界线很远进入自由区。仿真中的噪声也会使相点向左或向右远离分割线,这意味着会产生全或无放电。当相点从右向左穿过分离线时会产生动作电位,从左向右跨越分割线则会消失。但是当分割路径逐渐消失在自由区的上方时,阈值现象会变得不快速,可见Fitzhugh-Nagumo模型与HH模型类似。

除此之外,Fitzhugh-Nagumo模型也可以阐释激发块现象。例如,当刺激电流的振幅增加时重复脉冲中止;当刺激电流较弱或者为0时,电压零线与恢复变量零线的交点在电压零线的左支上,此时模型处于静息态。增加刺激电流会移动零线,从而使平衡点上移至电压零线的中间部分,模型展现出周期簇放电现象。继续增加刺激电流幅值,会移动平衡点至电压零线的右支,此时振荡受阻。精确的数学机制涉及极限环吸引子。

Fitzhugh-Nagumo模型也可以阐述回弹脉冲,当刺激电流为负时,静息态左移。当系统从超极化释放时,轨迹从远远小于静息态的一点起始,形成一个较大的偏移,激发暂时放电然后恢复静息态。Fitzhugh-Nagumo模型在揭示HH模型的放电调整动力学机制方面也有贡献。当刺激电流强度缓慢增加时,神经元仍然保持静息态,并且Fitzhugh-Nagumo模型的静息平衡点缓慢向右移动,系统的状态逐渐出现不放电现象。相反地,当刺激迅速增加时,即使增加得很小,轨迹不能直接到达静息态,而是诱发暂时放电。此现象与回弹响应类似,更重要的是,Fitzhugh-Nagumo模型表达式成为反应扩散系统(Reaction-Difusion Systems)的最优方程,改善后的方程如下:

img

(2.36)

式中,f(y)表示y的三次立方项;I表示输入电流;Vxx表示空间变量x的二阶微分。该方程易于分析,可以得到在不借助计算机模拟时传输脉冲的许多重要特性,可以模拟激发介质波形的传输,如心脏组织或神经纤维。

4.Hindmarsh-Rose模型

蜗牛等无脊椎动物的脑部神经元初始时是静止的,但是可以由短暂的电流脉冲引发去极化,从而产生比刺激电流大得多的簇放电现象。Thompon及其合作者在多种软体动物的破裂神经元中也观察到相似类型的反应,并且发现需要超极化才能阻止这种簇放电。当这些持续的超极化细胞由短暂电流脉冲去极化时,会产生一个动作电位,并且跟随一个缓慢的去极化后电势[74]。为了对这一现象进行解释,Hindmarsh和Rose建立了一个簇放电神经元的多项式模型[75],即Hindmarsh-Rose模型,该模型的表达式如下:

img

(2.37)

式中,x表示膜动作电位;y表示恢复变量;z表示慢适应电流;abcdsx1都是常数。(x1, y1)是Hindmarsh-Rose模型在没有外部刺激电流时最左侧的坐标。Hindmarsh-Rose模型对很小的去极化电流脉冲的响应依赖给定的常数rs。当r=0.001、s=1、a=1、b=3、c=1、d=5时,此模型会产生独立的动作电位簇,如图2.17(a)所示。该现象与实验中观察到的现象相似,如图2.17(b)所示。当r=0.001、s=4时,可以得到图2.17(c)中的现象,与Thompson等人观察到的现象[74]相似。类似地,记录到的簇放电趋向于一个调整之前的初始加速度[见图2.17(b)]。这些稳定下降后极化或放电频率现象在模型中没有发现,其中一个有趣的特色是簇放电之后模型缓慢超极化到某个恢复变量的值,因此膜电压比静息态电势的恢复变量(x1)更小。当完全超极化时,我们发现恢复变量很缓慢地恢复到初始值x1,如图2.18(a)所示。在图2.18中a=1、b=3、c=1、d=5、r=0.001、s=4,但是对于这3幅图的刺激电流的幅值分别为I=0.4mA[见图2.18(a)]、I=2mA[见图2.18(b)]、I=4mA[见图2.18(c)]。其中,图2.18(b)在刺激开始后700ms开始记录,图2.18(c)在刺激开始后1000ms开始记录。

观察x−y相平面,如图2.19可以分析这些现象出现的原因。严格地讲,由于Hindmarsh-Rose模型中存在3个变量,即xyz,因此应该观察三维相平面。然而当r=0.001时,变量z相较于变量xy变化比较缓慢。因此在式(2.37)中xy的一阶微分函数内,我们将z视为缓慢变化参数。考虑神经元模型在x1的初始状态,每次动作电位产生时调整电流会增加,这也会导致x零线向上移动至连续周期。放电频率会下降,这是因为相平面左侧的狭窄通道更加狭窄。当x零线向上移动时平衡点也会发生变化,鞍点向上移动,同时最左侧的平衡点向下移动。最终鞍点会向上移动直至极限环轨迹穿越鞍点分界线。放电终止,相点进入低于鞍点的狭窄通道并缓慢向最左侧的平衡点移动。这样的缓慢移动导致缓慢超极化波形的出现,最左侧的平衡点已经向下移动,它接近一个低于x1的值。当调整电流恢复到初始状态时,当最左侧的平衡点移动到原始位置时,x的值会重新接近x1

img

图2.17 模型的适应性[75]

img

图2.18 簇放电

我们对Hindmarsh-Rose模型的簇放电现象进行阐述,其中刺激电流保持恒定的正数。在不同刺激电流下的数值仿真结果如图2.18所示,其他参数及记录方式如上所述。当刺激电流幅值较低时(I=0.4mA),模型产生单个簇放电,随后是后超极化波形,这意味着膜电压缓慢恢复到初始状态。当刺激电流增加时(I=2mA),模型出现长时间簇放电,调整和终止产生周期簇放电的情况如图2.18(b)所示。当刺激电流幅值更大时(I=4mA),出现持续高频放电,放电频率在电流刺激起始阶段逐渐下降,直至稳定状态,此时神经元重复放电。神经元模型xy相平面如图2.19所示。簇放电机制可以用图2.20阐释,图2.20(a)为簇放电周期的机制解释;图2.20(b)和图2.20(c)为α时刻和β时刻对应的相平面图;图2.20(d)和图2.20(e)刻画了γ时刻和ε时刻的x−y相平面图。当处于α时刻时,调整电流刚刚上升,当刺激电流将3个平衡点改变为一个不稳定的平衡点通过时,电压零线向下移动,如图2.20(b)所示。相应地,由于极限环的出现,模型出现粗放电行为。当调整电流上升,电压零线向上移动时,一个鞍点随之产生,并随着调整电流的持续增加而缓慢向上移动;与此同时,调整电流持续上升。最终,相点穿过鞍点分割线消失,进入低于鞍点的狭窄通道,簇放电现象消失,如图2.20(c)所示。在β时刻,相点在狭窄通道下方向最左侧平衡点移动,如图2.20(c)所示。当后簇放电超极化完成时,相点到达最左侧的平衡点,调整电流开始下降,向电压零线的下方移动;与此同时,最左侧的平衡点和鞍点互相靠近,如图2.20(d)所示,更小的平衡点缓慢向上运动,造成在γ时刻的簇间起搏器去极化。在δ时刻两个平衡点互相消融,最后模型在ε时刻开始返回一个平衡点,如图2.20(f)和图2.20(b)所示,产生另一个簇。

img

图2.19 神经元模型x−y相平面[75]

注:a=1、b=3、c=1、d=5,其中ABC分别是稳点节点、不稳定鞍点、不稳定螺旋,相应的电压电流趋势在相平面上方展示。

img

图2.20 簇放电机制[75]

除了簇放电,Hindmarsh-Rose模型还能出现随机簇放电现象及后抑制回弹现象。更重要的是,该模型可以用来解释软体动物神经元簇放电实验现象[76],如在内部簇间隔中,上升的去极化会激活阈下内流电流,该电流主要由Ca2+流动形成,从而激发簇放电。在簇放电过程中,K+缓慢地去极化及动作电位造成Ca2+内流,从而激活缓慢外流Ca2+依赖K+电流。此电流使后簇超极化(Post-Burst Hyperpolarization)循环终止。因而Ca2+依赖K+电流下降,细胞逐渐去极化,开始另外的循环。在这一过程中,Ca2+依赖K+电流相当于Hindmarsh-Rose模型中的缓慢适应电流。除此之外,此模型也可以预测兔子第5层皮层锥体神经元在随机输入信号的响应,并且准确度高于其他已知的神经元模型,当模型参数调整到适合阈下测量的现象时,此模型仍然可以很好地表现电特性[88]