程元芬，劉深泉

華南理工大學(xué)數(shù)學(xué)學(xué)院，廣東廣州510640

前言

Hodgkin等［1］在1952年通過電生理實(shí)驗(yàn)研究發(fā)現(xiàn)槍烏賊軸突動作電位的產(chǎn)生與4種非線性的、時(shí)間依賴變量的相互作用有關(guān)，這4個變量分別是膜電位V、鈉離子激活變量m、鈉離子失活變量h以及鉀離子激活變量n。以Hodgkin-Huxley模型為基準(zhǔn)，Mccormick等［2］通過實(shí)驗(yàn)研究發(fā)現(xiàn)，影響人類及哺乳動物新皮層神經(jīng)元動作電位產(chǎn)生的因素更加復(fù)雜，至少包括12種不同離子電流的相互作用。除此之外，還包括低閾值的鈣離子電流It、慢變后超極化鉀離子電流IAHP等多種電流的相互作用，使得皮層神經(jīng)元的發(fā)放模式呈現(xiàn)出多種多樣的形態(tài)。這些離子通道相互作用產(chǎn)生的復(fù)雜性使得對皮層神經(jīng)元的分析變得困難，所以人們開始尋求用簡單的數(shù)學(xué)模型來對其進(jìn)行研究。Morris等［3］在1981年提出著名的Morris-Lecar模型，該模型在不同的參數(shù)下能得到不同的神經(jīng)元放電模式；Hindmarsh等［4］在1984年提出一個簡單且易于分析的三維神經(jīng)元模型，但這兩個模型的缺點(diǎn)在于模型并不遵守歐姆定律；后來，Wilson［5］提出一個四維神經(jīng)元模型，它雖然改進(jìn)了Hindmarsh-Rose模型中不遵守歐姆定律的這一缺點(diǎn)，但由于該模型是四維的，這給系統(tǒng)的分析帶來了不便。

不同的放電模式對應(yīng)著神經(jīng)元不同的信息編碼方式，利用動力系統(tǒng)分岔理論和快慢動力學(xué)方法來研究神經(jīng)元模型不同類型的放電活動，被認(rèn)為是行之有效的方法［6-10］，推動了神經(jīng)動力學(xué)的發(fā)展。近年來，利用快慢動力學(xué)的方法，楊卓琴等［11］對Chay神經(jīng)元模型的簇放電模式進(jìn)行了研究；Wang等［12-13］研究了胰腺β細(xì)胞神經(jīng)元模型及呼吸神經(jīng)元模型的簇放電模式；Zhan等［14-15］對Purkinje神經(jīng)元模型及垂體模型的放電模式進(jìn)行了研究。可見快慢動力學(xué)對于時(shí)間尺度差異很大的系統(tǒng)的分析是非常有效的，本研究也將采用該方法來對模型的簇放電模式進(jìn)行研究。

為了保持三維系統(tǒng)易于分析的優(yōu)點(diǎn)，又能遵守生物的生理特性，Zhao等［16］提出一個哺乳動物新皮層神經(jīng)元模型，該模型是將Wilson模型中的快變量與Hindmarsh-Rose模型中的慢反饋?zhàn)兞肯嘟Y(jié)合得到的一個三維模型，所得模型既保持了三維模型易于分析的優(yōu)點(diǎn)，又在一定程度上遵守了生物的生理特性，同時(shí)具有豐富的發(fā)放現(xiàn)象。本研究以該模型為基礎(chǔ)，分析該系統(tǒng)的不同簇放電模式；利用快慢動力學(xué)分析方法，作出快變子系統(tǒng)關(guān)于慢變量的平衡點(diǎn)分岔曲線，分析得出快變子系統(tǒng)的簇發(fā)放類型；將該模型與Morris-Lecar神經(jīng)元模型進(jìn)行電耦合，分析耦合強(qiáng)度及交流電頻率對神經(jīng)元動力學(xué)行為的影響。

1 模型描述

為了保持三維神經(jīng)元模型易于分析的特性，同時(shí)在一定程度上遵守歐姆定律，文獻(xiàn)［16］提出了如下的三維神經(jīng)元模型：

其中，

模型（1）中的第一個方程描述的是膜電壓的變化。其中，V為膜電壓；Cm為膜電容；g(V)為鈉離子激活函數(shù)；gR為鉀離子電導(dǎo)；H表示調(diào)節(jié)電流；gH代表調(diào)節(jié)電流所對應(yīng)的電導(dǎo)；VNa和VK分別為鈉離子和鉀離子的平衡電壓；Iapp為外界刺激電流；R為鉀離子激活變量；R∞為激活變量平衡態(tài)；τR為R松弛時(shí)間變量；Vh表示反轉(zhuǎn)電壓；τH為H松弛時(shí)間變量。詳細(xì)描述見參考文獻(xiàn)［16］。

