哺乳动物皮层神经元模型的放电模式研究

程元芬，刘深泉

华南理工大学数学学院，广东广州510640

前言

Hodgkin等［1］在1952年通过电生理实验研究发现枪乌贼轴突动作电位的产生与4种非线性的、时间依赖变量的相互作用有关，这4个变量分别是膜电位V、钠离子激活变量m、钠离子失活变量h以及钾离子激活变量n。以Hodgkin-Huxley模型为基准，Mccormick等［2］通过实验研究发现，影响人类及哺乳动物新皮层神经元动作电位产生的因素更加复杂，至少包括12种不同离子电流的相互作用。除此之外，还包括低阈值的钙离子电流It、慢变后超极化钾离子电流IAHP等多种电流的相互作用，使得皮层神经元的发放模式呈现出多种多样的形态。这些离子通道相互作用产生的复杂性使得对皮层神经元的分析变得困难，所以人们开始寻求用简单的数学模型来对其进行研究。Morris等［3］在1981年提出著名的Morris-Lecar模型，该模型在不同的参数下能得到不同的神经元放电模式；Hindmarsh等［4］在1984年提出一个简单且易于分析的三维神经元模型，但这两个模型的缺点在于模型并不遵守欧姆定律；后来，Wilson［5］提出一个四维神经元模型，它虽然改进了Hindmarsh-Rose模型中不遵守欧姆定律的这一缺点，但由于该模型是四维的，这给系统的分析带来了不便。

不同的放电模式对应着神经元不同的信息编码方式，利用动力系统分岔理论和快慢动力学方法来研究神经元模型不同类型的放电活动，被认为是行之有效的方法［6-10］，推动了神经动力学的发展。近年来，利用快慢动力学的方法，杨卓琴等［11］对Chay神经元模型的簇放电模式进行了研究；Wang等［12-13］研究了胰腺β细胞神经元模型及呼吸神经元模型的簇放电模式；Zhan等［14-15］对Purkinje神经元模型及垂体模型的放电模式进行了研究。可见快慢动力学对于时间尺度差异很大的系统的分析是非常有效的，本研究也将采用该方法来对模型的簇放电模式进行研究。

为了保持三维系统易于分析的优点，又能遵守生物的生理特性，Zhao等［16］提出一个哺乳动物新皮层神经元模型，该模型是将Wilson模型中的快变量与Hindmarsh-Rose模型中的慢反馈变量相结合得到的一个三维模型，所得模型既保持了三维模型易于分析的优点，又在一定程度上遵守了生物的生理特性，同时具有丰富的发放现象。本研究以该模型为基础，分析该系统的不同簇放电模式；利用快慢动力学分析方法，作出快变子系统关于慢变量的平衡点分岔曲线，分析得出快变子系统的簇发放类型；将该模型与Morris-Lecar神经元模型进行电耦合，分析耦合强度及交流电频率对神经元动力学行为的影响。

1 模型描述

为了保持三维神经元模型易于分析的特性，同时在一定程度上遵守欧姆定律，文献［16］提出了如下的三维神经元模型：

其中，

模型（1）中的第一个方程描述的是膜电压的变化。其中，V为膜电压；Cm为膜电容；g(V)为钠离子激活函数；gR为钾离子电导；H表示调节电流；gH代表调节电流所对应的电导；VNa和VK分别为钠离子和钾离子的平衡电压；Iapp为外界刺激电流；R为钾离子激活变量；R∞为激活变量平衡态；τR为R松弛时间变量；Vh表示反转电压；τH为H松弛时间变量。详细描述见参考文献［16］。

模型部分参数取值为：Cm=30 μF/m2，VNa=0.048 V，VK=-0.095 V，V3=-0.038 V，VH=-0.075 04 V，τR=0.56 ms，τH=100 ms，v0=178.1 Am-2V-1，v1=4 758 Am-2V-1，v2=33 800 Am-2V-1，r1=12.9 V-1，r2=330 V-1。

