电磁感应下神经元模型的放电模式及分岔分析

赵晴晴，刘深泉**，孟 盼

(1.华南理工大学 数学学院，广东 广州 510641；2.广东药科大学 医药信息工程学院，广东 广州 510006)

生物神经系统是由大量神经元连接而成的多级信息网络，其基本结构是神经元，神经信息的表达主要依赖于神经元电生理活动的实现[1].神经元模型的建立为研究神经元的放电模式提供了很好的途径.1952年Hodgkin 和Huxley 建立的Hodgkin-Huxley(HH)离子通道模型是神经元理论研究的基础[2].在HH 模型的基础上先后建立了FitzHugh-Nagumo 模型[3-4]、Morris-Lecar 模型[5]、Hindmarsh-Rose 模型[6]、Chay 模型等[7].许多神经科学家利用分岔理论和快慢动力学方法对神经元模型的放电模式和动力学特性进行了研究[8-12].

传统神经元模型的研究主要考虑离子通道电流、外界刺激对细胞膜电位的影响，实际上细胞内外带电离子的分布是不均匀的，根据法拉第电磁感应定律，细胞内外的带电离子穿越细胞膜会产生电磁感应效应，进而改变神经元的电生理活动.因此，在研究神经元模型的放电活动时应考虑穿过细胞膜的磁通和电磁效应.近年来，许多神经科学家开始研究电磁感应对神经元和神经网络的影响.Lyu 等[13]在Hindmarsh-Rose 神经元模型中引入磁通变量和膜电位反馈电流研究其电活动.Ma 等[14]基于FitzHugh-Nagumo 神经元模型，研究了外部电磁辐射对心肌细胞电活动的影响.Zhan 等[15]在Morris-Lecar 模型中同时考虑了电磁感应和高斯白噪声，发现电磁感应可以将神经元膜电位从放电状态转变为静止状态.Duan 等[16]对呼吸神经元模型在电磁感应下的动力学性质进行了研究.乔帅等[17]研究了磁通 e-HR 神经元模型的放电特性和分岔现象，发现该系统存在隐藏放电行为.李宁等[18]研究了胰腺 β细胞在电磁感应下的放电模式和同步现象.Tabi 等[19]构建了电磁感应影响下的Hindmarsh-Rose 神经网络，研究了磁耦合与电耦合对神经网络的影响以及神经网络同步的原因.

簇模式是神经元中主要的放电模式，不同簇模式间的差异性表明了神经信息的多样性，掌握神经元的各种放电节律模式的产生机理以及不同放电模式之间的转迁过程对理解神经系统中信息的编码和传递是极其重要的.研究表明，簇放电的产生机理可以定性地归因于系统内部不同时间尺度分隔而成的2 个子系统相互作用的结果：一个是快子系统，表现出静息或周期性振动状态；另一个是慢子系统，负责调控静息态和活跃态的转迁.因此可以将快变量与慢变量分离开来，分析快子系统的发放行为，这种方法被称为快慢动力学分析.Rinzel 等[20-21]首次将快慢动力学分析应用到对β 细胞的研究中，随后系统地分析了神经元簇放电模式的产生机理.利用该方法，科学家对抛物型、椭圆型、圆锥型等簇放电模式的产生机理进行了分析[22-24].但是这种分类方法只能将簇放电模式大致地分为不同的类型，却很难分辨出其他具有更细微区别的簇放电模式.因此Izhikevich[25]在前人的研究成果之上考虑了与簇放电模式有关的分岔点，对簇放电模式进行了更加详细的分类.

