下载:
-
神经元是神经系统结构和功能的基本单位. 神经系统信息的产生、整合和传输是由神经元产生动作电位引起放电行为来实现. 在不同的外界电流刺激下，神经元呈现出丰富多样的放电模式，如静息态、峰放电、簇放电和混沌放电等；而由于系统参数的改变所造成不同放电模式之间的转换将会产生分岔行为，如加周期分岔、倍周期分岔、逆倍周期分岔等行为[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)分岔图以及双参分岔图分别分析该系统在不同时滞影响下的动力学行为. 通过对该模型在不同时滞下放电模式以及放电模式之间相互转换所引起的加周期分岔行为的动力学分析，不但能够揭示多种有趣的动力学现象和各组织细胞的运行规律，而且还可以更深入地认识脑神经系统的信息活动机理和感觉、学习、记忆和思维等认知功能的原理和复杂性，并对神经科学实验和医学诊断提供重要的参考依据.
全文HTML
-
Ghostburster模型[9-11]包括胞体与树突两部分，有神经元胞体与树突膜电流电压等6个变量，可以准确解释弱电鱼(ELL)锥体神经元细胞中细胞体与树突膜电位信息的转变过程. 基于原Ghostburster神经元模型，引入磁通量实现电磁感应电流对膜电位的调制[13]. 假设由胞体与树突之间的耦合由简单电子分布产生，从胞体流向树突或从树突流向胞体的介质传递过程中存在时间延迟，考虑在原Ghostburster模型的基础上加入时滞，得到神经元模型如下：
其中：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；其中，
是关于系统稳态函数的激活或者失活变量，且其中Vy，ky为该稳态函数的中点电压和斜率系数；ENa，EK，EL分别为Na+，K+和泄漏电流的反转电位，其单位均为mV；φ为通过神经元细胞膜的磁通；k1，k2，k3为反馈增益；W(φ)为由磁通量控制忆阻器的记忆电导，其表达式为W(φ)=
$\frac{{{\rm{d}}q(\varphi )}}{{{\rm{d}}\varphi }}$ =α+3βφ2(其中α，β为确定的参数)；τ为时滞，其单位为ms.
-
设系统(1)的非零平衡点为χ*=(Vs*，ns*，Vd*，hd*，nd*，pd*，φ*)，利用变量代换
$\tilde V$ s(t)=Vs(t)-Vs*，$\tilde n$ s(t)=ns(t)-ns*，$\tilde V$ d(t)=Vd(t)-Vd*，$\tilde h$ d(t)=hd(t)-hd*，$\tilde n$ d(t)=nd(t)-nd*，$\tilde p$ d(t)=pd(t)-pd*，$\tilde {\varphi} $ (t)=φ(t)-φ*，则系统(1)可化为其中：
此时，系统(1)的平衡点χ*=(Vs*，ns*，Vd*，hd*，nd*，pd*，φ*)转化成系统(2)的平衡点(0，0，0，0，0，0，0). 则系统(2)在平衡点(0，0，0，0，0，0，0)处所对应的特征方程为
其中：Ai(i=0，1，…，5)，Bj(0，1，…，4)，Ck(k=1，2，…，7)均为特征方程(3)相应的系数.
根据时滞值的不同情形，分以下两种情形讨论.
情形1 当τ=0时，原特征方程可以化简表示为
其中：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；
(iii) M14M15M16+M12M16>M13M162+M142；
(iv) M16 M15(1-M12M15-M10)+M13(2M12-M13M10) +M14(M12M15+M10-M11M16)-M14M13M14-M122>0；
(v)
$\prod\limits_{i = 2}^6 {{M_{1i}}} $ -M162(M13(M10M15-M11M14+M12M13)+M11(M11M16-M12M15-2M10))+ M15M16(M10M14M15-M11M142-M122M15-M10M12)-M142(M10M15+M12M13)+ M122(2M13M16+M14M15-M12)+M14M12(-3M11M16+2M10)+M11M143-M16M102>0；(vi)
$\prod\limits_{i = 1}^6 {{M_{1i}}} $ -M10M15M16(M13(3M11M16-M12M15+M13M14)+M15(M10-2M11M14))+M112M162 (-M11M16+2M12M15+M13M14)+M11M16(-M11M142M15-M122M152-M12M132M16)+ M102M132M162+M102M15(3M13M16+M14M15-M12)+M10M142(-2M11M15+M132)+ M10M12M13(-2M13M16-M14M15)-M112M14(3M12M16-M142)+M11M122(2M13M16+M14M15)- M11M12M14·(M13M14-3M10)-M102(3M11M16+2M13M14-M10)+M122(M10M13-M11M12)- M10M11M16(M12M15+M13M14-3M11M16)>0；(vii) M10>0.
情形2 当τ>0时，特征方程(3)左右两边乘以eλτ得
将λ=iω(ω>0)代入方程(5)得
其中：
由式(6)可得：
因此，可得关于ω的方程：
设方程(7)存在一个正实数根为ω0，则可以得
由于方程(5)是隐函数，在特征方程(5)两边对τ求导得
将λ=iω0代入式(9)有
其中：N11，N12，N21，N22的表达式可由λ=iω0代入式(9)整理得出.
那么，满足
${\mathop{\rm Re}\nolimits} \left[ {{{\left. {{{\left( {\frac{{{\rm{d}}\lambda }}{{{\rm{d}}\tau }}} \right)}^{ - 1}}} \right|}_{\lambda = {\rm{i}}{\omega _0}}}} \right] \ne 0$ 的充分条件为定理2 当系统(1)满足以下条件：
(i) 满足定理1中的条件；
(ii) 对于方程(7)存在正实根；
(iii) 满足条件(11).
则当时滞τ∈[0，τ0)时，系统(1)的平衡点χ*是局部渐近稳定的；当τ≥τ0时，系统(1)在平衡点χ*处是不稳定的；当τ=τ0时，系统(1)在平衡点χ*处发生Hopf分岔.
-
对于特征方程(7)的一个特征根ω0则存在一个τ0. 在每一个τ0处系统(1)发生Hopf分岔.下面利用中心流形定理与范式理论来判定分岔周期解的稳定性及分岔方向.
设τ=τ0+μ，x1(t)=Vs(t)-Vs*，x2(t)=ns(t)-ns*，x3(t)=Vd(t)-Vd*，x4(t)=hd(t)-hd*，x5(t)=nd(t)-nd*，x6(t)=pd(t)-pd*，x7(t)=φ(t)-φ*，x(t)=(x1(t)，x2(t)，x3(t)，x4(t)，x5(t)，x6(t)，x7(t))T∈
$\mathbb{R}$ 7. 将系统(1)在定义域C=C([-1, 0]，$\mathbb{R}$ 7)上转化为泛函微分方程：定义(12)式中的线性算子Lμ：C
$\mathbb{R}$ 7如下：其中：ϕ(t)=(ϕ1(t)，ϕ2(t)，ϕ3(t)，ϕ4(t)，ϕ5(t)，ϕ6(t)，ϕ7(t))T∈C，
定义(12)式中非线性算子F：C
$\mathbb{R}$ 7如下：其中：F1=gNa，sm∞，s2ϕ10(ENa-ϕ10)(1-ϕ20)+gDr，sϕ202(EK-ϕ10)-3k1β(ϕ10ϕ202-ϕ31ϕ702)，F2=n∞，s2ϕ10σNS-1，F3=gNa，dm∞，d2ϕ30ϕ40(ENa-ϕ30)+gDr，dϕ60ϕ502(EK-ϕ30)-3k1βϕ702(ϕ30-ϕ11)，F4=h∞，dϕ30σHD-1，F5=n∞，dϕ30σND-1，F6=p∞，dϕ30σPD-1，ϕi0=ϕi(0)，ϕj1=ϕj(-1)，i，j=1，2，…，7.
由Riesz表示定理[22]知，存在有界变差函数η(θ，μ)，(-1≤θ≤0)使得
其中ϕ(θ)∈C，η(θ，μ)=(τ0+μ)η0ϕ(θ)-(τ0+μ)η-1ϕ(θ+1). 定义：
系统(12)等同于
其中xt(θ)=x(t+θ)，θ∈[-1，0]. 对于ψ(l)∈C([0, 1]，
$\mathbb{R}$ 7)，定义A*(0)和双线性内积如下所示：其中：η(θ)=η(θ，0)，ψ(0)是ψ(0)的共轭；A(0)和A*(0)是伴随算子. 显然，±iω0τ0也是A*(0)的特征值. 设q(θ)，q*(l)分别是A(0)，A*(0)对应特征值iω0τ0，-iω0τ0的特征向量，有A(0)q(θ)=iω0τ0q(θ)，A*(0)q*(l)=-iω0τ0q*(l).
令q(θ)=(1，v1，v2，v3，v4，v5，v6)Teiω0τ0θ，q*(l)=G(1，v1*，v2*，v3*，v4*，v5*，v6*)Teiω0τ0l，利用
可得出vi，vi*(i=1，2，…，6)的具体表达式.
由式(19)可得
其中，q*(0)为q*(0)的共轭，vi*(i=1，2，…，6)是vi*的共轭.
为了满足〈q*(l)，q(θ)〉=1，〈q*(l)，q(θ)〉=0. 选择G=
${\left( {\left( {1 + \sum\limits_{i = 1}^6 {{v_i}\bar v_i^*} } \right) + {\tau _0}{{\rm{e}}^{ - {\rm{i}}{\omega _0}{\tau _0}}}\left( {{J_{31}}v_2^* + {J_{13}}} \right)} \right)^{ - 1}}$ .接下来，计算μ=0时中心流形C0的坐标：
其中：
从而，计算可得
由以上讨论可以得出：
1) μ2的符号决定Hopf分岔的方向：当μ2>0(μ2 < 0)，则Hopf分岔是超临界分岔(亚临界分岔).
2) β2控制分岔周期解的稳定性：当β2 < 0(β2>0)，则可得分岔周期解是稳定(不稳定)的.
3) 系统分岔周期解的周期由T2决定：当T2>0(T2 < 0)成立，则可得分岔周期解的周期是递增(递减)的.
-
本节将通过数值模拟进一步研究该系统的动力学行为. 当系统(1)受到外部刺激电流时，神经元的放电模式会发生变化. 系统(1)中的放电活动与系统中参数的选择密切相关. 取系统相关参数[11]：
图 1给出了系统(1)在不同时滞下的时间序列图，该4幅时间序列图呈现了在不同时滞作用下尖峰态的放电模式. 图 1中的虚线框内14，15，16表示第14，15，16个峰值. 在相同的时间t下随着时滞τ的增大，比较(a)，(b)，(c)和(d)知第14，15，16峰值出现延迟，在(c)，(d)中同一时间点未达到第16峰值，其出现的峰值减少，整体尖峰放电态出现了延迟现象.
当参数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)比较可知：最后发生的尖峰放电态的持续时间逐渐增加，这表明加入时滞后系统原来的放电模式发生改变.
现实中，系统参数可能相互关联. 图 3研究了系统(1)在Is和gNa，d双参数下的分岔结构，并通过改变时滞来分析时滞对系统放电行为产生的影响. 其他参数与时间序列图中参数相同. 图 3中不同的颜色表示系统(1)不同的放电状态，右边颜色栏的数值表示相应的周期簇放电态，如数字2表示周期2簇放电态，白色区域表示周期大于20簇放电态或者混沌放电态.
由图 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种不同，在原有的两两周期簇放电态之间穿插了尖峰放电态，时滞的增大导致放电模式紊乱.
-
本文在原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神经元系统放电模式的一个重要因素，也很好地解释了电磁辐射作用下具有延迟效应的神经元系统的异常放电行为.