含时滞的磁通Ghostburster神经元模型的Hopf分岔分析①

南梦冉， 张建刚， 魏立祥， 张美娇

兰州交通大学 数理学院， 兰州 730070

神经元是神经系统结构和功能的基本单位． 神经系统信息的产生、 整合和传输是由神经元产生动作电位引起放电行为来实现． 在不同的外界电流刺激下， 神经元呈现出丰富多样的放电模式， 如静息态、 峰放电、 簇放电和混沌放电等； 而由于系统参数的改变所造成不同放电模式之间的转换将会产生分岔行为， 如加周期分岔、 倍周期分岔、 逆倍周期分岔等行为[1-3]． 因此， 研究神经元放电活动的动力学特性具有一定的生理学意义． 目前针对不同的神经元细胞已构建了多种神经元模型， 如Hodgkin-Huxley(HH)模型[4]、 Morris-Lecar(ML)模型[5]、 Chay模型[6]和Hindmarsh-Rose(HR)模型[7]等． 1994年， 文献[8]研究了弱电鱼电敏感侧线叶(electrosensory lateral line lobe， ELL)锥体细胞中的胞体和树突间离子通道信息的传导模式以及放电规律； 2002年， 文献[9-10]讨论了在外电场作用下ELL椎体细胞的非线性特性， 将多室结构降维至二室并命名为Ghostburster模型； 随后在不同的参数空间内有关Ghostburster神经元模型的分岔行为与放电模式的报道不断涌现[2，10-12]．

引入忆阻器和磁通等元素后神经元将呈现出更为丰富的放电特性． 文献[13]研究了带有磁通量的神经元模型的动力学特性， 通过改变初始状态观察到多种模式的放电行为； 文献[14]通过建立具有磁控忆阻器的四维集成电路系统， 发现忆阻耦合可以有效地控制混沌电路， 也可以成功地再现神经元活动的动力学行为； 文献[15]在ML神经元模型上加入磁控忆阻器， 采用数值模拟的方法来研究电磁感应作用下由化学突触和电突触耦合的神经元系统的同步行为． 此外， 由于信号传输速度的有限性和递质释放的滞后性， 使得神经元在信息表达与传递的过程中出现了延迟现象． 文献[16-17]研究了两类具有忆阻耦合的神经元系统， 考虑神经元间信号传输过程中的时滞效应， 通过平衡点的稳定性分析获得与时滞相关的稳定性条件； 文献[18]等研究了含时滞磁通神经元模型的稳定性及Hopf分岔； 文献[19]研究了时滞对静息状态下耦合后的HR神经元动力学行为的影响， 发现时滞可以诱发丰富的尖峰模式．

基于以上研究与讨论， 本文在原Ghostburster模型的基础上分别加入忆阻器及时滞项， 并进一步研究了时滞对磁通Ghostburster神经元模型的放电行为的影响． 首先， 利用Routh-Hurwitz和Hopf分岔分析方法[20]给出系统平衡点的稳定性以及Hopf分岔存在的条件； 其次， 利用中心流形定理和范式理论研究了Hopf分岔的方向和周期解的稳定性； 最后， 选取树突膜钠离子电流最大电导值和胞体膜注入电流值作为分岔参数， 利用时间序列图、 峰峰间期(interspike interval， ISI)分岔图以及双参分岔图分别分析该系统在不同时滞影响下的动力学行为． 通过对该模型在不同时滞下放电模式以及放电模式之间相互转换所引起的加周期分岔行为的动力学分析， 不但能够揭示多种有趣的动力学现象和各组织细胞的运行规律， 而且还可以更深入地认识脑神经系统的信息活动机理和感觉、 学习、 记忆和思维等认知功能的原理和复杂性， 并对神经科学实验和医学诊断提供重要的参考依据．

1 模型的建立

Ghostburster模型[9-11]包括胞体与树突两部分， 有神经元胞体与树突膜电流电压等6个变量， 可以准确解释弱电鱼(ELL)锥体神经元细胞中细胞体与树突膜电位信息的转变过程． 基于原Ghostburster神经元模型， 引入磁通量实现电磁感应电流对膜电位的调制[13]． 假设由胞体与树突之间的耦合由简单电子分布产生， 从胞体流向树突或从树突流向胞体的介质传递过程中存在时间延迟， 考虑在原Ghostburster模型的基础上加入时滞， 得到神经元模型如下：