本文以三维的类Hindmarsh-Rose 神经元模型[26]为研究对象，利用分岔理论和快慢动力学方法探讨了电磁感应作用下新四维神经元模型的分岔现象及放电模式.为了防止混淆特别说明文献[13]是在传统的Hindmarsh-Rose 模型的基础上进行的研究，而本文是在Tsaneva-Atanasova 基于Hindmarsh-Rose 模型的结构构建的类Hindmarsh-Rose 模型的基础上进行的研究.本文安排如下：第1 节简单介绍了模型；第2 节利用余维1 分岔图、峰峰间距分岔图探究了不同参数下系统的动力学性质，着重计算了用于判断Hopf 分岔稳定性的 L yapunov系 数，通过 Lyapunov指数谱确定系统的混沌区域；在第3 节中利用快慢动力学分析得到了4 种簇放电模式；第4 节中利用快变子系统的双参数分岔理论，分析了不同簇放电模式间的转迁过程，接着通过全系统的余维2 分岔得出神经元处于放电状态的参数区域，并给出余维2 分岔点附近的拓扑规范形；在第5 节给出结论.

1 模型介绍

为了研究内分泌细胞模型中的两类簇放电模式，Tsaneva-Atanasova[26]基于Hindmarsh-Rose 模型的结构构建了三维的类Hindmarsh-Rose 模型.我们在该三维模型的基础上引入磁通变量和感应电流（磁通量和磁场变化引起的电流），得到以下新四维模型：

式中：变量x、y、z、φ分别表示膜电位、电压依赖的钾离子通道门控开通的概率、细胞内钙离子浓度、磁通量， ϕ、ε是控制时间尺度分离的时间常数，因为钙离子浓度变化的时间尺度比其他变量小几个数量级，所以本文取 ϕ =1,ε=0.02.-s(-ax3+x2)表 示内向钙离子电流， -y表 示外向电压依赖的钾离子电流， -bz表示外向钙敏感的钾离子电流， -k0ρ(φ)x表示感应电流，k0为 感应电流的反馈系数.ρ(φ)=α+3βφ2为磁通控制记忆电阻器的记忆电导，用来描述神经元膜电位与磁通的耦合关系，其中 α和 β 为固定参数值 ( α=1, β=0.001)，k1x和k2φ分别表示膜电位上电磁感应变化的磁通量和泄露磁通量，具体描述见文献[13].模型中部分参数取值为：a=0.5,b=1,a1=-0.1,b1=-0.05,k=0.2,k1=1,k2=3.

本文使用Matlab 和Matcont 软件包对模型进行求解和图像处理，算法采用步长为0.01 的四阶龙格-库塔算法.

2 单参数分岔分析

利用余维1 分岔和峰峰间距(ISIs)分岔图研究单参数对系统动力学的影响.

首先分析参数s对系统动力学的影响.利用Matcont 计算得到系统的余维1 分岔点，如图1 所示，其中蓝色和黑色曲线分别代表稳定和不稳定的平衡点.这条分岔曲线上有1个中性鞍点H和1 个Hopf 分岔点H1，此时s=-1.402 134.当s＜-1.402 134时 ，系统(1)只有1 个不稳定的平衡点，神经元处于放电状态.当s=-1.402 134时 ，系统经由Hopf 分岔点H1结 束放电状态，当s＞-1.402 134时系统有1 个稳定的平衡点，此时神经元处于静息状态.神经元静息状态与放电状态的转迁和平衡点曲线上的分岔点有关.

图1 参数 s 的余维1 分岔图Fig.1 One-parameter bifurcation versus s

现在验证图1 中分岔点的动力学性质，通过计算Hopf 分岔的第1 个 Lyapunov 系 数，判断 Hopf分岔是亚临界还是超临界.当第1 个 L yapunov系数为正(或负)时，Hopf 分岔为亚临界(或超临界).

接下来计算Hopf 分岔点H1的第1 个Lyapunov 系数.为了求得平衡点，我们将系统(1)改写为：

式中：F1=-s(-ax3+x2)-y-bz-k0ρ(φ)x，F2=ϕ(x2-y)，F3=ε(sa1x+b1-kz)，F4=k1x-k2φ.式（2）中各参数的定义与第1 节相同.

首先需要计算系统(2)的雅可比矩阵，雅可比矩阵可以表示为：

