基于麥克斯韋電磁場理論的神經元動力學響應與隱藏放電控制*

安新磊 喬帥 張莉

1) (蘭州交通大學數理學院, 蘭州 730070)

2) (蘭州理工大學電氣工程與信息工程學院, 蘭州 730050)

3) (蘭州工業學院基礎學科部, 蘭州 730050)

鈣、鉀、鈉等離子在細胞內連續泵送和傳輸時產生的時變電場不僅會影響神經元的放電活動, 而且會誘導時變磁場去進一步調節細胞內離子的傳播.根據麥克斯韋電磁場理論, 時變的電場和磁場在細胞內外的電生理環境中會相互激發而產生電磁場.為了探究電磁場影響下的神經元放電節律轉遷, 本文在三維Hindmarsh-Rose (HR)神經元模型的基礎上, 引入磁通變量和電場變量, 建立了一個五維HR 神經元模型(簡稱EMFN 模型).首先, 結合Matcont 軟件分析了EMFN 模型的平衡點分布與全局分岔性質, 發現并分析了該模型存在的亞臨界Hopf 分岔、隱藏放電及其周期放電與靜息態共存等現象.其次, 利用雙參數及單參數分岔、ISI 分岔和最大Lyapunov 指數等工具進行數值仿真, 詳細分析了EMFN 模型存在的伴有混沌及無混沌的加周期分岔結構、混合模式放電和共存模式放電等現象, 同時揭示了電場和磁場強度影響其放電節律的轉遷規律.最后,利用Washout 控制器將EMFN 模型的亞臨界Hopf 分岔轉化為超臨界Hopf 分岔, 使其在分岔點附近的拓撲結構發生改變, 由此達到消除其隱藏放電的目的.本文的研究結果證實了新建神經元模型具有豐富的放電節律, 將影響神經元的信息傳遞和編碼, 為完善神經元模型, 揭示電磁場對生物神經系統的影響, 以及探求一些神經性疾病的致病機理提供了思路.

1 引 言

神經系統是生物體內起主導作用的系統, 支配著感覺、運動和記憶等功能.生物神經系統的基本結構和功能單位是神經元, 神經元通過復雜的放電活動來實現非線性功能.神經元的放電活動具有一定的節律性, 這與神經信息的產生、轉遞和解碼密切相關, 因此, 研究神經元的放電活動的動力學特性具有重要的生理學和神經工程意義[1?4].神經元模型的建立為定量分析神經元放電活動的規律奠定了基礎, 同時也極大地促進了神經科學的發展[5?8].實際上, 神經元模型可以分為兩大類: 一類是基于電導依賴或生物物理學模型, 它強調了各離子通道對神經元膜電位的影響, 其中典型的模型有Hodgkin-Huxley(HH)模型[9]、Morris-Lecar(ML)模型[10]、Chay 模型[11]和Hindmarsh-Rose (HR)模型[12]等, 通常情況下, 這類模型的維數較高并且數值分析困難, 不適合應用于耦合系統和神經元網絡的研究; 另一類是現象學或者非電導依賴模型,這類模型注重神經元的輸入和輸出關系, 并不考慮具體的生理機制對膜電壓的影響, 例如Izhikevich模型[13]、Fitzhugh-Nagumo(FHN)模型[14]和True-North 模型[15]等, 這類模型結構簡單并且計算簡便, 從而被廣泛應用于神經元動力學分析及其網絡仿真的研究[16?18].

近些年, 一些新的元素, 如憶阻器、磁通和隱藏吸引子等被引入神經元模型, 使得這些常規的神經元模型展現出更豐富的動力學和放電特性.Babacan 等[19]提出了一種具有閾值切換機制的新型憶阻器仿真器, 并建立了一個能夠產生峰和簇放電行為的神經元電路.Bao 等[20]報道了一種帶有憶阻器的三維HR 神經元模型, 該模型不存在任何平衡點, 卻存在周期吸引子與混沌吸引子共存的隱藏動力學現象, 并且基于電路實驗進行了數值模擬驗證.Usha 和Subha[21]分析了帶有磁通變量的HR神經元模型分岔模式, 并研究了電突觸和化學突觸耦合下同步、去同步和尖峰放電態與靜息態共存振蕩行為.Zhao 等[22]討論了具有磁通耦合的擴展HR神經元的相位同步問題, 發現在相同的激勵電流下, 隨著耦合相位的增大, 耦合強度達到相位同步的閾值逐漸變小.Pham 等[23]研究了具有記憶性突觸權值的神經網絡, 它可以在不存在平衡點的情況下產生超混沌隱藏吸引子.

近年來研究發現, 神經系統中各離子的波動或者外界電磁場的變化都會引起神經元膜電位的變化, 同時, 神經元放電活動的變化也會引起細胞內外的電場和電磁場分布的改變.因此, 考慮電磁感應效應的神經元放電活動成為最近幾年被廣泛研究的神經元課題之一[24?28].Ma 和Tang[24]開創性地考慮在離子穿越細胞膜或持續的電磁輻射刺激下, 細胞膜內外會產生電磁感應而改變細胞膜上的磁通量, 進而在HR 神經元模型中引入磁通量去描述上述電磁感應現象, 并對其進行新角度下的放電行為初探.隨后, Lü等[25]討論了考慮磁通量的神經元模型的動力學特性, 通過改變初始狀態, 可以觀察到多種模式的放電活動, 反映了該神經元的記憶效應.Wu 等[26]通過對上述神經元施加加性相位噪聲來檢測神經元在膜態中的動態響應和相變,對電活動的動力學特性進行了檢測和討論, 并觀察到雙相干共振行為.此外, 在神經元模型的膜電位采樣時間序列中可以觀察到多種模式的電活動.Kafraja 等[27]提出了一個三變量具有記憶性的Izhikevich 模型來描述神經元在電磁感應和噪聲作用下的動力學行為, 該模型可以描述內外磁場對神經元的影響.An 和Zhang[28]基于雙參數數值仿真發現了考慮磁通量的神經元模型具有豐富的分岔結構, 例如混合模式振蕩和加周期分岔, 并設計一個Hamilton 能量反饋控制器, 將其控制到不同的周期簇放電狀態.上述文獻的研究表明, 電磁感應現象的引入, 使原有的神經元模型更具有現實和物理意義, 同時也能展現出更豐富的放電模式.此外,這些模型也可用于進一步對電磁輻射影響下的生物神經系統的研究.