(1)

其中：Vs，Vd分别为胞体、 树突跨膜电压， 其单位均为mV；Is为流入体细胞的外部刺激电流， 其单位为mA·cm-2；ns，nd分别为胞体、 树突细胞滞后K+电流的激活变量；hd为树突细胞Na+电流失活变量；pd为树突细胞滞后K+电流失活变量； 细胞胞体与树突两个部分之间通过简单的耗散电流耦合， (INa，s，IDr，s，INa，d，IDr，d)代表离子电流， 均受最大电导值gmax即gNa，s，gDr，s，gNa，d，gDr，d的影响，gL表示泄露电导率， 其单位均为mS·cm-2；σ是一个时间通道常数即σNS，σHD，σND，σPD， 其单位为ms； 其中，

2 平衡点的局部稳定性与Hopf分岔分析

(2)

其中：

(3)

其中：Ai(i=0，1，…，5)，Bj(0，1，…，4)，Ck(k=1，2，…，7)均为特征方程(3)相应的系数．

根据时滞值的不同情形， 分以下两种情形讨论．

情形1 当τ=0时， 原特征方程可以化简表示为

λ7+M16λ6+M15λ5+M14λ4+M13λ3+M12λ2+M11λ+M10=0

(4)

其中：M16=C1，M15=A0+C2，M14=B0+A1+C3，M13=B1+A2+C4，M12=B2+A3+C5，M11=B3+A4+C6，M10=B4+A5+C7．

由Routh-Hurwitz判据[21]可得，

定理1系统(1)在平衡点χ*处是局部渐近稳定的当且仅当该特征方程(4)的所有根具有负实部， 即特征方程(4)的系数满足如下条件：

(i)M16>0；

(ii)M15M16>M14；

(iv)M16[M15(1-M12M15-M10)+M13(2M12-M13M10)]+M14(M12M15+M10-M11M16)-

M10M11M16(M12M15+M13M14-3M11M16)>0；

(vii)M10>0．

情形2 当τ>0时， 特征方程(3)左右两边乘以eλτ得

(5)

将λ=iω(ω>0)代入方程(5)得

(6)

其中：

M11(ω)=-C1ω6+(C3+A1)ω4-(C5+A3)ω2+(C7+A5)，

M12(ω)=ω7+(A0-C2)ω5+(C4-A2)ω3+(A4-C6)ω，M13(ω)=-B0ω4+B2ω2-B4，

M21(ω)=-ω7+(A0+C2)ω5-(C4+A2)ω3+(A4+C6)ω，

M22(ω)=-C1ω6+(C3-A1)ω4+(A3-C5)ω2+(C7-A5)，M23(ω)=B1ω3-B3ω．

由式(6)可得：

因此， 可得关于ω的方程：

(7)

设方程(7)存在一个正实数根为ω0， 则可以得

(8)

由于方程(5)是隐函数， 在特征方程(5)两边对τ求导得

(9)

将λ=iω0代入式(9)有

(10)

其中：N11，N12，N21，N22的表达式可由λ=iω0代入式(9)整理得出．

N11N21+N12N22≠0

(11)

定理2当系统(1)满足以下条件：

(i) 满足定理1中的条件；

(ii) 对于方程(7)存在正实根；

(iii) 满足条件(11)．

则当时滞τ∈[0，τ0)时， 系统(1)的平衡点χ*是局部渐近稳定的； 当τ≥τ0时， 系统(1)在平衡点χ*处是不稳定的； 当τ=τ0时， 系统(1)在平衡点χ*处发生Hopf分岔．

3 Hopf分岔方向与周期解的稳定性

对于特征方程(7)的一个特征根ω0则存在一个τ0． 在每一个τ0处系统(1)发生Hopf分岔． 下面利用中心流形定理与范式理论来判定分岔周期解的稳定性及分岔方向．

(12)

Lμ(φ(t))=(τ0+μ)·(η0φ(0)+η-1φ(-1))

(13)

其中：φ(t)=(φ1(t)，φ2(t)，φ3(t)，φ4(t)，φ5(t)，φ6(t)，φ7(t))T∈C，

F(μ，φ(t))=(τ0+μ)·(F1，F2，F3，F4，F5，F6， 0)T

(14)

由Riesz表示定理[22]知， 存在有界变差函数η(θ，μ)， (-1≤θ≤0)使得

