計算相位響應曲線的方波擾動直接算法?

謝勇程建慧

(西安交通大學航天航空學院,機械結構強度與振動國家重點實驗室,西安 710049)

計算相位響應曲線的方波擾動直接算法?

謝勇?程建慧

(西安交通大學航天航空學院,機械結構強度與振動國家重點實驗室,西安 710049)

(2016年11月25日收到;2017年1月6日收到修改稿)

通過相位響應曲線可對具有極限環周期運動的動力系統的性質有更為深入的理解.神經元是一個典型的動力系統,因此相位響應曲線提供了一種研究神經元重復周期放電行為的新思路.本文提出一種求解相位響應曲線的方法,即方波擾動的直接算法,通過Hodgkin-Huxley,FitzHugh-Nagumo,Morris-Lecar和Hindmarsh-Rose神經元模型驗證該算法可計算周期峰放電、周期簇放電的相位響應曲線.該算法克服了其他算法在運用過程中的局限性.利用該算法計算結果表明:周期峰放電的相位響應曲線類型是由其分岔類型所決定;在Morris-Lecar模型中發現一種開始于Hopf分岔終止于鞍點同宿軌道分岔的閾上周期振蕩,其相位響應曲線屬于第二類型.通過大量的相位響應曲線的計算發現相位響應的相對大小及正負性僅取決于擾動所施加的時間,而且周期簇放電的相位響應曲線比周期峰放電的相位響應曲線更為復雜.

相位響應曲線,峰放電,簇放電,分岔

1 引 言

神經元和神經元所組成的神經系統具有高度復雜的動力學行為,呈現明顯的非線性現象,因此從非線性動力學角度研究和理解神經元和神經系統的行為規律是十分必要的.其研究的一般思路是建立動力學模型,然后利用動力學系統理論解讀神經元和神經系統活動規律及其動力學機理.在神經元建模方面,英國生理學家Hodgkin和Huxley［1]建立了Hodgkin-Huxley(HH)模型,它能很好地描述槍烏賊的神經元峰放電(spiking)行為.這一開創性工作成為定量研究可興奮生物細胞電生理特性的里程碑.為了便于理論分析,在1961年和1962年,FitzHugh和Nagumo等分別獨立地導出了一個由多項式表達的二維FitzHugh-Nagumo(FHN)模型［2,3].該模型盡管形式非常簡潔,卻抓住了神經元興奮的本質特征.1981年,Morris和Lecar［4]結合HH模型和FHN模型各自的優點提出了一個描述藤壺(一種甲殼綱的小動物)肌肉纖維電生理特征的二維Morris-Lecar(ML)模型,現在也被作為神經元模型而廣泛引用.隨后,Hindmarsh和Rose［5]為了解釋“尾電流反轉(tail current reversal)”現象提出了Hindmarsh-Rose(HR)模型.HR模型可以呈現峰放電和簇放電(bursting)行為.目前針對不同生物和不同類型的神經元已經建立了大量的模型［6].有了神經元模型,接下來就可以通過非線性動力學理論分析神經元復雜而豐富的放電行為,并揭示其發生的動力學機理.已有大量的理論或數值分析確定神經元模型中分岔、混沌和分形等非線性現象［7?17].噪聲無處不在,對隨機共振現象的研究表明噪聲可能增強神經元、神經元網絡以及神經系統檢測微弱信號的能力［18?22].

通過對神經元模型的數值計算很容易再現神經元各種復雜的放電模式,比如周期峰放電、周期簇放電、混沌峰放電、混沌簇放電和整數倍放電等.本文主要關注周期峰放電和周期簇放電兩種放電模式.這兩種放電模式從非線性動力學角度講都是極限環振子,而描述極限環的最簡單的方式就是采用振子的相位.通過相位變換就可以把復雜的狀態空間模型映射到一維的相位模型,這樣有助于獲得振子系統的解析解.處于極限環運動狀態的振子系統對外界刺激的響應特征可以通過相位響應曲線(phase response curve,PRC)進行刻畫［23,24],相位響應曲線有時也稱為相位重置曲線(phase resetting curve,PRC).PRC是一個關于相位的函數,它是指極限環振子在不同相位處由外界擾動引起振蕩周期的暫態變化.

