含時滯的磁通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神經元系統放電模式的一個重要因素， 也很好地解釋了電磁輻射作用下具有延遲效應的神經元系統的異常放電行為．