神經元胞體內外包含大量的離子, 其跨膜運動會形成跨膜電流, 從而會引起膜電位的波動.同時,細胞膜表面可視為一個具有均勻分布電荷的帶電表面, 于是細胞內外可產生感應電場.另外, 當神經元處于外界電場下時, 細胞內部電場的分布也會改變.因此, 在膜電位的連續波動過程中, 電場分布的影響更需要考慮, 但現有神經元模型中只有僅少數考慮電場的影響.Ma 等[29]考慮電荷分布和極化變化引起的電場效應, 在HR 神經元的基礎上,建立了一個具有感應電場的神經元模型, 并對其電活動的動力學行為進行了分析和討論, 證實了在不同的電場作用下, 神經元的電活動會發生不同的模態躍遷.Du 等[30]研究了感應電場作用下正弦IEF和隨機相位對單個神經元和神經元網絡信息編碼和傳輸的聯合作用機制, 得到了單個神經元放電模式的躍遷狀態, 包括簇放電、峰放電和閾下振蕩,并提出應該采用多種方法改進或重建神經元系統電磁感應模型.另外, 通過電容器耦合, 非線性電路之間也可以建立時變的電場.Wang 等[31]綜述了憶阻器在神經動力學、時延、同步模式, 特別是電磁感應與輻射、場效應下神經元模型的建立、場耦合對非線性電路間信號傳播的影響.同時指出,兩個電路的輸出可以通過在極板上連續充電和放電, 激活耦合電容器中的時變電場.Oliveira 等[32]討論了兩個電容耦合Van der Pol 振蕩器的同步問題, 并通過數值仿真和理論分析指出此耦合方式不會增加額外的噪聲源和電力消耗.Xu 等[33]指出, 通過電容耦合的兩個非線性電路被激活時, 可以產生一個時變電場, 耦合電容器的能量流傳播可進一步調節輸出和動力學耦合電路.電場和電容器耦合為神經元之間的信號編碼提供了新的認識, 即使在突觸耦合被抑制或未激活的情況下, 電場效應也有利于神經元之間的信號傳播.

綜上, 這些研究與對傳統的離子通道電流對細胞膜電位影響的研究不同, 它們為神經元之間信號的編碼和傳播提供了新視角.本文認為, 在神經元細胞內外的電生理環境中會因帶電離子的傳輸而產生電場效應, 繼而激發出磁場.根據麥克斯韋電磁場理論, 變化的電場和磁場相互激發, 形成電磁場進而共同調節神經細胞的電活動.因此, 從物理學的角度同時考慮電場和磁場對神經元放電行為的影響更具有意義, 這對神經元模型的優化設計將具有重要的參考價值.

基于上述討論, 考慮到電場和磁場都能對神經元放電活動產生不可忽略的影響, 本文在三維HR神經元模型的基礎上, 同時引入了電場變量和磁通變量分別來描述電場和磁場共同作用下的HR 神經元模型, 改進了一個電磁場影響下的五維神經元模型(EMFN 模型).通過Matcont 軟件, 分析EMFN模型的全局分岔和第一Lyapunov 系數, 揭示了該模型存在亞臨界Hopf 分岔、雙穩態以及隱藏放電現象.利用雙參數分岔、單參數分岔、最大Lyapunov指數和時間響應圖等數值仿真, 發現在不同的參數平面上, 新模型普遍存在伴有混沌的加周期分岔、無混沌的加周期分岔、混合模式放電和共存模式放電等動力學行為.最后, 采用Washout 控制器對EMFN模型的Hopf 分岔類型進行控制, 在不改變平衡點位置的前提下, 使其亞臨界Hopf 分岔轉化為超臨界Hopf 分岔, 以達到消除隱藏放電的目的.本文對設計的EMFN 模型進行數學和物理學視角下的動力學行為及其放電特性初探, 以期為神經元的信息傳遞和編碼以及振動模態提供理論指導.

2 模型描述

馬軍研究組[24,25]認為神經元膜電位的持續波動或者改變將會在細胞周圍產生磁場, 由此通過憶阻器引入磁通變量來模擬磁場與膜電位之間的相互作用, 建立了一個四變量的HR 神經元模型:

式中, 變量 x , y, z, φ, I 分別表示膜電位、與鈉離子和鉀離子相關的恢復快電流、與鈣離子相關的自適應慢電流、磁通變量和外界刺激電流; a, b, c, d,r, s, χ0是模型重要的參數; Ie表示磁場產生的電磁感應電流, 其數學表達式為Ie=?k0W(φ)x=?k0(α+3βφ2)x , W (φ) 為磁控憶阻器, 參數 k0調節膜電位和磁通之間的相互作用, 可視為磁場的反饋強度.由于該模型引入了磁通變量, 使該模型具有更多的分岔參數, 并且膜電位的放電模式將更加豐富, 從而被廣泛應用于神經元的分岔分析、膜電位的遷移控制和場耦合下神經網絡的同步行為研究[34,35].眾所周知, 神經元的放電活動會受各離子通道作用的影響, 離子通道能改變細胞內外離子濃度的分布, 同時也決定了細胞的動作電位.Ma 等[29]認為外電場參與了細胞內各離子的傳輸和交換, 從而引起膜電位的變化, 通過引入電場變量, 詳細研究了電場作用下HR 神經元模型的放電特性.

目前, 大多數神經元模型主要關注膜內的離子通道電流.實際上, 細胞膜上離子的反復交換所產生的電場和磁場效應會對膜電位的變化和帶電離子的交換有一定的反饋作用.因此, 根據麥克斯韋電磁場理論, 從物理的角度對電磁場影響下的神經元模型進行研究很有必要.上述的研究啟發我們引入電場變量和磁通變量分別刻畫感應電場和感應磁場, 建立如下的電磁場影響下的HR 神經元模型(EMFN 模型):

