程元芬，劉深泉

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

前言

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

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

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

1 模型描述

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

其中，

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

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

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

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

和系統(tǒng)

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

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

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

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

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

圖1 當gR=260 Am-2V-1，gH=13 Am-2V-1時，產(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”型簇放電

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

圖2 當gR=105 Am-2V-1，gH=40 Am-2V-1時，產(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分岔方向，本研究將計算該Hopf分岔的一階Lyapunov系數(shù)。當gR=105 Am-2V-1，H=1.472 899 Am-2時，重寫該系統(tǒng)的快變子系統(tǒng)：

其中，

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

其中，

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

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

考慮系統(tǒng)

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

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

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

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

其中，

因此，有：

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

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

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

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

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

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

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

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

圖3 當gR=16.5 Am-2V-1，gH=120 Am-2V-1時，產(chǎn)生“fold/fold”點-點滯后環(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)元耦合強度，gL為泄露離子通道電導，Iapp1和Iapp2分別為加入到兩個神經(jīng)元中直流電，Iext1和Iext2為添加到兩個神經(jīng)元中的交流電刺激，Cn表示ML神經(jīng)元模型中的膜電容，gK和gca表示鉀離子和鈣離子的電導，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 耦合強度對耦合神經(jīng)元膜電位的影響

固定Iapp1=0.35 mA/cm2，Iapp2=0.30 mA/cm2，不添加交流電刺激，在此基礎上研究耦合強度變化對MORRIS-LECAR神經(jīng)元膜電位發(fā)放模式的影響。圖4a～d為不同耦合強度所對應的膜電位時間歷程圖。當耦合強度為gc=ms/cm2時，耦合模型中的MORRIS-LECAR膜電位只有一種類型的簇；當gc=0.1 ms/cm2時，膜電位的一個周期內(nèi)包含兩個簇，其中左邊的簇包含12個峰數(shù)量，右邊的簇中包含2個峰數(shù)量；將gc增大到0.2 ms/cm2時，構成一個周期的兩個簇中峰的數(shù)量都為2，但左邊簇的峰值明顯比右邊簇的峰值大；繼續(xù)增大gc到0.3 ms/cm2，結果顯示，膜電位周期性的出現(xiàn)一種簇。顯然，耦合強度的變化對耦合神經(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為不同交流電頻率所對應的膜電位時間歷程圖。當ω=0 rad/ms時，膜電位的一個周期包含兩個不同的簇，左邊簇的峰數(shù)量為3，右邊簇的峰數(shù)量為2；當ω=0.028 rad/ms時，周期內(nèi)兩種不同簇的間距減小，并且兩種簇中所含峰的數(shù)量都明顯增加；當ω增大到0.086 rad/ms時，膜電位時間歷程圖中只出現(xiàn)了一種簇；繼續(xù)增大ω到0.14 rad/ms時，該歷程圖一個周期內(nèi)出現(xiàn)了4種不同的簇。顯然，交流電刺激頻率的改變對膜電位有明顯的影響。

4 小結

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

圖4 當Iapp1=0.35 ms/cm2，Iapp2=0.3 ms/cm2時，不同的耦合強度所對應的膜電壓圖Fig 4 Membrane potential diagrams in different coupling strengths,withIapp1=0.35 ms/cm2,Iapp2=0.3 ms/cm2

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