Hopf Bifurcation Analysis of Magnetic Flux Ghostburster Neuron Model with Time Delay

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

其中：V_s，V_d分别为胞体、树突跨膜电压，其单位均为mV；I_s为流入体细胞的外部刺激电流，其单位为mA·cm^-2；n_s，n_d分别为胞体、树突细胞滞后K⁺电流的激活变量；h_d为树突细胞Na⁺电流失活变量；p_d为树突细胞滞后K⁺电流失活变量；细胞胞体与树突两个部分之间通过简单的耗散电流耦合，(I_Na，s，I_Dr，s，I_Na，d，I_Dr，d)代表离子电流，均受最大电导值g_max即g_Na，s，g_Dr，s，g_Na，d，g_Dr，d的影响，g_L表示泄露电导率，其单位均为mS·cm^-2；σ是一个时间通道常数即σ_NS，σ_HD，σ_ND，σ_PD，其单位为ms；其中，

是关于系统稳态函数的激活或者失活变量，且其中V_y，k_y为该稳态函数的中点电压和斜率系数；E_Na，E_K，E_L分别为Na⁺，K⁺和泄漏电流的反转电位，其单位均为mV；φ为通过神经元细胞膜的磁通；k₁，k₂，k₃为反馈增益；W(φ)为由磁通量控制忆阻器的记忆电导，其表达式为W(φ)= $\frac{{{\rm{d}}q(\varphi )}}{{{\rm{d}}\varphi }}$=α+3βφ²(其中α，β为确定的参数)；τ为时滞，其单位为ms.

设系统(1)的非零平衡点为χ^*=(V_s^*，n_s^*，V_d^*，h_d^*，n_d^*，p_d^*，φ^*)，利用变量代换$\tilde V$_s(t)=V_s(t)-V_s^*，$\tilde n$_s(t)=n_s(t)-n_s^*，$\tilde V$_d(t)=V_d(t)-V_d^*，$\tilde h$_d(t)=h_d(t)-h_d^*，$\tilde n$_d(t)=n_d(t)-n_d^*，$\tilde p$_d(t)=p_d(t)-p_d^*，$\tilde {\varphi} $(t)=φ(t)-φ^*，则系统(1)可化为

其中：

此时，系统(1)的平衡点χ^*=(V_s^*，n_s^*，V_d^*，h_d^*，n_d^*，p_d^*，φ^*)转化成系统(2)的平衡点(0，0，0，0，0，0，0). 则系统(2)在平衡点(0，0，0，0，0，0，0)处所对应的特征方程为

其中：A_i(i=0，1，…，5)，B_j(0，1，…，4)，C_k(k=1，2，…，7)均为特征方程(3)相应的系数.

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

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

其中：M₁₆=C₁，M₁₅=A₀+C₂，M₁₄=B₀+A₁+C₃，M₁₃=B₁+A₂+C₄，M₁₂=B₂+A₃+C₅，M₁₁=B₃+A₄+C₆，M₁₀=B₄+A₅+C₇.

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

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

(i) M₁₆>0；

(ii) M₁₅M₁₆>M₁₄；

(iii) M₁₄M₁₅M₁₆+M₁₂M₁₆>M₁₃M₁₆²+M₁₄²；

(iv) M₁₆ M₁₅(1-M₁₂M₁₅-M₁₀)+M₁₃(2M₁₂-M₁₃M₁₀) +M₁₄(M₁₂M₁₅+M₁₀-M₁₁M₁₆)-M₁₄M₁₃M₁₄-M₁₂²>0；

(v) $\prod\limits_{i = 2}^6 {{M_{1i}}} $-M₁₆²(M₁₃(M₁₀M₁₅-M₁₁M₁₄+M₁₂M₁₃)+M₁₁(M₁₁M₁₆-M₁₂M₁₅-2M₁₀))+ M₁₅M₁₆(M₁₀M₁₄M₁₅-M₁₁M₁₄²-M₁₂²M₁₅-M₁₀M₁₂)-M₁₄²(M₁₀M₁₅+M₁₂M₁₃)+ M₁₂²(2M₁₃M₁₆+M₁₄M₁₅-M₁₂)+M₁₄M₁₂(-3M₁₁M₁₆+2M₁₀)+M₁₁M₁₄³-M₁₆M₁₀²>0；