(15)

其中φ(θ)∈C，η(θ，μ)=(τ0+μ)η0φ(θ)-(τ0+μ)η-1φ(θ+1)． 定义：

(16)

系统(12)等同于

(17)

其中xt(θ)=x(t+θ)，θ∈[-1， 0]． 对于ψ(l)∈C([0， 1]， R7)， 定义A*(0)和双线性内积如下所示：

(18)

(19)

由式(19)可得

接下来， 计算μ=0时中心流形C0的坐标：

其中：

H3=(H31，H32，H33，H34，H35，H36， 0)T，H4=(H41，H42，H43，H44，H45，H46， 0)T，

H5=(H51，H52，H53，H54，H55，H56，H57)，H6=(H61，H62，H63，H64，H65，H66，H67)，

H51=(J11，J21，J31e-2iω0τ0， 0， 0， 0，J71)T，H52=(J12，J22， 0， 0， 0， 0， 0)T，

H53=(J13e-2iω0τ0， 0，J33，J43， 0， 0，J73)T，H54=(0， 0，J34，J44， 0， 0， 0)T，

H55=(0， 0，J35，J55， 0， 0， 0)T，H56=(0， 0，J36， 0， 0，J66， 0)T，

H57=(J17， 0，J37， 0， 0， 0，J77)T，H61=(J11，J21，J31， 0， 0， 0，J71)T，

H62=(J12，J22， 0， 0， 0， 0， 0)T，H63=(J13， 0，J33，J43，J53，J63，J73)T，

H64=(0， 0，J34，J44， 0， 0， 0)T，H65=(0， 0，J35， 0，J55， 0， 0)T，

H66=(0， 0，J35， 0， 0，J66， 0)T，H67=(J17， 0，J37， 0， 0， 0，J77)T，

从而， 计算可得

由以上讨论可以得出：

1)μ2的符号决定Hopf分岔的方向： 当μ2>0(μ2<0)， 则Hopf分岔是超临界分岔(亚临界分岔)．

2)β2控制分岔周期解的稳定性： 当β2<0(β2>0)， 则可得分岔周期解是稳定(不稳定)的．

3) 系统分岔周期解的周期由T2决定： 当T2>0(T2<0)成立， 则可得分岔周期解的周期是递增(递减)的．

4 数值模拟

本节将通过数值模拟进一步研究该系统的动力学行为． 当系统(1)受到外部刺激电流时， 神经元的放电模式会发生变化． 系统(1)中的放电活动与系统中参数的选择密切相关． 取系统相关参数[11]：

VMS=-40 mV；kMS=3 mV；VMD=-40 mV；kMD=5 mV；VNS=-40 mV；kNS=3 mV；

VHD=-52 mV；kHD=-5 mV；VND=-40 mV；Is=6.96 mA·cm-2；kND=5 mV；VPD=-65 mV；

kPD=-6 mV；gNa，s=32 mS·cm-2；gDr，s=20 mS·cm-2；gDr，d=5 mS·cm-2；gL=0.18 mS·cm-2；

ENa=52 mV；EK=-88.5 mV；gNa，d=5 mS·cm-2；EL=-37 mV；σNS=0.39 ms；σHD=1 ms；

σND=0.91 ms；σPD=5.8 ms；α=0.02；β=0.3；k1=0.2；k2=0.8；k3=0.5．

图1给出了系统(1)在不同时滞下的时间序列图， 该4幅时间序列图呈现了在不同时滞作用下尖峰态的放电模式． 图1中的虚线框内14，15，16表示第14，15，16个峰值． 在相同的时间t下随着时滞τ的增大， 比较(a)，(b)，(c)和(d)知第14，15，16峰值出现延迟， 在(c)，(d)中同一时间点未达到第16峰值， 其出现的峰值减少， 整体尖峰放电态出现了延迟现象．

图1 不同时滞下系统(1)的时间序列图