(2)式中的變量和參數與(1)式中相同.變量 E 表示感應電場, 考慮到與自適應慢電流 z 相比, 快電流 y 對電場的變化更敏感, 由此通過對變量 y 施加k1E 來描述電場的影響, k1用來調節電場和快電流之間的相互作用, 可視為電場的反饋強度.對于神經元系統而言, 參數 k1, k2, k4的取值很大程度上取決于神經元的固有的興奮特性.此外 ? k3φ 和?k5E 分別表示神經元對磁場和電場的自適應調節項.最初的HR 模型是一個無量綱的模型, 因此,模型(2)的參數均為無量綱參數.這里, 部分參數取Hindmarsh 和Rose 在1984 年建立HR 神經元模型時設定的基準值: a =1.0 , b =3.0 , c =1.0 ,d=5.0 , s =4.0 , r =0.006 , χ0=?1.61.I 一般 在[?10,10] 內取值, 本文取 I =3.其余參數取為: α=0.2, β =0.03 , k0=0.1 , k1=0.1 , k2=0.3 , k3=0.5 ,k4=0.2 , k5=0.3.EMFN 模型在神經元模型(1)的基礎上引入電場作為第五個變量去描述神經元內部的電場效應, 以此來研究神經元內部電場和磁場相互激發而形成的電磁場對神經系統的影響.

3 分岔分析與隱藏放電

3.1 平衡點穩定性與Hopf 分岔分析

考慮電場和磁場效應的EMFN 模型(2)將具有更加豐富的動力學特性, 為了探索其放電和分岔模式, 首先需要分析其平衡點的分布及其穩定性.由模型(2)的動力學方程可得零線方程組:

通過計算方程組(3)的實數解(x0,y0,z0,φ0,z0), 可獲得模型(2)平衡點為P0=(x0,y0,z0,φ0,z0), 并且由方程組(3)可得

當各參數取基準參數值時, 方程(4)的根依賴于 I 和 ki(i=1,2,3,4,5) 的取值, 并且其實數根就是平衡點處的 x0.由方程組(3)可確定平衡點處其余的坐標:

模型(2)在平衡點 P0的穩定性是由相應線性化矩陣 A 的特征根所決定:

通過計算線性化矩陣的特征根, 不僅可以直接得到系統的穩定性, 而且也可以分析其分岔規律.這里, 使用Matcont 軟件去避免復雜繁瑣的計算,對于連續系統而言, 它是強大分岔分析工具.由此通過Matcont 編程分析可得到如圖1 所示的EMFN模型(2)關于參數 I 和膜電壓 x 的全局分岔圖(圖中黑色曲線為平衡點處的 x0, 紅色星點為Hopf 分岔點 H , 綠色區域為穩定的靜息態, 紅色區域為周期2 簇放電態, 黃色區域為周期1 尖峰放電態, 淺藍色區域為小振幅振蕩發散態), 并且在Hopf 分岔點 H 處的位置和相應的特征根為

此外, 利用Matcont 還能計算出相應的第一Lyapunov 系數 LH=0.00024971>0.由此可知,EMFN 模型(2)在分岔點 H 發生亞臨界Hopf 分岔, 其分岔的方向為參數 I 減小的方向.通過觀察圖1 可知, EMFN 模型(2)在Hopf 分岔點 H 分岔前后都存在放電模式不同的吸引域, 即在分岔點H附近, 模型(2)的放電模式不僅與參數 I 取值相關,而且與初值 x 也密切相關.因此, 揭示神經元模型(2)在分岔點 H 附近的放電規律, 可為對神經元異常放電行為的探討提供參考.

圖1 EMFN 模型(2)在 x ∈[-20,20] 時的吸引域Fig.1.The attractive basins of EMFN model (2) when x ∈[-20,20].

3.2 共存放電

在非線性系統參數不變的情況下, 改變初始狀態, 系統運行可能漸近趨向的定點、混沌或周期、準周期等不同穩定狀態的現象稱為多穩定性或共存吸引子[36,37], 多穩定性表現出非線性動力系統豐富的穩定狀態多樣性, 并為系統提供了很大的靈活性.本文所設計的EMFN 模型(2)具有豐富的共存放電現象.

當I=I1=1.172>IH, 相應的平衡點為

則矩陣(6)對應的特征根為

明顯地, 平衡點 P1為不穩定的焦結點.若 I1不變, 取初值為

EMFN 模型(2)的放電波形圖如圖2(a)中藍色曲線所示, 圖2(b)為圖2(a)中藍色曲線的放大.易知, 模型(2)在6000 個時間單位內處于發散趨勢的小振幅振蕩態(經過更長的時間仿真, 此發散趨勢的振蕩態將趨于周期2 簇放電態, 這里將不再展示).其原因是平衡點 P1的特征根0<故當初值在 P1附近取值時, 神經元呈現出圖2(b)所示的發散態.明顯地, 神經元模型(2)在平衡點 P1附近的小振幅振蕩區域比較小,而周期2 簇放電區域比較大.

圖2 EMFN 模型(2)的時間響應圖和吸引子 (a) P 1 附近的共存振蕩; (b) 圖(a)中藍色放大圖; (c) P 1 附近的周期2 極限環;(d) P 2 附近的共存振蕩; (e) 圖(d)中藍色放大圖; (f) P 2 附近的周期2 極限環; (g) P 3 附近的共存振蕩; (h) 圖(g)中藍色放大圖;(i) P 3 附近的周期1 極限環Fig.2.Time responses and attractors of EMFN model (2): (a) Coexistence oscillation near P 1 ; (b) enlargement of the blue in (a);(c) limit cycle with period-2 near P 1 ; (d) coexistence oscillation near P2 ; (e) enlargement of the blue in (d); (f) limit cycle with period-2 near P 2 ; (g) coexistence oscillation near P 3 ; (h) enlargement of the blue in (g); (i) limit cycle with period-1 near P 3.

當初值在離 P1較遠處取值時, 如

EMFN 模型(2)的放電波形圖如圖2(a)中紅色曲線所示的周期2 簇放電狀態, 其振蕩幅度遠大于處在發散趨勢的藍色曲線所示的振蕩幅度, 圖2(c)為對應的極限環吸引子.

當外界刺激電流選取如下的不同值時, 會出現不同振幅的共存振蕩現象.

1)I =I2=1.152

相應的平衡點為

則矩陣(6)對應的特征根為

明顯地, 平衡點 P2為穩定的焦結點.若 I2不變, 取

EMFN 模型(2)的放電波形圖如圖2(d), (e)中藍色曲線所示, 圖2(e)為圖2(d)中藍色曲線的放大.易知, 神經元模型此時處于收斂趨勢的小振幅振蕩態.其原因是平衡點 P2的特征根故當初值在 P2附近取值時, 模型(2)呈現出如圖2(d),(e)中所示的靜息態.