(vi) $\prod\limits_{i = 1}^6 {{M_{1i}}} $-M₁₀M₁₅M₁₆(M₁₃(3M₁₁M₁₆-M₁₂M₁₅+M₁₃M₁₄)+M₁₅(M₁₀-2M₁₁M₁₄))+M₁₁²M₁₆² (-M₁₁M₁₆+2M₁₂M₁₅+M₁₃M₁₄)+M₁₁M₁₆(-M₁₁M₁₄²M₁₅-M₁₂²M₁₅²-M₁₂M₁₃²M₁₆)+ M₁₀²M₁₃²M₁₆²+M₁₀²M₁₅(3M₁₃M₁₆+M₁₄M₁₅-M₁₂)+M₁₀M₁₄²(-2M₁₁M₁₅+M₁₃²)+ M₁₀M₁₂M₁₃(-2M₁₃M₁₆-M₁₄M₁₅)-M₁₁²M₁₄(3M₁₂M₁₆-M₁₄²)+M₁₁M₁₂²(2M₁₃M₁₆+M₁₄M₁₅)- M₁₁M₁₂M₁₄·(M₁₃M₁₄-3M₁₀)-M₁₀²(3M₁₁M₁₆+2M₁₃M₁₄-M₁₀)+M₁₂²(M₁₀M₁₃-M₁₁M₁₂)- M₁₀M₁₁M₁₆(M₁₂M₁₅+M₁₃M₁₄-3M₁₁M₁₆)>0；

(vii) M₁₀>0.

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

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

其中：

由式(6)可得：

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

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

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

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

其中：N₁₁，N₁₂，N₂₁，N₂₂的表达式可由λ=iω₀代入式(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，τ₀)时，系统(1)的平衡点χ^*是局部渐近稳定的；当τ≥τ₀时，系统(1)在平衡点χ^*处是不稳定的；当τ=τ₀时，系统(1)在平衡点χ^*处发生Hopf分岔.

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

设τ=τ₀+μ，x₁(t)=V_s(t)-V_s^*，x₂(t)=n_s(t)-n_s^*，x₃(t)=V_d(t)-V_d^*，x₄(t)=h_d(t)-h_d^*，x₅(t)=n_d(t)-n_d^*，x₆(t)=p_d(t)-p_d^*，x₇(t)=φ(t)-φ^*，x(t)=(x₁(t)，x₂(t)，x₃(t)，x₄(t)，x₅(t)，x₆(t)，x₇(t))^T∈ $\mathbb{R}$⁷. 将系统(1)在定义域C=C([-1, 0]，$\mathbb{R}$⁷)上转化为泛函微分方程：

定义(12)式中的线性算子L_μ：C $\mathbb{R}$⁷如下：

其中：ϕ(t)=(ϕ₁(t)，ϕ₂(t)，ϕ₃(t)，ϕ₄(t)，ϕ₅(t)，ϕ₆(t)，ϕ₇(t))^T∈C，

定义(12)式中非线性算子F：C$\mathbb{R}$⁷如下：

其中：F₁=g_Na，sm_∞，s²ϕ₁₀(E_Na-ϕ₁₀)(1-ϕ₂₀)+g_Dr，sϕ₂₀²(E_K-ϕ₁₀)-3k₁β(ϕ₁₀ϕ₂₀²-ϕ₃₁ϕ₇₀²)，F₂=n_∞，s²ϕ₁₀σ_NS^-1，F₃=g_Na，dm_∞，d²ϕ₃₀ϕ₄₀(E_Na-ϕ₃₀)+g_Dr，dϕ₆₀ϕ₅₀²(E_K-ϕ₃₀)-3k₁βϕ₇₀²(ϕ₃₀-ϕ₁₁)，F₄=h_∞，dϕ₃₀σ_HD^-1，F₅=n_∞，dϕ₃₀σ_ND^-1，F₆=p_∞，dϕ₃₀σ_PD^-1，ϕ_i0=ϕ_i(0)，ϕ_j1=ϕ_j(-1)，i，j=1，2，…，7.

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