PRC是由文獻[25,26]在研究生物節律的重置過程中引入的,在很多方面都有極其重要的作用,隨后被應用到神經元系統中.在相空間中PRC表示擾動后周期振子相位的提前或滯后.如果周期振子的周期縮短了,則稱為相位提前(phase advance);相反,如果周期變長了,則稱為相位滯后(phase delay).在研究生物節律時,PRC已成為研究的標準工具.比如人體對光刺激的PRC是調整睡眠時間療法的理論基礎.在心臟病學中,PRC可用來理解各種心律失常的發生［27,28].在神經科學中,相位簡化模型可用來簡化復雜的神經科學的生物模型［29?31].如果知道一個振子系統的PRC,就可以很容易地預測該振子在受到外界刺激時的行為.另外,PRC可用來預測周期刺激的拖帶(entrainment)以及耦合振子的相鎖(phase locking),還可以用來理解相干振蕩器、行波(traveling waves)等動力學行為.耦合神經元的相互作用函數可以通過PRC數值計算得到［32].文獻[33]詳細介紹了神經科學中PRC及其應用.已有的研究表明周期峰放電的神經元的PRC可分為兩大類［23]:第一類型PRC是指PRC取值全為正或者全為負,這樣興奮性擾動將導致峰放電提前或者滯后;而第二類型PRC是指PRC取值是兩相的,部分為正,部分為負,這樣將導致當擾動施加在不同的相位處時,部分相位處將引起放電提前,而部分相位處引起放電滯后.

而PRC的類型與神經元從靜息到周期峰放電所跨越的分岔類型密切相關.具體地講,具有第一類型PRC的神經元跨越了不變圓上的鞍結分岔(saddle-node on invariant circle bifurcation,SNIC),而具有第二類型PRC的神經元則經歷了Hopf分岔.因此PRC的類型與神經元的興奮性類型一一對應.

周期簇放電的PRC與周期峰放電的PRC定義類似,是指簇放電神經元在受到外界擾動時下一個簇到來時刻與未擾動的簇之間的相位差.對于周期簇放電的PRC與其分岔類型之間的關系,現在的研究還比較少,它們之間的關系還不清楚.不過,要弄清這一關系的前提就是首先要計算獲得周期簇放電的PRC.

PRC可通過實驗或理論分析得到.在進行理論分析時,需要知道動力系統具體的表達式.對于一些簡單的模型,其PRC可通過相位簡化方法得到其解析解;但對復雜的、非線性的神經元模型,往往無法解析求解,因此有必要通過數值計算得到其PRC.目前計算PRC的基本思路有兩種［32]:一種是根據PRC的定義,即用神經元振子對短的脈沖刺激的響應直接計算,該思路比較簡單,但擾動是有限幅值擾動,導致計算結果不是很精確;另一種是將神經元振子在極限環附近線性化,求解線性伴隨方程,然而求解伴隨方程需要向后積分且不穩定,因此該算法并不簡單.不過該算法是一些軟件求解PRC的依據,比如XPPAUT［34]和MatCont［35].對于前一種思路所對應的直接算法中,通常在膜電位上施加典型的峰擾動,這一典型峰是取自簇放電本身的一個尖峰,該擾動形式更接近實際的生理過程.但該方法是針對簇放電的,而且典型的峰擾動所需確定的參數較多,比如對于HR模型,需要確定5個參數,方程增加一維,計算費時［36].無限小擾動的改進算法也是直接算法的一種［37],該算法計算速度快.但該算法本質在于求解特征矩陣,在計算簇放電神經元模型的PRC時,容易溢出,出現錯誤.例如該算法應用于Chay模型和Wang模型時,就出現了溢出錯誤.因此基于上述算法的局限性,本文提出了一種基于PRC定義的直接算法,以方波擾動作為擾動形式,該算法可適用于峰放電和簇放電.

本文內容安排如下:第2部分首先介紹我們的直接算法,通過HH模型驗證該算法的可行性,然后將該算法運用到FHN模型中,揭示FHN神經元的PRC類型;第3部分利用本文的直接算法通過ML模型探究周期峰放電模型的PRC的表現特征,揭示分岔機理與PRC類型之間的關系;第4部分通過HR模型驗證該算法同樣適用于計算周期簇放電的PRC;第5部分是結論.

2 PRC的直接算法及算例

2.1 方波擾動的直接算法

對一個自治動力系統,其模型可定義為