模型部分參數(shù)取值為：Cm=30 μF/m2，VNa=0.048 V，VK=-0.095 V，V3=-0.038 V，VH=-0.075 04 V，τR=0.56 ms，τH=100 ms，v0=178.1 Am-2V-1，v1=4 758 Am-2V-1，v2=33 800 Am-2V-1，r1=12.9 V-1，r2=330 V-1。

文獻(xiàn)［16］通過數(shù)值模擬說明了該模型在不同的參數(shù)下具有豐富的發(fā)放現(xiàn)象。本研究利用龍格-庫塔算法，結(jié)合MATCONT軟件，得出快變子系統(tǒng)在不同參數(shù)下的幾種分岔類型。同時(shí)理論計(jì)算了Hopf分岔點(diǎn)的一階Lyapunov系數(shù)，從而確定Hopf分岔的方向。最后，將該模型與Morris-Lecar神經(jīng)元模型進(jìn)行電耦合，研究通過耦合得到的新模型的發(fā)放模式。

2 快變子系統(tǒng)分岔曲線上的簇放電類型

在神經(jīng)元模型（1）中，由于H的變化要比R慢很多，因而本研究將H看作慢變量。由此，模型（1）可以看成是由

和系統(tǒng)

構(gòu)成，將慢變量H看作快變子系統(tǒng)［式（3）］的分岔參數(shù)，從而對系統(tǒng)［式（1）］進(jìn)行快慢動力學(xué)分析。

2.1 經(jīng)由“fold/homoclinic”滯后環(huán)的“fold/homoclinic”型簇放電

當(dāng)gR=260 Am-2V-1，gH=13 Am-2V-1時(shí)，神經(jīng)元模型（1）中的膜電位時(shí)間歷程圖如圖1a所示，運(yùn)用快慢動力學(xué)和MATCONT軟件來研究該系統(tǒng)的簇放電模式，得到圖1b。

分析圖1b可知，快變子系統(tǒng)（2）的平衡點(diǎn)曲線為一條Z型曲線，由3部分構(gòu)成。其中穩(wěn)定的焦點(diǎn)和不穩(wěn)定焦點(diǎn)以Hopf分岔點(diǎn)為分界點(diǎn)構(gòu)成了曲線的上支，即隨著分岔參數(shù)H的增大，曲線上的穩(wěn)定焦點(diǎn)經(jīng)由Hopf分岔點(diǎn)變成了不穩(wěn)定焦點(diǎn)，同時(shí)出現(xiàn)穩(wěn)定的極限環(huán)；鞍點(diǎn)構(gòu)成了該曲線的中支；曲線的下支由穩(wěn)定的結(jié)點(diǎn)構(gòu)成。將圖1a中的簇發(fā)放軌線圖疊加在快變子系統(tǒng)平衡點(diǎn)曲線上，分析在該參數(shù)條件下的簇發(fā)放類型。

隨著控制變量H的減小，平衡點(diǎn)分岔曲線下支相應(yīng)于靜息狀態(tài)的穩(wěn)定結(jié)點(diǎn)消失，向上轉(zhuǎn)遷到相應(yīng)于放電狀態(tài)的穩(wěn)定極限環(huán)周圍。經(jīng)由Hopf分岔點(diǎn)分支出的穩(wěn)定極限環(huán)隨著分岔參數(shù)H的增大，逐漸逼近由鞍點(diǎn)組成的中支，最后碰到鞍點(diǎn)變成鞍點(diǎn)同宿軌，經(jīng)由鞍點(diǎn)同宿軌線回到相應(yīng)于靜息狀態(tài)的穩(wěn)定結(jié)點(diǎn)處，放電狀態(tài)結(jié)束。因此，靜息狀態(tài)與發(fā)放狀態(tài)相互轉(zhuǎn)遷的放電類型為“fold/homoclinic”型簇放電。曲線的下狀態(tài)經(jīng)由鞍結(jié)點(diǎn)躍遷至曲線的上狀態(tài)，上狀態(tài)經(jīng)由鞍點(diǎn)同宿軌回到下狀態(tài)。因此滯后環(huán)產(chǎn)生的簇放電為“fold/homoclinic”型簇放電，根據(jù)文獻(xiàn)［10］提出的放電類型分類方法，此簇放電的模式為經(jīng)由“fold/homoclinic”滯后環(huán)的“fold/homoclinic”型簇放電。