當初值在離 P2較遠處取值時, 模型(2)也表現出初值敏感特性:

EMFN 模型(2)的放電波形圖如圖2(d)中紅色曲線所示的周期2 簇放電狀態, 其相軌跡為穩定的周期2 極限環吸引子, 如圖2(f)所示.明顯地, 模型(2)在兩組初值下呈現小振幅振蕩和周期2 簇放電的共存狀態.

2)I =I3=1.086

相應的平衡點為

則矩陣(6)對應的特征根為

明顯地, 平衡點 P3為穩定的焦結點.若 I3不變,取初值

EMFN 模型(2)的放電波形圖如圖2(g)和圖2(h)中藍色曲線所示, 圖2(h)為圖2(g)中藍色曲線的放大.可以看出, 模型(2)此時處于靜息態.對比圖2(e)和圖2(h)可知, I =I3時的膜電壓的收斂速度遠比 I =I2時要快, 其原因是平衡點 P3距離亞臨界Hopf 分岔點 H 比平衡點 P2較遠, 并且

當初值取

時, EMFN 模型(2)的放電波形圖呈現如圖2(g)中紅色曲線所示的周期1 峰放電狀態, 其相軌跡為穩定的周期1 極限環吸引子, 如圖2(i)所示.此時,神經元的雙穩態吸引域如圖3(e)和圖3(f)所示(紅色區域為周期2 簇放電, 黃色區域為周期1 峰放電, 淺藍色區域為具有收斂趨勢的小振幅振蕩,綠色區域為靜息態, 黑色紅星為平衡點), 相比圖3(c)和圖3(d)所示的吸引域, 其拓撲結構發生了顯著的變化, 即代表靜息態的綠色區域變大, 而代表周期1 峰放電態的黃色區域比處于周期2 簇放電的紅色區域小.明顯地, EMFN 模型(2)在兩組初值下呈現靜息態和周期1 峰放電的共存狀態.

通過上述分析, EMFN 模型(2)在不同的初始狀態下發現不同吸引子共存的現象, 說明該系統對初始狀態具有很強的敏感性, 這與磁控憶阻器的引入而導致的神經元模型非線性加強有關.

3.3 隱藏振蕩

從上述討論可知, 當外界刺激電流取不同的值時, EMFN 模型(2)在不同的初始狀態下展現不同的放電狀態, 如靜息態、周期1 放電和周期2 放電等共存的現象.另外, 還可以發現EMFN 模型(2)具有另一個有趣的現象, 即隱藏吸引子.隱藏吸引子的吸引域不包含平衡點的鄰域[38](無平衡點或者僅有穩定平衡點的動力系統產生的吸引子都屬于隱藏吸引子).

圖3 EMFN 模型(2)的吸引域 (a), (c)和(e)是外界刺激電流分別取 I 1, I2, I3 時x-y 上的吸引域; (b), (d)和(f)是外界刺激電流分別取 I 1, I2, I3 時x-φ 上的吸引域Fig.3.The attractive basins of EMFN model (2): (a), (c) and (e) are the attractive basins in x-y plane under I 1, I2, I3 , respectively;(b), (d) and (f) are the attractive basins in x-φ plane under I 1, I2, I3 , respectively.

當 I =I1時, 通過分析發現EMFN 模型(2)在其對應的不穩定平衡點 P1的自激振蕩作用下產生周期2 的簇放電現象.從圖3(a)和圖3(b)的紅色區域可以看出, 周期2 簇放電的吸引域與平衡點的小領域不相交.明顯地, 此時的周期2 簇放電屬于隱藏放電范疇.

當 I =I2和 I3時, EMFN 模型(2)對應的平衡點為穩定的焦結點, 產生的周期2 簇放電現象則屬于隱藏放電范疇.其機制為模型(2)在分岔點 H 處發生了亞臨界Hopf 分岔, 相應的平衡點的穩定性由不穩定的放電狀態變為穩定的放電狀態, 并產生了不穩定的極限環, 同時, 在這個不穩定的極限環外存在一個穩定的周期2 極限環吸引子.因此, 當參數取值不變, 在穩定的周期2 極限環所產生吸引域內取值時, 模型(2)處于周期為2 的簇放電狀態,如圖3(c),(d)直觀所示, 在x-y 和x-φ平面上的吸引域上都存在著大范圍的隱藏周期2放電區域.

4 基于雙參數的分岔分析

4.1 伴有混沌和無混沌的加周期分岔

EMFN 模型(2)能比較精確地描述電磁場下神經元的放電活動, 并且具有更多的分岔參數.當外界電磁場變化時, 很難保持新系統的各參數取值不變, 因此, 研究多個參數同時變化下神經元模型的放電節律轉遷更具有實際參考價值.本節通過計算雙參數分岔圖、單參數分岔圖、最大Lyapunov指數和時間響應圖等來模擬EMFN 模型(2)隨各參數組合下的分岔結構.此外, 在數值模擬過程中,采用四階龍格-庫塔算法, 初值固定為

圖4 EMFN 模型(2)關于x 的雙參數分岔圖 (a) I 和 b 對應的分岔圖; (b) I 和 d 對應的分岔圖; (c) c 和 d 對應的分岔圖;(d) I 和 r 對應的分岔圖; (e) s 和 r 對應的分岔圖; (f) χ 0 和 r 對應的分岔圖; (g) I 和 k1對應的分岔圖; (h) I 和 χ 0 對應的分岔圖;(i) I 和 s 對應的分岔圖Fig.4.Two-parameter bifurcation diagrams of EMFN model (2) versus x: (a) Bifurcation diagram versus I and b ; (b) bifurcation diagram versus I and d ; (c) bifurcation diagram versus c and d ; (d) bifurcation diagram versus I and r ; (e) bifurcation diagram versus s and r ; (f) bifurcation diagram versus χ 0 and r ; (g) bifurcation diagram versus I and k1 ; (h) bifurcation diagram versus I and χ 0 ; (i) bifurcation diagram versus I and s.