当s=-1.402 134时 ，系统(2)的平衡点为 (0.379 145, 0.143 751, 0.015 806, 0.126 382 0).该平衡点处的雅可比矩阵为：

矩阵A|H1有 一对共轭纯虚根 λ3,4=±ωi， 其中 ω =0.105 7.我们取：

其中q为 特征值 λ3的 特征向量.它满足Aq=iωq，p为 伴随特征向量满足ATp=-iωp和 正交化条件 〈p,q〉=1.〈p,q〉=p1q1+p2q2+p3q3+p4q4是 C4中的标准数量积.

将平衡点转换到坐标原点，进行以下转换：

式中：x0、y0、z0、φ0为平衡点的坐标.通过这个变换，系统(2)被转换成：

考虑系统：

其中在x=0的邻域，A=A|H1,F(x)=B(x,y)和C(x,y,z)是多线性向量函数，x=(x1,x2,x3,x4)T，y=(y1,y2,y3,y4)T，z=(z1,z2,z3,z4)T，在坐标系中：

式中： ξ =(ξ1,ξ2,ξ3,ξ4).

由文献[27]中(5.62)可以求得H1的 第1 L yapunov系数为：

因此H1是一个超临界的Hopf 分岔，可以分支出稳定的极限环，即神经元通过超临界Hopf 分岔实现从放电状态到静息状态的转迁.

利用ISIs 分岔图探究神经元模型产生的更加详细的分岔现象来描述其丰富的动力学性质，如图2(a)所示.从图2(a)中可以看出，当 -2.6 ＜s＜-1.402 134 时 神经元模型表现出复杂的动力学行为.当-2.6 ＜s＜-2.2时 表现为逆加周期分岔，见图2(b).在s= -2.2 附 近逆加周期分岔消失，随着s的进一步增加神经元模型经由倍周期分岔进入混沌状态，最后在s=-1.402 134左右结束放电状态，即系统经由Hopf 分岔结束放电状态.

图2 参数s 的ISIs 分岔图Fig.2 The bifurcation diagram of ISIs with respect to s

由图3(a)可以清楚地看出s在 [-1.68,-1.62]范围内变化时系统从简单的单峰放电开始，经历了周期2、4 及周期8 等放电模式，随后系统经由倍周期分岔进入混沌状态.动力系统的 Lyapunov指数谱是用来判断动力系统是否混沌的指标之一.如果最大 L yapunov指 数为正而第2 个 L yapunov指数为0，则系统是混沌的.我们计算系统的第1 和第2 L yapunov指数来说明混沌现象，如图3(b)所示（红色和蓝色曲线分别表示系统的第1 和第2 L yapunov指数）.这与上面的描述是一致的.

图3 倍周期分岔Fig.3 The period-doubling bifurcation

其次研究感应电流反馈系数k0引 起的分岔现象.由图4(a)可以看出，当k0较小时系统处于混沌状态.随着k0的 增加系统经历周期3、周期6 等放电模式，随后经由倍周期分岔进入混沌状态.当k0增加到约0.054 5 时，系统通过逆倍周期分岔回到周期3 放电.随着k0进一步增大到0.088 时，系统进入混沌状态，随后通过2 次逆倍周期分岔最后回到规律的峰放电.我们计算第1 和第2 L yapunov系数来说明混沌行为，如图4(b)所示（红色和蓝色曲线分别表示系统的第1 和第2 L yapunov指数），与我们前面描述的行为一致.显然，感应电流对神经元膜电位有着明显的影响.

图4 参数 k0的ISIs 分岔图Fig.4 The bifurcation diagram of ISIs with respect to k0

3 放电模式分析

钙离子浓度变化速率比其他变量小几个数量级，即系统(1)中变量z的变化速率比其他3 个变量慢.因此将变量z作 为慢变量，x、y、φ看作快变量.由此得到快变子系统：

和慢变子系统：

我们将模型进行快慢分解，进而探究不同簇放电模式的产生机制.通过对快子系统的数值计算，发现在不同的参数s，感应电流反馈系数k0下 神经元可以产生4 种不同的簇放电模式.将慢变量z作为快速子系统(6)的1 个分岔参数，用Matcont 软件给出快子系统(6)分岔分析的结果.