其中,X=(x1,x2,...,xN)為N 維相空間的狀態變量.記該系統的解為X(t)=Ψ(X0,t),其中X0=X(0)是初始值.假定該系統存在一個穩定的極限環,其周期為T.在極限環上選取一點作為初始點,即,則可得到一個沿著極限環移動的周期為T的跡線,即.而位于極限環上的一點的相位可以由兩種不同的方式進行定義:一種是空間相位,按照極限環的長度均勻增加進行度量;另一種以一狀態變量沿著極限環時間均勻增加進行度量.根據相位簡化方法,本文采用后一種定義,因為該定義可得到一個簡單的相位公式.如果認為極限環上的初始點相位為0,則極限環上其他點的相位為

對系統在極限環上一點x0處施加一擾動?F,若該擾動僅使得該點沿極限環偏移,則該擾動的作用僅是增加或減少動力系統的周期,或者說延遲或提前擾動后的跡線回到點x0處的時間,則該點的相位響應為

其中t0為未加擾動從點x0處出發第一次回到點x0處的時間,0為擾動后第一次回到x0處的時間,為擾動后的周期.PR表示相位響應(phase response).

小幅值的擾動?F通常用來模擬一個瞬時的增加或減小后突觸神經元的電流作用,而忽略關于離子通道參數的影響.在計算此類擾動的PRC時,先數值積分未擾動的系統F到施加擾動相位(或時刻),再積分擾動后的系統,積分區間長度為擾動持續的時間長度,然后繼續數值積分原系統F,最后相位響應根據(3)式求解,則PRC由施加在[0,1)區間上的擾動所對應的相位響應構成.這就是本文計算PRC的直接算法.在計算PRC過程時,常常假定擾動后的跡線會在一次放電環后回到極限環周期軌道上,并以下一次放電環上的參考相位(周期峰放電一般選取膜電位最高的尖峰時刻)來測量相位響應.此外還需要注意的是,選取未擾動的極限環軌道上的一點作為參考點表明開始時間,記為t=0,相位θ=0.本文中的擾動為方波擾動,其擾動幅值為Ip,持續時間為tp,即在擾動持續時間tp內擾動幅值保持常數Ip,因此稱為方波擾動.該擾動施加在表示神經元膜電位變化的微分方程右端項上,對有生理意義的神經元模型而言,可以表示為,其中C為膜電容.記膜電位在受到外界擾動時的PRC為PRC1,由于本文后面的計算都是在膜電位上施加外界擾動,因此我們用PRC來表示PRC1,不再進行區分.關于方波擾動的參數Ip,tp的選擇,對周期峰放電而言,Ip的大小取決于膜電位的變化幅值以及模型中外加電流的大小(若存在),可取膜電位幅值或者外加電流大小的1/20—1/5;tp的值取決于周期的大小,并且盡量持續時間較短,一般情況下選取不超過峰放電的脈沖寬度.對周期簇放電而言,tp的值一般不超過簇內一個尖峰動作電位的持續時間;Ip的取值類似周期峰放電,但對于周期簇放電,要注意在擾動作用下簇內的尖峰個數不能發生變化.具體的算法實現過程為:

1)計算神經元周期峰放電或周期簇放電所對應的極限環及其周期,記周期為T0,對神經元系統進行時間尺度變換,將周期變換為1,可參考文獻[37];

2)在極限環上選取參考相位,記為t0(對周期峰放電,一般選取膜電位最大的時刻作為參考相位;對周期簇放電,一般選擇第一個尖峰膜電位最大的時刻作為參考相位);

3)將極限環按時間等間距分為n個區間,n=1/tp,在不存在擾動時,按未擾動系統數值積分,在存在擾動的[(m?1)tp,mtp]時間范圍內,按擾動系統數值積分,其中1≤m≤n;記錄第一次回到參考相位的時間,記為0;根據(3)式計算每個節點處的PRC值,若大于0則為相位提前,若小于0表示相位滯后.

2.2 HH神經元模型周期峰放電的PRC

為了與現有的文獻計算結果比較來說明本文提出的直接算法的正確性,特意選擇文獻[37]中的HH神經元模型算例.HH神經元模型的具體形式如下:

其中V為膜電位;C為膜電容;m,h為Na離子通道的門控變量;n為K離子通道的門控變量;gK,gNa和gL分別為K離子通道、Na離子通道和漏電流的最大電導;VK,VNa和VL分別為它們各自的反轉電位;Iapp為外加電流強度.an(V),bn(V),am(V),bm(V),ah(V)和bh(V)滿足:

選取參數