需要說明的是, 為了避免短暫的振蕩行為, 舍棄前2 × 105次迭代, 用后續的3 × 107步計算最大Lyapunov 指數.隨后繼續進行20 × 105步的迭代, 并且記錄迭代過程中膜電壓x 存在的極值(最大值或最小值), 可通過判斷各組脈沖是否重復,由此獲得周期態或者混沌態[39?41].EMFN 模型(2)的膜電壓x 在雙參數平面中的分岔圖如圖4 所示,其仿真圖采用了 5 00×500 等距參數點覆蓋, 此外,用不同的顏色繪制膜電壓的不同放電狀態, 并且圖右側的顏色應運用相應的數字進行標記.需要說明的是, 當膜電壓x 的周期大于等于20 時, 統一運用白色區域表示.圖4(a)—(f)展示了EMFN 模型(2)存在的“梳”狀混沌結構.如圖4(a)所示, 在參數I ∈[2.1,4.1], b ∈[2.8,3.5]的參數平面上, 模型(2)存在豐富的周期放電和混沌放電狀態.觀察圖4(a)中最大的顏色區域(即黃色區域), 膜電壓x 處于周期1 尖峰放電態, 并且通過倍周期分岔結構與混沌區域相連接.類似地, 其余周期區域也存在倍周期分岔結構, 并且以“舌形”嵌入混沌區域中[41].此外, 通過觀察發現許多自相似和遞歸性的分岔模式.當保持參數 b =?0.3485I +4.2318 不變, 參數I ∈[2.1,4.1]時(沿著圖4(a)中黑色直線所示的方向), 模型(2)的峰峰間期(ISI)分岔圖和膜電壓x 的最大值分岔圖分別如圖5(a)和圖5(b)所示,圖6 為圖5 相應的最大Lyapunov 指數.因此,EMFN 模型(2)表現出伴有混沌的加周期分岔模式, 即周期1 → 周期2 → 周期4 → ···→ 首次出現混沌窗口 → 周期2 → 周期4 → 周期8 → ···→ 再次出現混 沌 窗口 → 周期3 → 周期6 → 周期12 → ···.模型(2)按照上述的加周期分岔和混沌交替模式,并且每經過一次混沌放電后的周期數要比之前相應的周期數大1, 最終進入混沌放電態或者更高周期的簇放電態.在圖4(a)中, A 點所對應的參數點為 A (I,b)=(2.389,3.293) , 此時的模型(2)處于周期3 簇放電狀態, 其膜電壓x 的時間響應如圖7(a)所示.此外, 在圖4(a)中的B, C 和 D 點, 相應的參數取值分別為

圖5 關于 I 的ISI 分岔圖和單參分岔圖 (a) ISI 分岔圖; (b) 單參分岔圖Fig.5.ISI bifurcation and one-parameter bifurcation versus I : (a) ISI bifurcation; (b) one-parameter bifurcation.

圖6 對應于圖5 的最大 Lyapunov 指數圖Fig.6.The largest Lyapunov diagram corresponding to Fig.5.

時, 膜電壓x 相應的時間響應分別如圖7(b)—(d)所示.圖4(b)—(f)也存在類似的伴有混沌的加周期分岔結構, 本節不再詳細描述.

圖4(g)—(i)顯示了EMFN 模型(2)的無混沌加周期分岔結構, 即通過加周期分岔模式從某種周期狀態直接進入相鄰的周期狀態.當參數I ∈[1.5,3], k1∈[0,0.28]時, 如圖4(g)所示, 沿著左下到右上的方向, 模型(2)的分岔模式為周期1 → 周期2 → 周期3 → ···→ 周期18 → 更高的加周期分岔模式.此外, 還可以觀察到這些周期區域呈現帶條狀分布, 并且隨著周期數逐漸增大, 相應的顏色帶逐漸減小(例如, 周期2 的顏色帶明顯比周期3的顏色帶寬).類似地, 在圖4(h)和圖4(i)所示的參數平面上, 模型(2)也存在無混沌加周期分岔結構.

通過上述討論可知, EMFN 模型(2)在不同的參數區域內存在不同的加周期分岔.神經元的節律模式與控制參數呈現出明顯的對應關系, 因此, 在圖4(a)—(f)中的參數變化區域內出現伴有混沌的加周期分岔現象的產生機理是模型(2)在此組參數的變化下具有更高的興奮性和更強的節律性, 進而引起了激變現象, 即簇放電經合并激變進入混沌放電, 再經邊界激變進入更高周期的簇放電, 同時伴有內部激變現象.在圖4(g)—(i)中的參數變化區域內, 新建模型僅存在加周期分岔現象, 未見混沌放電存在, 呈現出比較規則的加周期分岔轉遷模式.

4.2 混合模式放電

混合模式振蕩是一類大振幅與小振幅交替出現復雜的動力學現象, 能夠產生混合模式振蕩的模型通常是包含了多重時間尺度的非線性微分方程組[42], 并且在電子電路、化學反應和生物數學等領域都有發現[43,44].神經元模型屬于多時間尺度的動力系統, 它不僅能產生 L0型簇放電(即周期為 L 的簇放電模式), 也能夠產生 LS型混合模式放電(其中 L ,S 分別表示一個振蕩周期內大振幅與小振幅振蕩出現的次數).考慮到EMFN 模型(2)具有更多的分岔參數, 因此該模型存在豐富的簇放電模式和分岔結構.如圖8 所示, 在參數k0∈[0,1.2],d ∈[4.2,5.4]平面上, 模型(2)存在大量的帶狀周期區域和一些無規則的混合窗口, 并且相應的周期都是通過倍周期分岔結構與混沌區域相連接.此外,對于連續系統(2)而言, 在圖中右側區域的周期顏色帶存在不連續的情況, 這些異常的分岔結構往往是由于系統分岔類型發生了突變, 并且常伴隨著共存吸引子的產生[39,41].能夠產生混合模式放電的最小神經元模型是包含多重時間尺度的三維非線性方程組, 對于高維的動力系統, 由于產生機理的相關理論仍不完備, 于是采用數值方法對模型產生的混合模式放電進行分析[45].

圖7 EMFN 模型(2)關于 I 和 b 的時間響應圖 (a) I =2.389,b=3.239 時的周期3 簇放電; (b) I =2.577,b=3.173 時的周期4 簇放電; (c) I =2.733,b=3.134 時的周期5 簇放電; (d) I =2.898,b=3.093 時的周期6 簇放電Fig.7.Time response diagram of EMFN model (2) versus I and b : (a) Bursting with period-3 when I =2.389,b=3.239 ;(b) bursting with period-4 when I =2.577,b=3.173 ; (c) bursting with period-5 when I =2.733,b=3.134 ; (d) bursting with period-6 when I =2.898,b=3.093.