当参数gNa，d在[0，12.5]上变化时， 系统单参分岔行为如图2所示． 图2(a)中标有2，3，…，10的数字分别代表2，3，…，10周期簇放电态． 观察图2(a)可知系统(1)刚开始处于2周期簇放电态， 随着钾离子电流最大电导值的增大， 经过加周期分岔模式从2周期逐渐增加， 最终通向混沌放电态或者更高周期簇放电态． 并且随着周期数的增加， 对应的周期簇放电态所持续的时间变短， 即周期窗口缩小． 随着时滞的不断增大， 整体上图2中系统加周期分岔模式未发生改变而是系统(1)初始状态2周期簇放电态的区域逐步减少， 可以明显看到图2(d)中系统(1)原来初始时刻处的2周期簇放电态已消失． 通过图2(b)，(c)，(d)比较可知： 最后发生的尖峰放电态的持续时间逐渐增加， 这表明加入时滞后系统原来的放电模式发生改变．

图2 不同时滞下系统(1)关于参数gNa，d峰峰间期(ISI)分岔图

现实中， 系统参数可能相互关联． 图3研究了系统(1)在Is和gNa，d双参数下的分岔结构， 并通过改变时滞来分析时滞对系统放电行为产生的影响． 其他参数与时间序列图中参数相同． 图3中不同的颜色表示系统(1)不同的放电状态， 右边颜色栏的数值表示相应的周期簇放电态， 如数字2表示周期2簇放电态， 白色区域表示周期大于20簇放电态或者混沌放电态．

图3 系统(1)在gNa，d和Is双参数下不同时滞的分岔图

由图3(a)可知， 当gNa，d∈[7，12.5](单位mS·cm-2)，Is∈[8，9.7](单位mA·cm-2)时， 系统(1)通过加周期分岔模式从2周期一直增加到19周期， 再进入混沌放电态， 最后进入1周期尖峰放电态． 并且从3周期开始， 随着周期数的不断增加， 每个周期数所代表的颜色带逐渐变窄． 如4周期颜色带比3周期颜色带窄， 5周期颜色带比4周期颜色带更窄． 图3(a)，(b)，(c)，(d)都呈现出加周期分岔模式的放电状态， 然而随着时滞τ的不断增大， 图中2周期簇放电态的颜色带逐渐变窄直至最后消失， 这与图2峰峰间期(ISI)分岔图的放电模式吻合． 系统(1)的整个放电状态随时滞的增大而发生改变， 这说明时滞对神经元的放电模式产生了一定的影响． 在图3(a)中当gNa，d取值区间分别为[8.3，9.1]， [9.5， 10]， [10.3， 10.6](单位 mS·cm-2)时， 3，4和5周期簇放电态中通过倍周期分岔分别出现了6，8和10周期簇放电态， 周期6簇放电态通过逆倍周期分岔出现了周期3簇放电态， 且周期6簇放电态与周期3簇放电态相互交替出现． 但随着时滞的增大， 在图3(b)，(c)和(d)中上述现象逐渐消失． 由图3(d)可知，τ=0.2 ms时的放电模式与前3种不同， 在原有的两两周期簇放电态之间穿插了尖峰放电态， 时滞的增大导致放电模式紊乱．

5 结论

本文在原Ghostburster模型的基础上加入磁控忆阻器， 并在树突和胞体之间加入时滞τ， 建立了一个含时滞的磁通神经元模型， 以更好地探讨Ghostburster神经元模型的动力学行为． 通过Routh-Hurwitz判据讨论了τ=0 ms时系统(1)在平衡点处局部渐近稳定的条件． 利用Hopf分岔定理证明了当τ∈[0，τ0)时， 系统(1)的平衡点是局部渐近稳定的； 当τ≥τ0时， 平衡点是不稳定的； 当τ=τ0时， 系统(1)在平衡点处存在Hopf分岔． 利用范式理论和中心流形定理进一步证明了μ2决定Hopf分岔的方向，β2控制分岔周期解的稳定性，T2决定分岔周期解的周期． 最后， 利用Matlab数值模拟得出系统(1)在不同时滞下时间序列图、 钠离子电流最大电导值的峰峰间期(ISI)分岔图以及关于外部刺激电流值和钠离子电流最大电导值的双参数分岔图． 在时间序列图中随着时滞的不断增大， 发现系统(1)尖峰放电行为逐步延迟． 在峰峰间期(ISI)分岔图与双参数分岔图中发现系统(1)的放电行为呈现出加周期分岔模式， 并随周期数的不断增大， 周期窗口逐渐变小， 且当时滞不断增大时， 放电模式发生了明显的改变． 研究结果验证了时滞是影响Ghostburster神经元系统放电模式的一个重要因素， 也很好地解释了电磁辐射作用下具有延迟效应的神经元系统的异常放电行为．