3.1 “fold/fold”型簇放电当s=-1.55,k0=0.05时，系统(1)产生如图5(a)所示的簇放电模式.运用快慢动力学分岔分析和Matcont 软件进行数值计算得到快变子系统(6)的平衡点关于慢变量z的分岔曲线为Z 型曲线，并将簇放电轨线叠加至图中，如图5(b)所示.该Z 型分岔曲线上有2 个鞍结分岔点L1和L2，蓝色的实线和虚线分别表示稳定和不稳定的的平衡点.

图5 当 s=-1.55,k0=0.05时，产生“fold/fold”型簇放电Fig.5 Fast-slow dynamics of "fold/fold" whens=-1.55,k0=0.05

从图5(b)可以看出，下静息状态通过鞍结分岔点L2转迁至Z 型分岔曲线上支对应的放电状态，随后经由鞍结分岔点L1跳跃进入分岔曲线下支的下静息状态，簇放电上下状态转迁对应的分岔点都是鞍结分岔.由文献[25]对放电模式的分类，这种簇放电模式应归类为“fold/fold”型簇放电.

3.2 “fold/Hopf”型簇放电当s=-1.55,k0=0.02时，系统(1)产生如图6(a)所示的第2 种簇放电模式.图6(b)为快慢动力学分岔分析图，该Z 型分岔曲线的上支有2 个超临界Hopf 分岔点H1和H2，可分支出稳定的极限环.L1和L2为鞍结分岔点，图中蓝色实线和虚线分别表示稳定和不稳定的平衡点、红色实线表示稳定极限环的最大值和最小值.

图6 当 s=-1.55,k0=0.02时，产生“fold/Hopf”型簇放电Fig.6 Fast-slow dynamics of "fold/Hopf" whens=-1.55,k0=0.02

从图6(b)可以看出，随着慢变量z的减小，下静息状态经由鞍结分岔L2转迁至分岔曲线上支稳定极限环对应的放电状态.轨线沿着稳定的极限环振荡后经由Hopf 分岔点H2进入较低的静息状态.由簇放电类型的分类可知该种簇放电模式为“fold/hopf”型簇放电.

3.3 “fold/homoclinic”型簇放电除上述2 种类型的簇放电模式外，通过改变参数s，系统(1)还能呈现出如图7(a)所示的“fold/homoclinic”型簇放电.图7(b)为快慢动力学分岔分析图，该Z 型分岔曲线上支有1个超临界Hopf 分岔点H1，H2为 中性鞍点，Hopf 分岔点H1分支出来的稳定极限环与分岔曲线的中支相交形成鞍同宿轨道，L1和L2为 鞍结分岔点，HC为鞍结同宿分岔，图中蓝色的实线和虚线分别表示稳定和不稳定的平衡点、红色实线表示稳定极限环的最大值和最小值.

如图7(b)所示，随着慢变量z的减小，下静息状态经由鞍结分岔点L1转迁至分岔曲线上支稳定极限环所对应的放电状态，随后经由鞍同宿轨道转迁至分岔曲线下支对应的静息状态.这类簇放电为“fold/homoclinic”型簇放电，也称为方波型簇放电.

3.4 “fold/subHopf”型簇放电通过改变参数s，系统(1)还能产生如图8(a)所示的“fold/subHopf”型簇放电.图8(b)为快慢动力学分岔分析图，该Z 型分岔曲线上支有1 个亚临界Hopf 分岔点Hsub，可以分支出不稳定的极限环.L1和L2为鞍结分岔点，图中蓝色的实线和虚线分别表示稳定和不稳定的平衡点、红色虚线表示不稳定极限环的最大值和最小值.

图8 当 s=-2.4,k0=0.02时，产生“fold/subHopf”型簇放电Fig.8 Fast-slow dynamics of "fold/subHopf" whens=-2.4,k0=0.02