圖8 EMFN 模型(2)關于 k0 和 d 的雙參數分岔圖Fig.8.Two-parameter bifurcation diagram of EMFN model(2) corresponding to k0 and d.

通過數值模擬發現, EMFN 模型(2)在圖8中左側的周期區域處于周期簇放電狀態, 但在右側區域處于混合模式放電狀態, 而在周期不連續區域,模型(2)處于周期簇放電與相應的混合模式放電共存狀態.例如, 若參數 ( k0,d)=(0.6587,4.2596) 保持不變, 當初值分別為 ( 9.0,0.1,0.1,0.1,0.1) 和(?9.0,0.1,0.1,0.1,0.1)時, 模型(2) 分別處于周期2簇放電和 21型混合模式放電, 其相應的時間響應圖分別如圖9(a)和圖9(b)所示.圖10(a)和圖10(b)展示的是x-φ 和y-E 平面上的吸引域, 其中黃色區域為周期2 簇放電, 紅色區域為 21型混合模式放電.由此可知, 模型(2)的放電模式不僅與其參數取值相關, 而且也與初值密切相關.通過觀察圖9(a)和圖9(b), 可以發現混合模式放電的大振幅周期數與簇放電的周期數相同, 即模型(2)可以通過初值的改變產生小振幅放電, 并且與原簇放電模式結合而形成新的混合放電模式.類似地, 當 k0和 d 分別取值 ( 0.7032,4.3758) 和 ( 0.7123,4.5032) 時, 模型(2)分別處于周期3 簇放電與周期 31共存模式放電狀態和周期4 簇放電與周期 42共存模式放電狀態,圖10(c)—(f)為相應平面上的吸引域(紅色區域表示周期簇放電, 黃色區域表示混合模式放電).此外, 在圖8 中不連續的分岔區域還存在著大量的共存模式放電, 這些區域對初值十分敏感, 對機體神經系統來說將可能會引起突發的心率失調, 因此在相關生物實驗中應盡量避免這些區域, 我們可以通過調節電磁反饋項 k0, 使神經元趨于穩定的周期簇放電狀態.

圖9 EMFN 模型(2)的時間響應圖 (a) 周期2 放電; (b) 周期 21 放電; (c) 周期3 放電; (d) 周期 31 放電; (e) 周期4 放電; (f) 周期 42 放電Fig.9.Time response diagrams of EMFN model (2): (a) Discharge with period-2; (b) discharge with period- 21 ; (c) discharge with period-3; (d) discharge with period- 31 ; (e) discharge with period-4; (f) discharge with period- 42.

當保持參數 d =1.9967k0?4.2091 不變, 并且k0∈[0,0.6]時, 即沿著圖8 中黑線所示的分岔方向,EMFN 模型(2)的ISI 分岔圖和關于膜電壓 x 的分岔圖分別如圖11(a)和圖11(b)直觀所示.模型(2)首先進行伴有混沌的加周期分岔模式: 周期2 → 混 沌 → 周 期3 → 混 沌 → 周 期4 → 混 沌 → 周 期5 → 混沌, 此后混沌窗口消失, 進入無混沌加周期分岔模式: 周期6 → 周期7 → ···→ 周期14 → 更高的簇放電狀態.圖12 展示了相應的最大Lyapunov指數, 很好地吻合了圖11 所示的分岔圖.此外, 在圖8 右側混合模式放電趨于白線所示的分岔方向,即當參數 d =2.9851k0?1.8119,k0∈[0.8,1.2] 時膜電壓 x 的發放數[43,44](在一個混合模式放電周期內,小振幅數占總振幅數的比例)的變化如圖13 所示.明顯地, 膜電壓 x 的發放數呈現階梯狀分布, 其變化規 律 為: 周 期 12→ 周 期 22→ 周 期 32→ 周 期42→周期 43→ 周期 53→ 周 期 63→ 周 期 73→ 周 期83→周 期 84→ 周 期 94→ 周 期 1 04→ 周 期 1 14→ 周 期124→ 周期 1 34→ 周期 1 35→ 周期 1 45→ 周期165→周期 1 75.雖然模型(2)的混合模式放電周期較高并且比較復雜, 但其變化具有一定的自相似規律,即隨著 k0的增加, 大振幅和小振幅都會隨之逐漸增加, 如圖13 所示的發放數呈現“四階梯”狀分布.

總之, EMFN 模型(2)普遍存在豐富且復雜的混合模式放電和共存模式放電, 這是由于高維系統的復雜性引起的, 本節的研究結果揭示了其分岔規律, 這將為神經元相關疾病的治療和控制提供有益的探討.

圖10 EMFN 模 型(2)的共存吸引域 (a) ( k0,d)=(0.6587,4.2596) 時 關 于x-φ 的 吸 引 域; (b) ( k0,d)=(0.6587,4.2596) 時 關于x-E 的吸引域; (c) ( k0,d)=(0.7032,4.3758) 時關于x-φ 的吸引域; (d) ( k0,d)=(0.7032,4.3758) 時關x-E 于的吸引域; (e) (k0,d)=(0.7123,4.5032) 時關于x-φ 的吸引域; (f) ( k0,d)=(0.7123,4.5032) 時關于x-E 的吸引域Fig.10.The coexisting attraction domains of EMFN model (2): (a) Attractive basins of x-φ plane when ( k0,d)=(0.6587,4.2596) ;(b) attractive basins of x-E plane when ( k0,d)=(0.6587,4.2596) ; (c) attractive basins of x-φ plane when (k0,d)=(0.7032,4.3758) ; (d) attractive basins of x-E plane when ( k0,d)=(0.7032,4.3758) ; (e) attractive basins of x-φ plane when(k0,d)=(0.7123,4.5032) ; (f) attractive basins of x-E plane when ( k0,d)=(0.7123,4.5032).

圖11 關于 k0 的ISI 分岔圖和單參分岔圖 (a) ISI 分岔圖; (b) 單參分岔圖Fig.11.ISI bifurcation and one-parameter bifurcation versus k0 : (a) ISI bifurcation; (b) one-parameter bifurcation.

圖12 對應于圖11 的最大 Lyapunov 指數圖Fig.12.The largest Lyapunov diagram corresponding to Fig.11.

圖13 膜電壓的發放數關于參數 k0 變化圖Fig.13.The change of spike count of membrane voltage versus parameter k0.