圖1 當(dāng)gR=260 Am-2V-1，gH=13 Am-2V-1時(shí)，產(chǎn)生經(jīng)由“fold/homoclinic”滯后環(huán)的“fold/homoclinic”型簇放電Fig 1 Fast-slow dynamics of"fold/homoclinic"bursting via"fold/homoclinic"hysteresis loop whengR=260 Am-2V-1,gH=13 Am-2V-1

2.2 經(jīng)由“fold/homoclinic”滯后環(huán)的“subHopf/homoclinic”型簇放電

當(dāng)gR=105 Am-2V-1，gH=40 Am-2V-1時(shí)，系統(tǒng)（1）［式（1）］有圖2a所示的膜電位時(shí)間歷程圖。將圖2a的簇發(fā)放軌線圖附于快變子系統(tǒng)平衡點(diǎn)分岔曲線上，得到圖2b。根據(jù)數(shù)值計(jì)算結(jié)果，Z型平衡點(diǎn)曲線上支存在一個Hopf分岔點(diǎn)，分岔點(diǎn)左邊部分為穩(wěn)定焦點(diǎn)構(gòu)成，右側(cè)由不穩(wěn)定焦點(diǎn)構(gòu)成。分岔曲線的中支和下支分別由鞍點(diǎn)和穩(wěn)定結(jié)點(diǎn)構(gòu)成。

圖2 當(dāng)gR=105 Am-2V-1，gH=40 Am-2V-1時(shí)，產(chǎn)生經(jīng)由“fold/homoclinic”滯后環(huán)的“subHopf/homoclinic”型簇放電Fig 2 Fast-slow dynamics of"subHopf/homoclinic"via"fold/homoclinic"hysteresis loop whengR=105 Am-2V-1,gH=40 Am-2V-1

為確定該參數(shù)下的Hopf分岔方向，本研究將計(jì)算該Hopf分岔的一階Lyapunov系數(shù)。當(dāng)gR=105 Am-2V-1，H=1.472 899 Am-2時(shí)，重寫該系統(tǒng)的快變子系統(tǒng)：

其中，

快變子系統(tǒng)的Jacobian矩陣為：

其中，

當(dāng)H=1.472 899 Am-2時(shí)，快變子系統(tǒng)的平衡點(diǎn)為(- 0.017 212,0.710 583)，可得該平衡點(diǎn)處的Jacobian矩陣為：

經(jīng)計(jì)算，該矩陣的特征值為λ1,2=±wi,w=3.123 0，由此可知，在該平衡點(diǎn)處產(chǎn)生了Hopf分岔。另外，λ對應(yīng)的特征向量為q=(0.037 565+0.065 698i,1)Τ，另取一個向量p使得p滿足AΤp=-wip，且p,q=1，計(jì)算得到：

考慮系統(tǒng)

其中，F(xiàn)(x)光滑。F關(guān)于x的Taylor展開F(x)=Ο(||x2||)至少從二次方開始，A是平衡點(diǎn)處的Jacobian矩陣，F(xiàn)(x)可寫成如下形式：

這里B(x,y)和C(x,y,z)是x,y,z∈R2的對稱多重線性向量函數(shù)。它們的表達(dá)式如下所示：

其中，ξ=(ξ1,ξ2)為系統(tǒng)的平衡點(diǎn)。

將快變子系統(tǒng)寫成如下形式：

其中，

因此，有：

可得一階Lyapunov系數(shù)為：

因此該Hopf分岔方向是亞臨界的。

根據(jù)圖2b，快變子系統(tǒng)平衡點(diǎn)曲線有兩處靜息態(tài)，分別對應(yīng)于平衡點(diǎn)曲線的下支及曲線上支Hopf分岔點(diǎn)的左側(cè)。隨著控制變量H的減小，下支靜息態(tài)經(jīng)由鞍結(jié)點(diǎn)轉(zhuǎn)遷到曲線上支的靜息狀態(tài)，在上支靜息態(tài)周圍發(fā)生小幅衰減振蕩后收斂于穩(wěn)定焦點(diǎn)。而后隨著H的增大，靜息態(tài)經(jīng)由Hopf分岔消失，轉(zhuǎn)遷為發(fā)放狀態(tài)，經(jīng)過幾個峰值不同的發(fā)放后重新轉(zhuǎn)遷回靜息狀態(tài)。