其中電導單位為mS/cm2,電壓單位為mV,電容單位為μF/cm2,電流單位為μA/cm2.

在上述參數取值條件下HH模型存在兩個穩定的極限環,分別稱為外極限環和內極限環.本文所取初始條件分別為(5.25,0.07,0.60,0.011)和(7.25,0.17,0.13,0.55),如圖1所示.

通過分岔軟件XPPAUT對HH模型做膜電位關于外加電流強度的分岔分析,發現有兩個區間存在兩個穩定極限環共存的情形,它們分別是[34.95,38.51]和[40.01,41.28],如圖2所示.圖2(b)是圖2(a)的局部放大.圖中粗實線表示穩定平衡態,點劃線表示不穩定平衡態,細實線表示穩定極限環的振蕩最大和最小值,而點線表示不穩定極限環的振蕩最大和最小值(本文中所有分岔圖都采用這種表達方式).HB表示Hopf分岔;FLC表示極限環折疊分岔(fold limit cycle bifurcation),有時也叫極限環鞍結分岔(saddle-node bifurcation of limit cycles).可以看出圖1中兩個共存的極限環正是在第二狹窄區間里取Iapp=41獲得的.

圖1 HH神經元模型兩個共存的穩定極限環在(V,m)相平面上的投影Fig.1.The(V,m)projections of two coexisting stable limit cycles of the HH neuron model.

接下來運用本文的直接算法計算HH模型內、外極限環的PRC.對內極限環所施加的擾動為Ip=10,tp=0.004,計算步長為0.001(周期歸一化為1,以下同此做法);對外極限環所施加的擾動為Ip=5,tp=0.004,計算步長為0.0001.這樣HH神經元模型的內、外極限環的PRC如圖3所示.我們的計算結果與文獻[30]的結果除了幅值不同以外幾乎完全一致,事實上,PRC相對值才具有意義.值得注意的是,由于在極限環上所取的初始參考點不一樣造成PRC開始的位置不同,但波動順序是完全一致的,從而說明本文的算法是可信的.

由圖3可知,無論內極限環還是外極限環,HH神經元模型的PRC值既有正值又有負值,這表明HH神經元模型的PRC為第二類型PRC.由前面的分岔分析可知HH神經元模型從靜息到放電的確經歷了一個次臨界Hopf分岔.

圖2 HH神經元模型膜電位對外加電流強度的分岔圖Fig.2.Bifurcation diagram of the membrane potential V versus the strength of the applied current Iappin the HH neuron model.

圖3 (a)內極限環的PRC;(b)外極限環的PRCFig.3.(a)The PRC of the inner limit cycle;(b)the PRC of the outer limit cycle.

2.3 FHN神經元模型周期峰放電的PRC

FHN神經元模型的放電行為和分岔機理都有相關的研究,但是對其PRC的計算還較為少見,因此本文采用方波擾動的直接算法計算其PRC.FHN神經元模型的形式如下［2,3]:

其中V表示膜電位;w表示恢復變量;a,b,c為系統參數;Iapp為外加電流強度.FHN模型是一個無量綱的神經元模型.

取參數a=0.139,b=2.54,c=0.008,分析FHN神經元模型隨外加電流強度的分岔機理,如圖4所示.

由圖4可知,當I=0.03469時,神經元經歷一個次臨界Hopf分岔由靜息態變為周期峰放電;當I=0.1509時,經歷另外一個次臨界Hopf點,由周期峰放電變為靜息態.運用直接算法計算該模型的PRC,所施加擾動Ip=0.002,tp=0.004,計算步長為0.001,其放電模式與PRC的對應關系如圖5所示.

圖4 FHN神經元模型的分岔圖Fig.4.Bifurcation diagram of the membrane potential versus the applied current strength in the FHN neuron model.

圖5 FHN神經元模型的放電形式及其對應的PRCFig.5.Periodic spiking and its corresponding PRC in the FHN neuron model.

從圖5可以看到FHN神經元模型的PRC取值有正有負,表現為第二類型的PRC.這與神經元所經歷的Hopf分岔密切相關.而且PRC的幅值最大值并不與神經元放電的尖峰時刻存在對應關系,而是取值取決于擾動所施加的時間.同時還可以看到PRC在前半周期取值幾乎接近于0,在后半周期才出現明顯的變化.這是由于在放電尖峰附近對外部擾動不敏感,而隨之出現的超極化區域所占時間又較長,因此導致了較長區域的PRC取值幾乎為0.