4.3 電場和磁場強度影響下的混合模式放電

與雙參數 ( k0,d) 影響下的動力學類似, 電場和磁場強度 ( k1,k0) 對EMFN 模型(2)也有很大的影響, 尤其是周期簇放電的產生與轉遷, 圖14 展示了 k1∈(0,0.25),k0∈(0,1.2) 時的雙參分岔圖.明顯地, 模型(2)在此組參數下存在大量的帶狀周期區域和一些無規則的混合窗口, 并且相應的周期也均是通過倍周期分岔結構與混沌區域相連接.

若EMFN 模型(2)的初值(0.1,0.1,0.1,0.1,0.1)保持不變, 參數為 ( k1,k0)=(0.01403,0.7214) 和(k1,k0)=(0.01603,0.7743)時, 模型(2)分別處于如圖15(a)和圖15(b)所示的周期6 和周期 63的混合模式放電狀態.類似地, 當 ( k1,k0) 分別取值(0.03808,0.7599) 和 ( 0.03809,0.8056) 時, 模型(2)分別處于如圖15(c)和圖15(d)所示的周期7 簇放電與周期73的混合模式放電狀態.當 ( k1,k0) 分別取值(0.06112,0.7841) 和 ( 0.05912,0.8441) 時, 模型(2)分別處于如圖15(e)和圖15(f)所示周期8 簇放電與周期 83的混合模式放電狀態.

圖14 EMFN 模型(2)關于 k0 和 k1 的雙參數分岔圖Fig.14.Two-parameter bifurcation diagram of EMFN model(2) corresponding to k0 and k1.

圖15 EMFN 模型(2)的時間響應圖 (a) 周期6 放電; (b) 周期 63 放 電; (c) 周期7 放電; (d) 周期 73 放 電; (e) 周期8 放電;(f) 周期 83 放電Fig.15.Time response diagrams of EMFN model (2): (a) Discharge with period-6; (b) discharge with period- 63 ; (c) discharge with period-7; (d) discharge with period- 73 ; (e) discharge with period-8; (f) discharge with period- 83.

圖16 關于 k0 的 ISI 分岔圖和單參分岔圖Fig.16.ISI bifurcation and one-parameter bifurcation versus k0.

從圖15 可以發現, 混合模式放電的大振幅周期數與簇放電的周期數相同, 即可以通過微調電場和磁場強度使模型(2)產生小振幅放電, 并且與原簇放電模式結合而形成新的混合放電模式.

當保持參數 k1=0.8333k0?0.3333 不變, 并且k0∈[0.4,0.7]時, EMFN 模型(2)的ISI 分岔圖和關于膜電壓 x 的分岔圖分別如圖16(a)和圖16(b)直觀所示.模型(2)展現出無混沌的加周期分岔模式: 周期5 → 周期6 → 周期7 → 周期8 → 周期9→···.此外, 在圖14 右側混合模式放電區域, 即當參數 k1=0.4762k0?0.3714(0.8 ≤k0≤1.2) 時膜電壓 x 的發放數變化如圖17 所示.明顯地, 膜電壓x的發放數呈現階梯狀分布, 其變化規律為: 周期63→ 周期 73→ 周 期 83→ 周期 93→ 周期 1 03→ 周期 113→ 周期 114→ 周 期 124→ 周期 134→ 周期144→ 周期 154→ 周期 164.

圖17 膜電壓的發放數關于參數 k0 變化圖Fig.17.The change of spike count of membrane voltage versus parameter k0.

上述混合模式放電周期多變且復雜, 但仍有一定的自相似規律, 即隨著 k0的增加, 大振幅和小振幅都會隨之逐漸增加, 如圖17 所示的發放數呈現“雙階梯”狀分布.本節揭示了電場和磁場影響下EMFN 模型(2)中的混合模式放電節律轉遷規律,以期有助于探尋電磁場對生物系統興奮性的影響機理.

4.4 電場方程對EMFN 模型(2)的影響

本文將第五維電場變量引入文獻[24]中考慮磁場時的神經元模型, (2)式中 k4和 k5是兩個關鍵參數, 體現了外界刺激電流對神經元內部電場的感應和反饋效應.圖18 描述了EMFN 模型(2)關于I和 k4, k5的雙參數分岔圖.

從圖18 中可以看出, 在參數I ∈[1,4], k4∈[0,1]平面上, EMFN 模型(2)存在周期和混沌振蕩區域.其中, 周期振蕩區域呈現“條狀”分布, 相應的周期數越大, 其“條形”區域越窄.此外, 在每個“條狀”周期的末端將通過倍周期分岔結構與混沌區域相連接, 并且, 圖中右上角存在大范圍白色區域,說明該系統在此參數區域內存在高周期簇放電狀態或者混沌態.類似地, 在參數I ∈[1,4],k5∈[0,1]平面上, EMFN 模型(2)呈現出有規律的周期分布, 即: 沿著由左上到右下的方向, 該模型將由靜息態 → 周期1 → 周期2 → ··· , 到高周期的簇放電態或混沌態.此外, 與混沌相接的周期區域呈現出“舌狀”分布, 并且具有明顯的倍周期分岔結構.圖19和圖20 描述了EMFN 模型(2) 關于 I 和 k4, k5的ISI 分岔圖和單參數分岔圖.

圖18 EMFN 模型(2)關于I 和 k 4, k5 的雙參數分岔圖Fig.18.Two-parameter bifurcation diagram of EMFN model (2) corresponding to I and k 4, k5.

圖19 當 I ∈[3.1, 4], k4 =0.4556I -1.4122 時, EMFN 模型(2)的 ISI 分岔圖(a)和單參分岔圖(b)Fig.19.(a) ISI bifurcation and (b) one-parameter bifurcation of EMFN model (2) when I ∈[3.1, 4], k4 =0.4556I -1.4122.

圖20 當 I ∈[2.9, 3.7], k5 =-0.9659I +3.7897 時, EMFN 模型(2)的 ISI 分岔圖(a)和單參分岔圖(b)Fig.20.(a) ISI bifurcation and (b) one-parameter bifurcation of EMFN model (2) when I ∈[2.9, 3.7], k5 =-0.9659I +3.7897.

5 隱藏動力學控制