此外，系統(tǒng)從下狀態(tài)經(jīng)由鞍結(jié)點(diǎn)躍遷到上狀態(tài)，經(jīng)由鞍點(diǎn)同宿軌線從上狀態(tài)躍遷到下狀態(tài)。因此，在該參數(shù)下的簇放電模式為經(jīng)由“fold/homoclinic”滯后環(huán)“subHopf/homoclinic”型簇放電。

2.3 “fold/fold”點(diǎn)-點(diǎn)滯后環(huán)型簇放電

當(dāng)gR=16.5 Am-2V-1，gH=120 Am-2V-1時(shí)，系統(tǒng)（1）產(chǎn)生如圖3a的發(fā)放模式?？熳冏酉到y(tǒng)的平衡點(diǎn)曲線呈現(xiàn)Z型，曲線的上支、中支、下支分別由穩(wěn)定焦點(diǎn)、鞍點(diǎn)、穩(wěn)定結(jié)點(diǎn)構(gòu)成。將圖3a的發(fā)放圖附于平衡點(diǎn)曲線上，得到快變子系統(tǒng)的快慢動力學(xué)分析圖3b。由于在該參數(shù)下，快變子系統(tǒng)分岔曲線不存在Hopf分岔點(diǎn)，因此不存在相應(yīng)于放電狀態(tài)的穩(wěn)定極限環(huán)，故不需要討論靜息狀態(tài)和放電狀態(tài)相互轉(zhuǎn)遷的分岔類型，只需要討論與滯后環(huán)的產(chǎn)生有關(guān)的分岔。從圖中可以看出下狀態(tài)經(jīng)由鞍結(jié)點(diǎn)轉(zhuǎn)遷到上狀態(tài)，上狀態(tài)又經(jīng)由鞍結(jié)點(diǎn)轉(zhuǎn)遷到下狀態(tài)。因此，在該參數(shù)下的系統(tǒng)簇放電模式為“fold/fold”點(diǎn)-點(diǎn)滯后環(huán)型簇放電。

3 耦合神經(jīng)元模型的發(fā)放特性

以該三維模型和Morris-Lecar模型為基礎(chǔ)，將這兩個具有豐富發(fā)放模式的神經(jīng)元模型進(jìn)行電耦合，得到一個新的六維神經(jīng)元模型，其動力學(xué)方程如下：

圖3 當(dāng)gR=16.5 Am-2V-1，gH=120 Am-2V-1時(shí)，產(chǎn)生“fold/fold”點(diǎn)-點(diǎn)滯后環(huán)型簇放電Fig 3 Fast-slow dynamics of“fold/fold”hysteresis loop bursting of point-point type forgR=16.5 Am-2V-1,gH=120 Am-2V-1

其中，V1表示Morris-Lecar模型神經(jīng)元膜電位，V2為三維哺乳動物神經(jīng)元膜電位，gc神經(jīng)元耦合強(qiáng)度，gL為泄露離子通道電導(dǎo)，Iapp1和Iapp2分別為加入到兩個神經(jīng)元中直流電，Iext1和Iext2為添加到兩個神經(jīng)元中的交流電刺激，Cn表示ML神經(jīng)元模型中的膜電容，gK和gca表示鉀離子和鈣離子的電導(dǎo)，Vca為鈣離子平衡電壓，w表示門控變量。

此外：

參數(shù)取值為VL=-0.5 mV，VK=-1.1 mV，VCa=1.0 mV，gL=0.5 ms/cm2，gk=2.0 ms/cm2，gCa=1.2 ms/cm2，v11=-0.01mV，v22=-0.15mV，v33=-0.10mV，v44=-0.05mV。

3.1 耦合強(qiáng)度對耦合神經(jīng)元膜電位的影響

固定Iapp1=0.35 mA/cm2，Iapp2=0.30 mA/cm2，不添加交流電刺激，在此基礎(chǔ)上研究耦合強(qiáng)度變化對MORRIS-LECAR神經(jīng)元膜電位發(fā)放模式的影響。圖4a～d為不同耦合強(qiáng)度所對應(yīng)的膜電位時(shí)間歷程圖。當(dāng)耦合強(qiáng)度為gc=ms/cm2時(shí)，耦合模型中的MORRIS-LECAR膜電位只有一種類型的簇；當(dāng)gc=0.1 ms/cm2時(shí)，膜電位的一個周期內(nèi)包含兩個簇，其中左邊的簇包含12個峰數(shù)量，右邊的簇中包含2個峰數(shù)量；將gc增大到0.2 ms/cm2時(shí)，構(gòu)成一個周期的兩個簇中峰的數(shù)量都為2，但左邊簇的峰值明顯比右邊簇的峰值大；繼續(xù)增大gc到0.3 ms/cm2，結(jié)果顯示，膜電位周期性的出現(xiàn)一種簇。顯然，耦合強(qiáng)度的變化對耦合神經(jīng)元膜電位有明顯的影響。