3 ML模型的PRC

由于ML模型在不同參數條件下隨外加電流強度有著不同的分岔行為,因此計算ML模型的PRC能更充分地揭示分岔與PRC形狀的關系.ML模型有三種電流,即鉀電流、鈣電流和漏電流,具體形式為［4]

其中V表示膜電位變量;n為恢復變量;C為膜電容;V1,V2,V3和V4為擬合電壓膜片鉗數據的參數;?為快慢時間尺度比率;gCa,gK和gL分別為Ca離子通道、K離子通道和漏電流的最大電導;VCa,VK和VL分別為它們各自的反轉電位;Iapp為外加電流強度.其中電導率單位為mS/cm2,電壓單位為mV,電容單位為μF/cm2,電流單位為μA/cm2.

盡管ML模型是一個簡單的、只含有兩個變量的動力系統,然而在不同的參數集下卻可以表現出豐富的放電行為［35,38].本文選取如下三組參數研究ML模型的動力學特性,如表1所列.

表1 ML模型的參數集Table 1.Parameter sets of the ML model.

在上述三組參數條件下計算了ML模型隨外加電流強度變化的分岔行為.在第一組參數條件下,ML模型從靜息態到周期峰放電在Iapp=38.76處經歷了一個SNIC,呈現第一類型興奮性,如圖6(a)所示.圖6(b)是在第二組參數條件下,ML從靜息態到周期峰放電在Iapp=45.47處經歷了一個次臨界Hopf分岔,呈現第二類型興奮性.特別地,在第三組參數條件下,隨著外加電流強度逐漸增大,ML模型在Iapp=39.96處靜息態通過一個平衡點的鞍結分岔消失,跳到一個由次臨界Hopf分岔而來的極限環上,呈現一個閾上的周期振蕩.這個極限環隨著外加電流強度逐漸減小,在Iapp=35.01處終止一個鞍點同宿軌道分岔(saddle homoclinic orbit bifurcation),如圖6(c)所示.從圖6(c)還可以看到存在一個由兩條豎直的短劃線所夾的狹窄區間,在此區間內,ML模型有兩個穩定的平衡態和一個穩定的極限環,呈現三穩態.左邊的短劃線位于次臨界Hopf分岔處,而右邊的短劃線位于平衡點鞍結分岔處.

針對這三組參數,取外加電流強度Iapp分別為39,45和39,可以發現前面兩種情形ML模型呈現周期峰放電,如圖7(a)和圖7(b)所示;而最后一種情形則表現為閾上周期振蕩,如圖7(c)所示.

圖6 ML模型在不同參數條件的分岔圖 (a)第一組參數;(b)第二組參數;(c)第三組參數Fig.6.Bifurcation diagrams of the membrane potential versus the applied current strength in the ML neuron model with di ff erent parameter sets:(a)The fi rst parameter set;(b)the second parameter set;(c)the third parameter set.

圖7 ML模型在不同參數條件的放電模式 (a)在第一組參數下當Iapp=39時ML模型的周期峰放電;(b)在第二組參數條件下,當Iapp=45時ML模型的周期峰放電;(c)在第三組參數條件下,當Iapp=39時ML模型的閾上周期振蕩Fig.7.Firing patterns of the ML neuron model with di ff erent parameter sets:(a)Periodic spiking with the fi rst parameter set when Iapp=39;(b)periodic spiking with the second parameter set when Iapp=45;(c)suprathreshold periodic oscillation with the third parameter set when Iapp=39.

運用直接算法計算ML模型在前面三組參數條件下當Iapp分別等于39,45和39時的PRC.在第一組參數條件下,當Iapp=39時所施加擾動為Ip=20,tp=0.004,計算步長為0.001,PRC如圖8(a)所示.在第二組參數條件下,當Iapp=45時所施加擾動為Ip=30,tp=0.004,計算步長為0.001,PRC如圖8(b)所示.圖8(c)所示的PRC為在第三組參數條件下當Iapp=39時的結果,此時所施加擾動為Ip=?10,tp=0.005,計算步長為0.001.

圖8 (a)在第一組參數下當Iapp=39時ML模型的周期峰放電及其對應的PRC;(b)在第二組參數條件下當Iapp=45時ML模型的周期峰放電及其對應的PRC;(c)在第三組參數條件下當Iapp=39時ML模型的閾上周期振蕩及其對應的PRCFig.8.Periodic oscillations and their corresponding PRCs in the ML neuron model with di ff erent parameter sets:(a)Periodic spiking and its corresponding PRC with the fi rst parameter set when Iapp=39;(b)periodic spiking and its corresponding PRC with the second parameter set when Iapp=45;(c)suprathreshold periodic oscillation and its corresponding PRC with the third parameter set when Iapp=39.