由图8(b)可以看出随着慢变量z的减小，下静息状态经由鞍结分岔点L1转迁至分岔曲线上支的放电状态，经由亚临界Hopf 分岔点Hsub转迁至分岔曲线下支的静息状态.因此，这种簇放电模式可归类为“fold/subHopf”型簇放电.

由上述的分析过程可以看出，引入了磁通变量和感应电流的神经元模型随着参数s和k0的变化表现出丰富的放电模式.

4 双参数分岔分析

本节利用快变子系统在 (z,k0) 和 (z,s)相平面的双参数分岔图解释上一节4 种不同的簇放电模式之间的转迁机制，同时利用动力系统的分岔理论和数值模拟研究神经元模型在 (k0,s)相平面上的余维2 分岔.

4.1 快速子系统双参数分岔我们使用Matcont 软件计算余维 2分岔点，如图9 所示，其中蓝色曲线h1和h2为Hopf 分 岔 曲 线，黑 色 曲 线f1、f2、f3和f4为 鞍 结 分 岔 曲 线.图9 中 的Gi(i=1,2)表 示 广 义Hopf 分 岔，Gi(i=1,2)表 示Cusp 分岔，Bi(i=1,2)表示Bogdanov-Takens 分岔，可参考文献[27]来理解这些分岔点的含义.表1 给出了每个分岔点所对应的参数值，特征值和规范性参数.

表1 分岔点的相关数据Tab.1 Data related to bifurcation points

图9 快变子系统的余维2 分岔图Fig.9 Codimension-two bifurcation of the fast subsystem

现分析随着参数k0的 变化不同簇放电模式之间的转迁机制.对于固定参数k0，快变子系统分岔曲线上的余维1 分岔点对应图9(a)中相应k0值 水平线与余维2 分岔曲线的交点.当k0＞0.033 33时，快变子系统的分岔曲线上没有Hopf 分岔点，因此系统呈现出的簇放电模式均为“fold/fold”型簇放电，类似于图5.k0=0.032 19时H opf 分岔曲线h1与 鞍结分岔曲线f2相 交于p1.当k0＜0.032 19时快变子系统分岔曲线上分岔点的位置发生变化，Hopf 分岔点位于鞍结分岔点的左侧，当轨线经过鞍结分岔点会转迁到稳定极限环所对应的放电状态，此时簇放电类型为“fold/ Hopf”型，类似于图6.考虑到实际情况，只分析k0＞0的部分.

接着我们分析参数s变化引起的不同簇放电模式的转迁动力学行为.对于固定的参数s，快变子系统分岔曲线上的余维1 分岔点对应图9(b)中相应s值水平线与余维2 分岔曲线的交点.当s＞-1.53时，快变子系统的分岔曲线上不存在 Hopf 分 岔点，因此系统表现出的簇放电类型均为“fold/ fold”型簇放电.随着s的减小，Hopf 分岔曲线h2与 鞍结分岔曲线f3相 交于点p3(-0.000 3,-1.530 4).当s＜-1.530 4时，分岔曲线上支分岔点的位置发生变化，Hopf 分岔点位于鞍结分岔点的左侧，当轨线经由鞍结分岔点时会转迁至稳定极限环对应的放电状态，此时簇放电模式为“fold/Hopf”型，类似于图6.

Hopf 分岔曲线h2和 鞍结分岔曲线f3相 切与B2点后，Hopf 分岔曲线的下支消失，图9(b)中Hopf 分岔曲线上的点变为中性鞍点.这表明快变子系统Z 型分岔曲线上的Hopf 分岔点与鞍结分岔点碰撞后消失，对应的极限环形成鞍点同宿轨道.从“fold/Hopf”型簇放电向“fold/homoclinic”型簇放电的转变正是因为穿过了同宿分支曲线[12].