其中ϕ(θ)∈C，η(θ，μ)=(τ₀+μ)η₀ϕ(θ)-(τ₀+μ)η_-1ϕ(θ+1). 定义：

系统(12)等同于

其中x_t(θ)=x(t+θ)，θ∈[-1，0]. 对于ψ(l)∈C([0, 1]，$\mathbb{R}$⁷)，定义A^*(0)和双线性内积如下所示：

其中：η(θ)=η(θ，0)，ψ(0)是ψ(0)的共轭；A(0)和A^*(0)是伴随算子. 显然，±iω₀τ₀也是A^*(0)的特征值. 设q(θ)，q^*(l)分别是A(0)，A^*(0)对应特征值iω₀τ₀，-iω₀τ₀的特征向量，有A(0)q(θ)=iω₀τ₀q(θ)，A^*(0)q^*(l)=-iω₀τ₀q^*(l).

令q(θ)=(1，v₁，v₂，v₃，v₄，v₅，v₆)^Te^{iω₀τ₀θ}，q^*(l)=G(1，v₁^*，v₂^*，v₃^*，v₄^*，v₅^*，v₆^*)^Te^iω₀τ₀l，利用

可得出v_i，v_i^*(i=1，2，…，6)的具体表达式.

由式(19)可得

其中，q^*(0)为q^*(0)的共轭，v_i^*(i=1，2，…，6)是v_i^*的共轭.

为了满足〈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时中心流形C₀的坐标：

其中：

从而，计算可得

由以上讨论可以得出：

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

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

3) 系统分岔周期解的周期由T₂决定：当T₂>0(T₂ < 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峰值，其出现的峰值减少，整体尖峰放电态出现了延迟现象.

当参数g_Na，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)在I_s和g_Na，d双参数下的分岔结构，并通过改变时滞来分析时滞对系统放电行为产生的影响. 其他参数与时间序列图中参数相同. 图 3中不同的颜色表示系统(1)不同的放电状态，右边颜色栏的数值表示相应的周期簇放电态，如数字2表示周期2簇放电态，白色区域表示周期大于20簇放电态或者混沌放电态.

由图 3(a)可知，当g_Na，d∈[7，12.5](单位mS·cm^-2)，I_s∈[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)中当g_Na，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，τ₀)时，系统(1)的平衡点是局部渐近稳定的；当τ≥τ₀时，平衡点是不稳定的；当τ=τ₀时，系统(1)在平衡点处存在Hopf分岔. 利用范式理论和中心流形定理进一步证明了μ₂决定Hopf分岔的方向，β₂控制分岔周期解的稳定性，T₂决定分岔周期解的周期. 最后，利用Matlab数值模拟得出系统(1)在不同时滞下时间序列图、钠离子电流最大电导值的峰峰间期(ISI)分岔图以及关于外部刺激电流值和钠离子电流最大电导值的双参数分岔图. 在时间序列图中随着时滞的不断增大，发现系统(1)尖峰放电行为逐步延迟. 在峰峰间期(ISI)分岔图与双参数分岔图中发现系统(1)的放电行为呈现出加周期分岔模式，并随周期数的不断增大，周期窗口逐渐变小，且当时滞不断增大时，放电模式发生了明显的改变. 研究结果验证了时滞是影响Ghostburster神经元系统放电模式的一个重要因素，也很好地解释了电磁辐射作用下具有延迟效应的神经元系统的异常放电行为.