由ML模型的PRC計算結果可知,如果從靜息態到周期峰放電經歷SNLC分岔,則屬于第一類型興奮性,其PRC的類型確為第一類型,PRC取值幾乎全部為正.當從靜息態到周期峰放電經歷Hopf分岔時,其興奮性屬于第二類型,而其PRC取值既有正又有負.在第三組參數條件下,當Iapp=39時ML模型的PRC既有正值又有負值,表現為第二類型的PRC.實際上,在第三組參數條件下,隨著外加電流強度減小,ML模型經歷一個次臨界Hopf分岔產生閾上周期振蕩,該周期振蕩的極限環終止于一個鞍點同宿軌道分岔.對此情形尚無相關的文獻給出周期振蕩所對應的PRC類型.但由其分岔機理可知,周期振蕩是由Hopf分岔引起的,因此其PRC表現為第二類型也是合乎情理的.同FHN神經元模型的PRC計算結果一樣,ML模型的膜電位變化的尖峰時刻與其PRC取值的局部極值并不對應,這表明PRC僅取決于施加擾動的時間.對ML模型而言,其在不同參數條件下表現出不同的分岔行為,因此決定其PRC可以表現為第一類型或第二類型,從而可通過調節參數來實現PRC類型的改變.

4 方波擾動直接算法計算周期簇放電PRC

為了驗證本文直接算法對周期簇放電PRC計算的有效性,我們選擇HR神經元模型進行PRC計算,其主要原因是目前已有的文獻中只有HR模型的PRC結果可資以參考.HR神經元模型形式如下［5]:

HR模型同FHN模型一樣,是無量綱的.式中x表示神經元膜電位;y是峰變量,代表Na和K等快速離子通道電流;而z是簇變量,代表其他慢離子通道電流;a,b,c,d,s和x0為系統參數;r為快慢時間尺度因子;Iapp為外加電流強度.本文選取a=1.0,b=3.0,c=1.0,d=5.0,s=4.0,x=?1.6,r=0.001.

與前面周期峰放電情形不同的是,在計算周期簇放電PRC時,要注意所施加方波擾動的強度和持續的時間.本文中僅考慮弱擾動,其具體表現為不改變簇放電簇內的尖峰個數.如果尖峰個數在擾動下發生變化,則稱該擾動為強擾動,如圖9所示,實線是未擾動時簇放電模式,虛線是施加擾動后的簇放電模式,顯然簇內尖峰個數發生了變化.當然在實際中強擾動亦有其意義,但本文中不討論這一情形.值得注意的是,在算法實現過程中,參考相位一般選取第一個尖峰所對應的時刻.

當Iapp=1.3時,HR神經元模型簇內尖峰個數為5,其放電模式和對應的極限環如圖10所示.

運用本文直接算法計算在該條件下HR神經元模型周期簇放電的PRC.其中Ip=0.05,tp=10,步長為0.0001,結果如圖11所示.

圖9 強擾動下周期簇放電簇內尖峰個數變化Fig.9.The variation of the number of spiking in a bursting under a strong perturbation.

圖10 HR神經元模型的周期簇放電模式及其對應的極限環Fig.10.Periodic bursting and its corresponding limit cycle in the HR neuron model.

運用無限小擾動的改進算法的計算結果如圖12所示.

對比方波直接算法的結果與無限小擾動的改進算法的結果可知,二者僅存在PRC幅值大小的不同,而形狀、波動趨勢和正負性幾乎完全相同.這表明本文方波擾動直接算法亦可適用于計算周期簇放電的PRC.同時通過HR神經元模型可以發現周期簇放電的PRC比周期峰放電的PRC有更為復雜的表現形式,最為明顯的一點是周期簇放電的PRC存在一個劇烈變化的區域,正好對應于簇內放電部分.

圖11 方波擾動直接算法計算HR神經元模型周期簇放電的PRCFig.11.PRC of periodic bursting in the HR neuron model using the direct method with square wave perturbation.

圖12 利用無限小擾動的改進算法時HR模型的PRCFig.12.PRC of periodic bursting in the HR neuron model using the modi fi ed direct method with in fi nitesimal perturbation.