广义 Hopf分岔的存在意味着实现了超临界Hopf 分岔与亚临界Hopf 分岔之间的转变.对于快变子系统Z 型分岔曲线来说，曲线上支的超临界Hopf 分岔变为亚临界 Hopf 分岔.当 -2.039 5 ＜s＜-1.673 5时，系统处于峰放电模式，如图2(a)所示.Hopf 分岔曲线h2与 鞍结分岔曲线f3相 交于p4(-0.000 02,-2.039 5)，随后快子系统Z 型分岔曲线上Hopf 分岔点与鞍结分岔点的位置发生变化.当s＜-2.039 5时，分岔曲线上支的亚临界Hopf 分岔点位于2 个鞍结分岔点之间，系统呈现“fold/subHopf”型簇放电，类似于图8.

本节利用快变子系统的双参数分岔理论对4 种簇放电模式间的转迁过程进行了详细分析，发现簇放电模式间转迁的关键在于快变子系统Z 型分岔曲线上的Hopf 分岔点、极限环和鞍结分岔点的位置关系.

4.2 全系统双参数分岔本节我们利用分岔理论与数值计算研究神经元模型关于参数k0和s的余维2 分岔，其他参数保持不变.利用Matcont 软件包计算得到全系统在 (k0,s)平面上的双参数分岔图，如图10 所示.其中h1和h2表 示Hopf 分岔曲线、Gi(i=3,4,5)代表广义Hopf分岔、HZ表示Zero-Hopf 分岔.表1 给出了每个分岔点所对应的参数值，特征值和规范性参数.

图10 全系统的余维2 分岔图Fig.10 Codimension-two bifurcation of the system

图10 中有3 个广义的Hopf 分岔点Gi(i=3,4,5)，1 个Zero-Hopf 分岔点HZ，并且大多数余维2 分岔点都在曲线h1上 ，2 条Hopf 分岔曲线h1和h2没 有交点.我们主要关注具有生理意义的参数区域 (k0＞0,s＜0).通过双参数分岔图确定神经元的放电区域大致位于Hopf 分岔曲线h1所围成的区域.

对于广义Hopf 分岔点来说，它的第1 个Lyapunov 系数为0.在Gi(i=3,4,5)岔点附近，系统(1)局部拓扑等价于规范形：

β1，β2∈R， 当实特征值为正时

对于Zero-Hopf 分岔点HZ来说，此时系统的雅克比矩阵有1 个0 特征值、1 对纯虚数特征值、1 个非零实部的特征值，但在分岔点附近没有固定的规范形.

5 结论

本文在三维的类Hindmarsh-Rose 模型中引入磁通变量和感应电流后构建了一个新四维神经元模型，利用分岔理论和快慢动力学分析研究该四维模型产生的分岔现象和放电模式.利用余维1 分岔研究了与神经元静息与放电状态转迁有关的分岔，着重计算Hopf 分岔点的第1L yapunov系数以确定其稳定性.通过数值计算系统的ISIs 分岔，发现ISIs 序列中存在逆加周期、逆倍周期及倍周期分岔导致的混沌现象，通过计算系统的 Lyapunov指数谱确定产生混沌的区域范围.利用快慢动力学分析得到了4 种簇放电模式及其产生机理，分别是当s=-1.55,k0=0.05时 系统表现为“fold/fold”型簇放电；s=-1.55,k0=0.02时系统表现为“fold/Hopf”型簇放电；s=-1.65,k0=0.02 时系统表现为“fold/homoclinic”型簇放电和s=-2.4,k0=0.02时的“fold/subHopf”型簇放电.接着利用快变子系统的双参数分岔理论对4 种簇放电模式间的转迁进行了详细分析，发现簇放电模式间转迁的关键在于快变子系统Z 型分岔曲线上Hopf 分岔点、极限环和鞍结分岔点的位置关系.最后利用全系统的余维2 分岔得到神经元处于放电状态的参数区域，求出全系统在参数平面(k0,s)上 的余维2 分岔点并给出分岔点附近的拓扑规范形.以上结果表明，参数s和感应电流反馈系数k0可以对神经元膜电位产生非常明显的影响.