文献［16］通过数值模拟说明了该模型在不同的参数下具有丰富的发放现象。本研究利用龙格-库塔算法，结合MATCONT软件，得出快变子系统在不同参数下的几种分岔类型。同时理论计算了Hopf分岔点的一阶Lyapunov系数，从而确定Hopf分岔的方向。最后，将该模型与Morris-Lecar神经元模型进行电耦合，研究通过耦合得到的新模型的发放模式。

2 快变子系统分岔曲线上的簇放电类型

在神经元模型（1）中，由于H的变化要比R慢很多，因而本研究将H看作慢变量。由此，模型（1）可以看成是由

和系统

构成，将慢变量H看作快变子系统［式（3）］的分岔参数，从而对系统［式（1）］进行快慢动力学分析。

2.1 经由“fold/homoclinic”滞后环的“fold/homoclinic”型簇放电

当gR=260 Am-2V-1，gH=13 Am-2V-1时，神经元模型（1）中的膜电位时间历程图如图1a所示，运用快慢动力学和MATCONT软件来研究该系统的簇放电模式，得到图1b。

分析图1b可知，快变子系统（2）的平衡点曲线为一条Z型曲线，由3部分构成。其中稳定的焦点和不稳定焦点以Hopf分岔点为分界点构成了曲线的上支，即随着分岔参数H的增大，曲线上的稳定焦点经由Hopf分岔点变成了不稳定焦点，同时出现稳定的极限环；鞍点构成了该曲线的中支；曲线的下支由稳定的结点构成。将图1a中的簇发放轨线图叠加在快变子系统平衡点曲线上，分析在该参数条件下的簇发放类型。

随着控制变量H的减小，平衡点分岔曲线下支相应于静息状态的稳定结点消失，向上转迁到相应于放电状态的稳定极限环周围。经由Hopf分岔点分支出的稳定极限环随着分岔参数H的增大，逐渐逼近由鞍点组成的中支，最后碰到鞍点变成鞍点同宿轨，经由鞍点同宿轨线回到相应于静息状态的稳定结点处，放电状态结束。因此，静息状态与发放状态相互转迁的放电类型为“fold/homoclinic”型簇放电。曲线的下状态经由鞍结点跃迁至曲线的上状态，上状态经由鞍点同宿轨回到下状态。因此滞后环产生的簇放电为“fold/homoclinic”型簇放电，根据文献［10］提出的放电类型分类方法，此簇放电的模式为经由“fold/homoclinic”滞后环的“fold/homoclinic”型簇放电。

图1 当gR=260 Am-2V-1，gH=13 Am-2V-1时，产生经由“fold/homoclinic”滞后环的“fold/homoclinic”型簇放电Fig 1 Fast-slow dynamics of"fold/homoclinic"bursting via"fold/homoclinic"hysteresis loop whengR=260 Am-2V-1,gH=13 Am-2V-1

2.2 经由“fold/homoclinic”滞后环的“subHopf/homoclinic”型簇放电

当gR=105 Am-2V-1，gH=40 Am-2V-1时，系统（1）［式（1）］有图2a所示的膜电位时间历程图。将图2a的簇发放轨线图附于快变子系统平衡点分岔曲线上，得到图2b。根据数值计算结果，Z型平衡点曲线上支存在一个Hopf分岔点，分岔点左边部分为稳定焦点构成，右侧由不稳定焦点构成。分岔曲线的中支和下支分别由鞍点和稳定结点构成。

图2 当gR=105 Am-2V-1，gH=40 Am-2V-1时，产生经由“fold/homoclinic”滞后环的“subHopf/homoclinic”型簇放电Fig 2 Fast-slow dynamics of"subHopf/homoclinic"via"fold/homoclinic"hysteresis loop whengR=105 Am-2V-1,gH=40 Am-2V-1

为确定该参数下的Hopf分岔方向，本研究将计算该Hopf分岔的一阶Lyapunov系数。当gR=105 Am-2V-1，H=1.472 899 Am-2时，重写该系统的快变子系统：

其中，

快变子系统的Jacobian矩阵为：