5 結 論

本文提出了一個基于PRC的定義的方波擾動直接算法,通過HH神經元模型驗證了該算法在計算周期峰放電的PRC是有效的.將該算法運用于FHN神經元模型,揭示了其PRC屬于第二類型,這與FHN神經元模型隨外加電流強度變化從靜息態到周期峰放電所跨越的Hopf分岔動力學機理完全符合.然后考察了在三組不同的參數集條件下ML模型隨外加電流強度變化的分岔行為,并計算了在各自周期峰放電或閾上周期振蕩情形的PRC,表明通過調節相應的參數可改變神經元興奮性類型,從而改變其PRC的類型,揭示了分岔機理與PRC類型之間的對應關系.最后通過HR神經元模型驗證了本文算法對于計算周期簇放電的PRC也是可行的.通過大量的計算表明該算法對不同類型的簇放電具有通用性.我們已經計算了多種周期簇放電的PRC,至于周期簇放電神經元的PRC與周期簇放電的分岔機理之間的關系將在今后的工作進一步展開.本文的方波擾動直接算法計算速度介于無限小擾動的改進算法和典型峰擾動算法之間;在數值積分時積分步長對計算結果有明顯的影響,當積分步長較大時,PRC不準確;當積分步長較小時,則會增加計算時間.方波擾動直接算法是從PRC定義出發設計的算法,計算過程中只需要確定方波擾動幅值Ip和持續時間tp.而典型峰擾動直接算法在計算時需要人為確定5個參數取值,并且方程還需增加一維.其他算法不是涉及伴隨方程求解就是特征方程求解,這一關鍵步驟在計算分岔機理比較復雜的周期簇放電的PRC時往往會求解失敗.這正是我們提出方波擾動直接算法的原因.

[1]Hodgkin A L,Huxley A F 1952J.Physiol.117 500

[2]FitzHugh R 1961Biophys.J.1 445

[3]Nagumo J,Arimoto S,Yoshizawa S 1962Proc.IRE50 2061

[4]Morris C,Lecar H 1981Biophys.J.35 193

[5]Hindmarsh J L,Rose R M 1984Proc.R.Soc.London Ser.B221 87

[6]Xu L F,Li C D,Chen L 2016Acta Phys.Sin.65 240701(in Chinese)[徐泠風,李傳東,陳玲 2016物理學報 65 240701]

[7]Holden A V,Fan Y S 1992Chaos Soliton.Fract.2 221

[8]Holden A V,Fan Y S 1992Chaos Soliton.Fract.2 349

[9]Holden A V,Fan Y S 1992Chaos Soliton.Fract.2 583

[10]Fan Y S,Holden A V 1993Chaos Soliton.Fract.3 439

[11]Izhikevich E M 2000Int.J.Bifurcat.Chaos10 1171

[12]Gong P L,Xu J X 2001Phys.Rev.E63 031906

[13]Ding X L,Li Y Y 2016Acta Phys.Sin.65 210502(in Chinese)[丁學利,李玉葉 2016物理學報 65 210502]

[14]Gu H G,Zhu Z,Jia B 2011Acta Phys.Sin.60 100505(in Chinese)[古華光,朱洲,賈冰 2011物理學報 60 100505]

[15]Jin Q T,Wang J,Wei X L,Deng B,Che Y Q 2011Acta Phys.Sin.60 098701(in Chinese)[金淇濤,王江,魏熙樂,鄧斌,車艷秋2011物理學報60 098701]

[16]Wang H X,Wang Q Y,Lu Q S 2011Chaos Soliton.Fract.44 667

[17]Yang Z Q,Guan T T,Gan C B,Zhang J Y 2011Acta Phys.Sin.60 110202(in Chinese)[楊卓琴,管亭亭,甘春標,張矯瑛2011物理學報60 110202]

[18]Longtin A 1993J.Stat.Phys.70 309

[19]Braun H A,Wissing H,Sch?fer K,Hirsch M C 1994Nature367 270

[20]Wiesenfeld K,Moss F 1995Nature373 33

[21]Yu Y,Wang W,Wang J,Liu F 2001Phys.Rev.E63 021907

[22]Liu F,Wang J,Wang W 1999Phys.Rev.E59 3453

[23]Ermentrout B 1996Neural Comput.8 979

[24]Gutkin B S,Ermentrout G B,Reyes A D 2005J.Neurophysiol.94 1623