隱藏吸引子的存在可能導致突發混沌或周期大幅度躍遷, 將會影響系統的正常運行而帶來災難性的后果[38,46].因此, 如何有效地抑制神經元系統中的隱藏放電具有重要的理論和現實意義.由于EMFN 模型(2)在分岔點 H 處發生亞臨界Hopf分岔, 并且在其附近發現周期1 和周期2 的隱藏極限環吸引子.為了有效消除隱藏放電現象, 本節基于Washout 控制器, 對分岔點 H 的Hopf 分岔類型進行控制, 在不改變平衡點位置的前提下, 使其亞臨界Hopf 分岔轉化為超臨界Hopf 分岔, 從而使EMFN 模型(2)在分岔點 H 處附近的拓撲結構發生改變, 由此達到消除其隱藏放電的目的.施加控制后的模型為

式中, v 為Washout 濾波器的狀態變量, m 為控制器的反饋增益, ξ 為濾波器時間常數的倒數, 本文選取參數 ξ =0.07 , 其余的參數與EMFN 模型(2)的基準參數取值相同.

把受控系統(7)寫出矩陣形式, 即

在(8)式中有

系統(7)關于外界刺激電流 I 分岔的位置是不變的, 即當參數 I =IH時, 受控系統也發生亞臨界Hopf分岔, 此時相應的平衡點為

此時, 系統(7)在平衡點 Pc處的線性化矩陣 Ac為

由此可得矩陣 Ac的特征根為

從而可知, 受控系統(7)在平衡點 Pc處存在一對實部為零的共軛特征根, 即該系統發生了Hopf分岔,下面分析系統(7)在平衡點 Pc處關于參數 m 的Hopf 分岔類型.

將控制系統(7)中的線性項與非線性項分開表示為

(8)式中 G (X)=F(X)?AcX 為非線性項,對系統(8)做線性變換 X =BY +P′c, 其中P′c是 Pc的轉置矩陣.矩 陣B =[Re(μ1),Im(μ1),μ3,μ4,μ5,μ6]的 前兩列分別為特征根λ1=ω0i=0.03230434i 對應的特征向量 μ1的實部和虛部, 第三列為特征根 λ3對應的特征向量 μ3, 第四列為特征根 λ4對應的特征向量 μ4, 第五列為特征根 λ5對應的特征向量 μ5, 第六列為特征根 λ6對應的特征向量 μ6.由此經計算可得矩陣 B 及其逆矩陣B?1分別為

變換后的系統(8)如下:

(10)式中 J =B?1AcB 是系統 的Jordan 矩陣,=B?1F(BY +P′c)?B?1AcBY 為非線性項, 其分量形式為

根據Hopf 分岔理論, 可知穩定性指標 η 表達式為

穩定性指標 η 決定Hopf 分岔周期解的穩定性, 當η <0時, 受控系統(10)分岔產生的極限環是穩定的, 即系統發生超臨界Hopf 分岔; 當 η >0 時, 受控系統(10)分岔產生的極限環是不穩定的, 即系統發生亞臨界Hopf 分岔, 其判別式(11)中各特征量計算表達式分別為

數值計 算 可得 η =0.00007054m+0.00249712 ,當 η <0 時, 即當控制參數m

由3.3 節可知, 當外界刺激電流 I

反饋增益 m 對受控系統(7)放電模式的影響如圖21(a)所示(紅色區域表示周期2 的隱藏放電,黃色區域表示周期1 的隱藏放電, 綠色區域表示穩定平衡點的吸引域, 黑色星點分別表示平衡點Pc1和 Pc2).當膜電壓 x ∈[?20,20] , 反饋增益m

圖21 反饋增益 m 對受控系統(7)的放電影響 (a) 當I =I2 時, 受控系統(7)放電演化圖; (b) 當 I =I3 時, 受控系統(7)放電演化圖Fig.21.The discharge influence of feedback gain m to controlled system (7): (a) Discharge evolution of the controlled system (7) when I =I2 ; (b) Discharge evolution of the controlled system (7) when I =I3.

反饋增益 m 對受控系統(7)的放電模式的影響如圖21(b)所示.通過數值仿真發現, 當反饋增益m

6 結 論

在神經細胞中, 穿過細胞膜的離子濃度的波動會產生時變的感應電場和感應磁場, 本文基于麥克斯韋電磁場理論, 討論了電磁場影響下的HR 神經元模型的放電節律轉遷及其隱藏放電控制.主要工作為: 1) 發現EMFN 模型存在亞臨界Hopf 分岔、隱藏放電及其周期放電與靜息態共存的現象; 2) 通過分析EMFN 模型的雙參數及單參數分岔、ISI分岔和最大Lyapunov 指數, 揭示其伴有混沌及無混沌的加周期分岔結構、混合模式放電和共存模式放電等現象, 分析了電場和磁場強度影響其放電節律的轉遷規律; 3) 利用Washout 控制器對EMFN模型進行控制, 通過改變其亞臨界Hopf 分岔, 消除不期望的隱藏放電現象.

EMFN 模型中的這些新的動力學性質及其現象背后的數學理論和生理意義值得去進一步的思考和研究.EMFN 模型比神經元模型(1)具有更豐富的放電現象, 如加周期模式、混合模式、共存模式和隱藏模式等放電現象.由此可見, 新模型在外部刺激下可以體現電磁場的感應現象及其反饋功能, 也可以利用電場和磁通變量作為接收外部電磁輻射的橋梁, 進而研究電磁輻射對神經元的影響,這將是本文的后續研究工作.

本文的研究是從數學和物理的角度出發, 將電場和磁場引入神經元模型.通過探討發現新建模型具有更多的分岔參數, 可以檢測出更豐富的電活動轉換模式, 如加周期模式、混合模式、共存模式和隱藏模式等放電現象.由于神經元不同的放電節律承載著不同的刺激信號, 對神經信息編碼具有重要的影響, 因此, 新建模型的放電節律會影響神經信息編碼, 進而影響生物神經系統.同時, 本文從電場和磁場的角度設計新的神經元模型, 豐富了神經元對感知外界刺激的類型和形式, 為完善節律轉遷的理論框架及理解神經編碼的機制提供了一定的依據.另外, 一些神經系統疾病(如帕金森氏癥和癲癇)是由于神經元放電節律異常導致的, 因此,本文的發現有助于探索電場和磁場刺激下神經元的放電節律及其產生機理, 為揭示電磁場對生物神經系統的影響, 以及探尋一些神經性疾病的致病機理提供了思路.