其中，

当H=1.472 899 Am-2时，快变子系统的平衡点为(- 0.017 212,0.710 583)，可得该平衡点处的Jacobian矩阵为：

经计算，该矩阵的特征值为λ1,2=±wi,w=3.123 0，由此可知，在该平衡点处产生了Hopf分岔。另外，λ对应的特征向量为q=(0.037 565+0.065 698i,1)Τ，另取一个向量p使得p满足AΤp=-wip，且p,q=1，计算得到：

考虑系统

其中，F(x)光滑。F关于x的Taylor展开F(x)=Ο(||x2||)至少从二次方开始，A是平衡点处的Jacobian矩阵，F(x)可写成如下形式：

这里B(x,y)和C(x,y,z)是x,y,z∈R2的对称多重线性向量函数。它们的表达式如下所示：

其中，ξ=(ξ1,ξ2)为系统的平衡点。

将快变子系统写成如下形式：

其中，

因此，有：

可得一阶Lyapunov系数为：

因此该Hopf分岔方向是亚临界的。

根据图2b，快变子系统平衡点曲线有两处静息态，分别对应于平衡点曲线的下支及曲线上支Hopf分岔点的左侧。随着控制变量H的减小，下支静息态经由鞍结点转迁到曲线上支的静息状态，在上支静息态周围发生小幅衰减振荡后收敛于稳定焦点。而后随着H的增大，静息态经由Hopf分岔消失，转迁为发放状态，经过几个峰值不同的发放后重新转迁回静息状态。

此外，系统从下状态经由鞍结点跃迁到上状态，经由鞍点同宿轨线从上状态跃迁到下状态。因此，在该参数下的簇放电模式为经由“fold/homoclinic”滞后环“subHopf/homoclinic”型簇放电。

2.3 “fold/fold”点-点滞后环型簇放电

当gR=16.5 Am-2V-1，gH=120 Am-2V-1时，系统（1）产生如图3a的发放模式。快变子系统的平衡点曲线呈现Z型，曲线的上支、中支、下支分别由稳定焦点、鞍点、稳定结点构成。将图3a的发放图附于平衡点曲线上，得到快变子系统的快慢动力学分析图3b。由于在该参数下，快变子系统分岔曲线不存在Hopf分岔点，因此不存在相应于放电状态的稳定极限环，故不需要讨论静息状态和放电状态相互转迁的分岔类型，只需要讨论与滞后环的产生有关的分岔。从图中可以看出下状态经由鞍结点转迁到上状态，上状态又经由鞍结点转迁到下状态。因此，在该参数下的系统簇放电模式为“fold/fold”点-点滞后环型簇放电。

3 耦合神经元模型的发放特性

以该三维模型和Morris-Lecar模型为基础，将这两个具有丰富发放模式的神经元模型进行电耦合，得到一个新的六维神经元模型，其动力学方程如下：

图3 当gR=16.5 Am-2V-1，gH=120 Am-2V-1时，产生“fold/fold”点-点滞后环型簇放电Fig 3 Fast-slow dynamics of“fold/fold”hysteresis loop bursting of point-point type forgR=16.5 Am-2V-1,gH=120 Am-2V-1

其中，V1表示Morris-Lecar模型神经元膜电位，V2为三维哺乳动物神经元膜电位，gc神经元耦合强度，gL为泄露离子通道电导，Iapp1和Iapp2分别为加入到两个神经元中直流电，Iext1和Iext2为添加到两个神经元中的交流电刺激，Cn表示ML神经元模型中的膜电容，gK和gca表示钾离子和钙离子的电导，Vca为钙离子平衡电压，w表示门控变量。

此外：

参数取值为VL=-0.5 mV，VK=-1.1 mV，VCa=1.0 mV，gL=0.5 ms/cm2，gk=2.0 ms/cm2，gCa=1.2 ms/cm2，v11=-0.01mV，v22=-0.15mV，v33=-0.10mV，v44=-0.05mV。

3.1 耦合强度对耦合神经元膜电位的影响