[25]Hastings J W,Sweeney B M 1958Biol.Bull.115 440

[26]Johnson C H 1999Chronobiol.Int.16 711

[27]Ikeda N 1982Biol.Cybern.43 157

[28]Tsalikakis D G,Zhang H G,Fotiadis D I,Kremmydas G P,Michalis L K 2007Comput.Biol.Med.37 8

[29]Ermentrout G B,Kopell N 1991J.Math.Biol.29 195

[30]Ermentrout G B 1992SIAM J.Appl.Math.52 1665

[31]Stiger T,Danzl P,Moehlis J,Neto ffT I 2010J.Med.Devices4 027533

[32]Shi X,Zhang J D 2016Chin.Phys.B25 060502

[33]Schultheiss N W,Prinz A A,Butera R J 2012Phase Response Curves in Neuroscience(New York:Springer)p3

[34]Ermentrout G B 2002Simulating,Analyzing,and Animating Dynamical Systems:a Guide to XPPAUT for Researchers and Students(Philadelphia:SIAM)p226

[35]Govaerts W,Sautois B 2006Neural Comput.18 817

[36]Sherwood W E,Guckenheimer J 2010SIAM J.Appl.Dyn.Syst.9 659

[37]Novicenko V,Pyragas K 2011Nonlinear Dynam.67 517

[38]Ermentrout G B,Terman D H 2010Mathematical Foundations of Neuroscience(New York:Springer)p51

PACS:05.45.–a,87.19.ll,87.19.lnDOI:10.7498/aps.66.090501

A direct algorithm with square wave perturbation for calculating phase response curve?

Xie Yong?Cheng Jian-Hui

(State Key Laboratory for Strength and Vibration of Mechanical Structures,School of Aerospace,Xi’an Jiaotong University,
Xi’an 710049,China)

25 November 2016;revised manuscript

6 January 2017)

Neuron is a typical dynamic system,therefore,it is quite natural to study the fi ring behaviors of neurons by using the dynamical system theory.Two kinds of fi ring patterns,i.e.,the periodic spiking and the periodic bursting,are the limit cycle oscillators from the point of view of nonlinear dynamics.The simplest way to describe the limit cycle is to use the phase of the oscillator.A complex state space model can be mapped into a one-dimensional phase model by phase transformation,which is helpful for obtaining the analytical solution of the oscillator system.The response characteristics of the oscillator system in the motion state of the limit cycle to the external stimuli can be characterized by the phase response curve.A phase response curve illustrates the transient change in the cycle period of an oscillation induced by a perturbation as a function of the phase at which it is received.Now it is widely believed that the phase response curve provides a new way to study the behavior of the neuron.Existing studies have shown that the phase response curve of the periodic spiking can be divided into two types,which are closely related to the bifurcation mechanism of neurons from rest to repetitive fi ring.However,there are few studies on the relationship between the phase response curve and the bifurcation type of the periodic bursting.Clearly,the fi rst prerequisite to understand this relationship is to calculate the phase response curve of the periodic bursting.The existing algorithms for computing the phase response curve are often unsuccessful in the periodic bursting.In this paper,we present a method of calculating the phase response curve,namely the direct algorithm with square wave perturbation.The phase response curves of periodic spiking and periodic bursting can be obtained by making use of the direct algorithm,which is veri fi ed in the four neuron models of the Hodgkin-Huxley,FitzHugh-Nagumo,Morris-Lecar and Hindmarsh-Rose.This algorithm overcomes the limitations to other algorithms in the application.The calculation results show that the phase response curve of the periodic spiking is determined by the bifurcation type.We fi nd a suprathreshold periodic oscillation starting from a Hopf bifurcation and terminating at a saddle homoclinic orbit bifurcation as a function of the applied current strength in the Morris-Lecar model,and its phase response curve belongs to Type II.A large amount of calculation indicates that the relative size of the phase response and its positive or negative value depend only on the time of imposing perturbation,and the phase response curve of periodic bursting is more complicated than that of periodic spiking.

phase response curve,spiking,bursting,bifurcation

10.7498/aps.66.090501

?國家自然科學基金(批準號:11272241,11672219)資助的課題.

?通信作者.E-mail:yxie@mail.xjtu.edu.cn

*Project supported by the National Natural Science Foundation of China(Grant Nos.11272241,11672219).

?Corresponding author.E-mail:yxie@mail.xjtu.edu.cn