3.2 交流電刺激對耦合神經(jīng)元膜電位的影響

固定Iapp1=0 mA/cm2，Iapp2=0.3 mA/cm2，Iext2=0 mA/cm2，gc=0.2 ms/cm2，在模型（4）的第一個方程中加入交流電刺激Iext1=0.5 sin(ωt)，通過改變交流電頻率，由此來研究交流電頻率變化對耦合神經(jīng)元中MORRIS-LECAR膜電位發(fā)放模式的影響。

圖5a～d為不同交流電頻率所對應(yīng)的膜電位時(shí)間歷程圖。當(dāng)ω=0 rad/ms時(shí)，膜電位的一個周期包含兩個不同的簇，左邊簇的峰數(shù)量為3，右邊簇的峰數(shù)量為2；當(dāng)ω=0.028 rad/ms時(shí)，周期內(nèi)兩種不同簇的間距減小，并且兩種簇中所含峰的數(shù)量都明顯增加；當(dāng)ω增大到0.086 rad/ms時(shí)，膜電位時(shí)間歷程圖中只出現(xiàn)了一種簇；繼續(xù)增大ω到0.14 rad/ms時(shí)，該歷程圖一個周期內(nèi)出現(xiàn)了4種不同的簇。顯然，交流電刺激頻率的改變對膜電位有明顯的影響。

4 小結(jié)

細(xì)胞膜上的離子通道是影響神經(jīng)元放電活動的重要因素，影響人類哺乳動物新皮層神經(jīng)元膜電壓變化的因素有很多，對其的分析也相當(dāng)復(fù)雜，文獻(xiàn)［16］所提出的簡化模型使得對哺乳動物新皮層神經(jīng)元放電模式的分析變得簡單。本研究采用快慢動力學(xué)分岔分析方法，對gR、gH取不同的值，以研究該三維哺乳動物新皮層神經(jīng)元模型的簇放電模式。對于該系統(tǒng)的快變子系統(tǒng)平衡點(diǎn)隨控制變量H的變化所形成的分岔曲線，主要出現(xiàn)了以下3種簇放電模式：（1）當(dāng)gR=260 Am-2V-1，gH=13 Am-2V-1時(shí)，系統(tǒng)的分岔表現(xiàn)為經(jīng)由“fold/homoclinic”滯后環(huán)的“fold/homoclinic”型簇放電；（2）當(dāng)gR=105 Am-2V-1，gH=40 Am-2V-1時(shí)，系統(tǒng)的分岔類型為經(jīng)由“fold/homoclinic”滯后環(huán)的“subHopf/homoclinic”型簇放電。同時(shí)，通過理論計(jì)算可知，快變子系統(tǒng)分岔曲線上Hopf分岔點(diǎn)對應(yīng)的一階Lyapunov系數(shù)為0.443 025，由此可判斷該Hopf分岔為亞臨界Hopf分岔；（3）當(dāng)gR=16.5 Am-2V-1，gH=120 Am-2V-1時(shí)，系統(tǒng)的簇放電模式“fold/fold”點(diǎn)-點(diǎn)為滯后環(huán)型簇放電。將該模型與MORRIS-LECAR神經(jīng)元模型進(jìn)行電耦合，得到一個新的六維神經(jīng)元模型。通過改變耦合強(qiáng)度和交流電頻率，膜電位出現(xiàn)豐富的變化。結(jié)果表明，耦合強(qiáng)度和交流電頻率的變化對神經(jīng)元膜電位產(chǎn)生非常明顯的影響。

圖4 當(dāng)Iapp1=0.35 ms/cm2，Iapp2=0.3 ms/cm2時(shí)，不同的耦合強(qiáng)度所對應(yīng)的膜電壓圖Fig 4 Membrane potential diagrams in different coupling strengths,withIapp1=0.35 ms/cm2,Iapp2=0.3 ms/cm2

圖5 當(dāng)Iext1=0.5sin(ωt)時(shí)，不同的交流電頻率所對應(yīng)的膜電位圖Fig 5 Membrane potential diagrams at different AC frequencies,withIext1=0.5sin(ωt)