固定Iapp1=0.35 mA/cm2，Iapp2=0.30 mA/cm2，不添加交流电刺激，在此基础上研究耦合强度变化对MORRIS-LECAR神经元膜电位发放模式的影响。图4a～d为不同耦合强度所对应的膜电位时间历程图。当耦合强度为gc=ms/cm2时，耦合模型中的MORRIS-LECAR膜电位只有一种类型的簇；当gc=0.1 ms/cm2时，膜电位的一个周期内包含两个簇，其中左边的簇包含12个峰数量，右边的簇中包含2个峰数量；将gc增大到0.2 ms/cm2时，构成一个周期的两个簇中峰的数量都为2，但左边簇的峰值明显比右边簇的峰值大；继续增大gc到0.3 ms/cm2，结果显示，膜电位周期性的出现一种簇。显然，耦合强度的变化对耦合神经元膜电位有明显的影响。

3.2 交流电刺激对耦合神经元膜电位的影响

固定Iapp1=0 mA/cm2，Iapp2=0.3 mA/cm2，Iext2=0 mA/cm2，gc=0.2 ms/cm2，在模型（4）的第一个方程中加入交流电刺激Iext1=0.5 sin(ωt)，通过改变交流电频率，由此来研究交流电频率变化对耦合神经元中MORRIS-LECAR膜电位发放模式的影响。

图5a～d为不同交流电频率所对应的膜电位时间历程图。当ω=0 rad/ms时，膜电位的一个周期包含两个不同的簇，左边簇的峰数量为3，右边簇的峰数量为2；当ω=0.028 rad/ms时，周期内两种不同簇的间距减小，并且两种簇中所含峰的数量都明显增加；当ω增大到0.086 rad/ms时，膜电位时间历程图中只出现了一种簇；继续增大ω到0.14 rad/ms时，该历程图一个周期内出现了4种不同的簇。显然，交流电刺激频率的改变对膜电位有明显的影响。

4 小结

细胞膜上的离子通道是影响神经元放电活动的重要因素，影响人类哺乳动物新皮层神经元膜电压变化的因素有很多，对其的分析也相当复杂，文献［16］所提出的简化模型使得对哺乳动物新皮层神经元放电模式的分析变得简单。本研究采用快慢动力学分岔分析方法，对gR、gH取不同的值，以研究该三维哺乳动物新皮层神经元模型的簇放电模式。对于该系统的快变子系统平衡点随控制变量H的变化所形成的分岔曲线，主要出现了以下3种簇放电模式：（1）当gR=260 Am-2V-1，gH=13 Am-2V-1时，系统的分岔表现为经由“fold/homoclinic”滞后环的“fold/homoclinic”型簇放电；（2）当gR=105 Am-2V-1，gH=40 Am-2V-1时，系统的分岔类型为经由“fold/homoclinic”滞后环的“subHopf/homoclinic”型簇放电。同时，通过理论计算可知，快变子系统分岔曲线上Hopf分岔点对应的一阶Lyapunov系数为0.443 025，由此可判断该Hopf分岔为亚临界Hopf分岔；（3）当gR=16.5 Am-2V-1，gH=120 Am-2V-1时，系统的簇放电模式“fold/fold”点-点为滞后环型簇放电。将该模型与MORRIS-LECAR神经元模型进行电耦合，得到一个新的六维神经元模型。通过改变耦合强度和交流电频率，膜电位出现丰富的变化。结果表明，耦合强度和交流电频率的变化对神经元膜电位产生非常明显的影响。

图4 当Iapp1=0.35 ms/cm2，Iapp2=0.3 ms/cm2时，不同的耦合强度所对应的膜电压图Fig 4 Membrane potential diagrams in different coupling strengths,withIapp1=0.35 ms/cm2,Iapp2=0.3 ms/cm2

图5 当Iext1=0.5sin(ωt)时，不同的交流电频率所对应的膜电位图Fig 5 Membrane potential diagrams at different AC frequencies,withIext1=0.5sin(ωt)