李揚(yáng) 劉先斌

(南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國(guó)家重點(diǎn)實(shí)驗(yàn)室,南京 210016)

引言

在生物物理學(xué)中,有越來(lái)越多的現(xiàn)象源于分段確定性動(dòng)力系統(tǒng)和連續(xù)時(shí)間馬氏過(guò)程的耦合作用,其數(shù)學(xué)模型為隨機(jī)混合系統(tǒng)(stochastic hybrid system).例如遺傳開(kāi)關(guān)現(xiàn)象[1-2],電力驅(qū)動(dòng)的細(xì)胞內(nèi)運(yùn)輸[3-4],隨機(jī)神經(jīng)元網(wǎng)絡(luò)[5-9]和隨機(jī)離子通道等等[10-12].此時(shí),由于連續(xù)狀態(tài)和離散狀態(tài)之間的耦合特性,使得混合系統(tǒng)比隨機(jī)微分方程更難處理.雖然可以通過(guò)采用擴(kuò)散近似方法得到相應(yīng)混合系統(tǒng)的擴(kuò)散過(guò)程近似,但是眾所周知,如果考慮亞穩(wěn)態(tài)動(dòng)力學(xué)行為,除了在接近分岔極限[13]的特殊情況之外,該方法通常會(huì)產(chǎn)生顯著的誤差[1,14-15].此外,恰恰是由于開(kāi)通道數(shù)的離散特性,混合系統(tǒng)比連續(xù)系統(tǒng)更符合實(shí)際物理模型.因此,為了揭示實(shí)際模型的內(nèi)在動(dòng)力學(xué)行為,直接處理原始的混合系統(tǒng)是必要且本質(zhì)的.

很多混合系統(tǒng)均可表現(xiàn)出可激性(excitability)特征,此時(shí)即便是很小的隨機(jī)擾動(dòng)也可能對(duì)其動(dòng)力學(xué)行為產(chǎn)生深遠(yuǎn)的影響,如自發(fā)性動(dòng)作電位(spontaneous action potential)現(xiàn)象.一般意義下,如果偏離靜息態(tài)的一個(gè)小擾動(dòng)能夠?qū)е缕潆娢辉诜祷仂o息態(tài)之前出現(xiàn)一個(gè)大的偏差,即稱產(chǎn)生了一個(gè)動(dòng)作電位[16].通常,自發(fā)性動(dòng)作電位可以分為兩個(gè)階段:初始階段和興奮階段.初始階段是指系統(tǒng)在噪聲驅(qū)動(dòng)下從靜息態(tài)到運(yùn)動(dòng)到某個(gè)閾值(threshold)的過(guò)程,因此可視為噪聲誘導(dǎo)的離出(exit)事件或逃逸事件.這個(gè)階段,由于噪聲強(qiáng)度較弱,將花費(fèi)相當(dāng)長(zhǎng)的時(shí)間.一旦達(dá)到閾值,系統(tǒng)將沿著確定性流演化,直至返回靜息態(tài),此為興奮階段.與前者相比,興奮階段是暫態(tài)的,其運(yùn)動(dòng)軌跡僅受確定性動(dòng)力學(xué)運(yùn)動(dòng)規(guī)律控制.因此,只要能夠準(zhǔn)確刻畫初始階段的離出行為,便可清晰地描述整個(gè)動(dòng)作電位.

事實(shí)上,Freidlin 和Wentzell[17]所建立的大偏差理論(large deviation theory)已被廣泛應(yīng)用于分析弱噪聲對(duì)具有多種極限集—-不動(dòng)點(diǎn)[18-19]、極限環(huán)[20-21]和奇怪吸引子[22-23]等系統(tǒng)的長(zhǎng)期動(dòng)力學(xué)行為的影響.依據(jù)大偏差理論中定義的作用(量)泛函(action functional)的概念,通過(guò)它在任一絕對(duì)連續(xù)函數(shù)處的值能夠得到任一隨機(jī)軌線通過(guò)該函數(shù)附近時(shí)其難易程度的量度.根據(jù)自吸引子出發(fā)所有路徑集合上作用泛函的駐值,便可以確定最優(yōu)的離出路徑,即系統(tǒng)以最大的概率沿著該路徑移動(dòng).

基于大偏差理論,人們對(duì)噪聲誘導(dǎo)的逃逸問(wèn)題進(jìn)行了大量的研究,但主要結(jié)果限于閾值具有明確定義的情況,一般為系統(tǒng)鞍點(diǎn)的穩(wěn)定流形,此時(shí)最優(yōu)路徑則為連接吸引子和鞍點(diǎn)的異宿軌道[24].然而遺憾的是,很多系統(tǒng)在相空間中并沒(méi)有鞍點(diǎn).因此為了研究這些離出問(wèn)題,必須定義合適的閾值.近期,Chen等[25]和Khovanov 等[26]通過(guò)選擇虛擬分界線作為FitzHugh-Nagumo 系統(tǒng)的閾值,描述了其離出方案和有限噪聲效應(yīng).Carson 等[27]研究了通道漲落引起的Hodgkin-Huxley 模型的放電事件,并通過(guò)一維郎之萬(wàn)(Langevin)近似計(jì)算了其平均放電頻率.James 等[28]將該模型簡(jiǎn)化為只保留簡(jiǎn)單的鈉離子通道噪聲的模型,并計(jì)算出平均首次離出時(shí)間的精確表達(dá)式.基于Hodgkin-Huxley 模型同時(shí)保留鈉離子和鉀離子通道噪聲的另一種簡(jiǎn)化系統(tǒng),此時(shí)其相應(yīng)的確定性模型只有一個(gè)穩(wěn)定的不動(dòng)點(diǎn),Brooks 等[12]發(fā)現(xiàn)了擬環(huán)(quasi-cycle)的出現(xiàn).Newby[11]利用WKB 近似方法對(duì)同一系統(tǒng)進(jìn)行了研究,得到了離出路徑滿足的Hamilton 表述.本文以他們的工作為基礎(chǔ),重點(diǎn)研究了該模型在自定義閾值下的自發(fā)性動(dòng)作電位,并討論了通道數(shù)對(duì)可激系統(tǒng)動(dòng)力學(xué)行為的影響.

迄今為止,用于計(jì)算擬勢(shì)的主要方法是WKB 近似,它被廣泛應(yīng)用于研究擴(kuò)散過(guò)程、生滅控制方程和隨機(jī)混合系統(tǒng)中的逃逸事件.此外,本文采用的數(shù)值方法包括幾何最小作用量法(geometric minimum action method,GMAM)[29]和有序逆風(fēng)法(ordered upwind method,OUM)[30].它們均基于幾何作用概念,即將作用泛函的時(shí)間參數(shù)化轉(zhuǎn)化為弧長(zhǎng),從而可將關(guān)于無(wú)窮時(shí)間的積分簡(jiǎn)化為有限長(zhǎng)度的積分.GMAM 可用于求解具有固定端點(diǎn)的最優(yōu)路徑和它的作用泛函,而OUM 可用于計(jì)算離散化的相空間節(jié)點(diǎn)上的擬勢(shì).

本文具體組織如下:第2 節(jié)介紹了隨機(jī)混合Morris-Lecar 模型,第3 節(jié)中給出了極值路徑滿足的輔助Hamilton 系統(tǒng)和改進(jìn)的混合系統(tǒng)的Monte Carlo模擬算法.通過(guò)定義虛擬分界線為閾值,第4 節(jié)研究了系統(tǒng)的自發(fā)性放電現(xiàn)象.在第5 節(jié)中,通過(guò)改變鈉離子和鉀離子通道數(shù),我們?cè)敿?xì)討論了不同噪聲組合對(duì)系統(tǒng)擬勢(shì)和最優(yōu)路徑的影響.結(jié)論見(jiàn)第6 節(jié).

1 隨機(jī)模型

本文所研究的模型最初由Hodgkin 和Huxley[31]提出并用于描述魷魚(yú)巨軸突的膜電位,Newby[11]將其簡(jiǎn)化為只考慮單通道鈉離子和鉀離子的情形,出于論文的完整性考慮,我們簡(jiǎn)述其模型結(jié)構(gòu).首先,膜電位v的演化受控于離子電流和外加刺激電流Iapp,即

其中,n=0,1,2,···,N和m=0,1,2,···,M分別表示鈉離子和鉀離子的開(kāi)通道數(shù),函數(shù)fi(v)=gi(vi-v),i=Na,K,leak 分別決定了鈉離子和鉀離子電流和漏泄電流.鈉離子和鉀離子開(kāi)閉狀態(tài)的切換取決于

其中,aNa(v)=exp[4(γNav+κNa)],bNa=1,aK(v)=exp(γKv+κK)和bK(v)=exp(?γKv-κK).由此可見(jiàn),電位v的演化受控于方程(1).在兩次相鄰跳躍之間,其值是確定性的,但會(huì)通過(guò)影響轉(zhuǎn)移率反過(guò)來(lái)修正隨機(jī)的跳躍時(shí)間,因此它們之間是耦合的.

根據(jù)耦合過(guò)程的馬氏性,其概率密度函數(shù)p(v,m,n,t)滿足下面 Chapman-Kolmogorov (CK)方程[9]

其中的算子定義如下

而轉(zhuǎn)移率分別為

假設(shè)βNa和M均很大且滿足此外,可將w=m/M視為第二個(gè)狀態(tài)變量,在弱噪聲極限ε→0 時(shí),系統(tǒng)恢復(fù)為確定性的Morris-Lecar 方程[32]

其中,x∞(v)={1+tanh[2(γNav+κNa)]}/2,w∞(v)=[1+tanh(γKv+κK)]/2 和τw(v)=2/βKcosh(γKv+κK).βK和βNa決定了鉀離子和鈉離子通道開(kāi)閉的快慢,通常有βK?βNa,意味著鈉離子通道的開(kāi)閉遠(yuǎn)快于鉀離子.

在研究隨機(jī)模型之前,需要先考察由方程(6)定義的確定性系統(tǒng)的動(dòng)力學(xué)行為,方程中的參數(shù)見(jiàn)附錄.如圖1(a)所示,棕色和粉紅色曲線分別代表v零傾線和w零傾線,它們相交于一個(gè)全局穩(wěn)定的平衡點(diǎn)SN,圖中的小箭頭指明了向量場(chǎng)的變化趨勢(shì).此外依據(jù)四階龍格?庫(kù)塔法,圖中給出了5 條具有代表性的確定性路徑,即圖中的綠色曲線.通過(guò)觀察,可以發(fā)現(xiàn),雖然所有的路徑最終都收斂于SN,但從不同的初始點(diǎn)出發(fā),它們之間存在著本質(zhì)上的區(qū)別.其中一些迅速而直接地衰減到SN,而另外的則先向右移動(dòng),直至到達(dá)v零傾線的右支,然后向左轉(zhuǎn),這樣便產(chǎn)生了一個(gè)動(dòng)作電位.圖中的紅色粗體曲線稱為虛擬分界線(ghost separatrix)[33],它是通過(guò)v零傾線的頂點(diǎn)P時(shí)間后向積分得到的.可以看出,初始點(diǎn)位于此分界線左側(cè)的路徑,運(yùn)動(dòng)的振幅較小,而右側(cè)對(duì)應(yīng)于大振幅運(yùn)動(dòng).在圖1(a)中,右側(cè)三條路徑的初始位置之間的差異僅為0.01,很難區(qū)分,但在這些路徑的中間位置上則出現(xiàn)了宏觀上的有限距離.因此,這條曲線起到了擬閾值的作用[16].

圖1 向量場(chǎng)、確定性軌線和隨機(jī)軌線.棕色和粉紅色曲線分別表示v 零傾線和w 零傾線,SN 為不動(dòng)點(diǎn),綠色曲線是確定性路徑,紅色粗體線指代分界線,小箭頭代表向量場(chǎng)的方向.一條代表性隨機(jī)軌道表示為灰色鋸齒線,里面包含一次自發(fā)性動(dòng)作電位Fig.1 Vector field,deterministic and stochastic trajectories.Brown and pink curves indicate the v-nullcline and w-nullcline respectively.SN denotes the fixed point,the green curves show the deterministic paths,the red bold line suggests the ghost separatrix,and the small arrows represent directions of vector field.A representative stochastic trajectory has been denoted by gray zigzag line,containing a spontaneous action potential

在圖1(b)中,灰色鋸齒形曲線是一條具有代表性的表示自發(fā)性動(dòng)作電位的隨機(jī)路徑.由于SN 的吸引性和較弱的噪聲強(qiáng)度,系統(tǒng)將長(zhǎng)期在SN 附近徘徊.隨后在某個(gè)偶然的時(shí)刻,它會(huì)躍過(guò)閾值,然后再次沿著確定的路徑返回到SN.這整個(gè)過(guò)程構(gòu)成一個(gè)自發(fā)性動(dòng)作電位.

2 理論分析與數(shù)值模擬方法

2.1 WKB 近似

隨機(jī)混合Morris-Lecar 的WKB 近似的結(jié)果已由Newby 導(dǎo)出[11].出于論文完整性的考慮,我們簡(jiǎn)述其主要結(jié)果以便于后面的理論分析及數(shù)值計(jì)算,并給出幾點(diǎn)重要的評(píng)述.

在弱噪聲條件下,假設(shè)平穩(wěn)概率密度具有下面近似的形式

這里C(v,w)是前因子項(xiàng),本文不予討論.r(n|v,w)是已知v,w后n的條件概率分布.W(v,w)稱為擬勢(shì),它度量了系統(tǒng)在隨機(jī)激勵(lì)下抵達(dá)相空間中某一位置的困難程度.將方程(7)代入平穩(wěn)CK 方程并合并所有ε 的最高階項(xiàng)系數(shù)得到Hamilton-Jacobi 方程

根據(jù)特征線法可知,這一方程的解由一個(gè)輔助Hamilton 系統(tǒng)生成

通過(guò)積分可以得到每條特征線上的擬勢(shì),即

對(duì)于WKB 近似的形式及結(jié)果,我們認(rèn)為有必要給出以下三個(gè)重要的評(píng)述:

首先,相比于擴(kuò)散過(guò)程,方程(7)的形式多了一項(xiàng)r(n|v,w),這是因?yàn)橄到y(tǒng)中存在離散變量.由此,概率密度的歸一化條件不能像擴(kuò)散過(guò)程那樣進(jìn)行直接積分且其積分值等于1,而必須先對(duì)離散變量求和然后再對(duì)連續(xù)變量進(jìn)行積分.

其次,若在方程(9)中令p=0,則輔助Hamilton方程恢復(fù)為確定性的Morris-Lecar 系統(tǒng),這意味著動(dòng)量p度量了離出的極值路徑與確定性流之間的距離.就是說(shuō),它刻畫了噪聲在驅(qū)動(dòng)系統(tǒng)離出的過(guò)程中所扮演的角色.因此,與隨機(jī)微分系統(tǒng)的情形類似,可將動(dòng)量值視為噪聲影響的一個(gè)指示器.

第三,作者已嚴(yán)格證明了混合系統(tǒng)一定不滿足細(xì)致平衡條件的事實(shí)[14],這一結(jié)論與擴(kuò)散過(guò)程的相關(guān)結(jié)論有著本質(zhì)的區(qū)別.對(duì)于混合系統(tǒng),離出路徑不會(huì)出現(xiàn)可以簡(jiǎn)單地視為確定性流的時(shí)間反向的情況.于是,尖點(diǎn)和焦散線等奇異性[34]的出現(xiàn)會(huì)更為普遍,這源于由極值路徑所追溯的流形出現(xiàn)折疊的特性,并因此導(dǎo)致擬勢(shì)的多值性[35].

2.2 Monte Carlo 模擬

基于Gillespie 算法[36]的Monte Carlo 模擬的基本原理是通過(guò)計(jì)算隨機(jī)跳躍時(shí)間和狀態(tài)轉(zhuǎn)移來(lái)演化樣本軌跡.Bressloff[8]和Newby[11]分別描述了他們算法的主要步驟,并應(yīng)用于神經(jīng)網(wǎng)絡(luò)和Morris-Lecar等模型.然而,他們的方法得以實(shí)施均有一個(gè)前提條件,即要求方程˙v=Iion(v,m,n)對(duì)固定的m,n有顯式解,因此很難推廣到更復(fù)雜的系統(tǒng).在此我們提出一種適用于幾乎所有隨機(jī)混合系統(tǒng)的新算法,并且此方法能夠計(jì)算出足夠精確的跳躍時(shí)間,從而可以近似到任何期望的精度.

為了計(jì)算下一次跳躍時(shí)間,我們需要考慮四種可能的轉(zhuǎn)移:n→n±1 和m→m±1.在開(kāi)通道數(shù)的兩次相鄰跳躍之間,電位的演化受控于確定性動(dòng)力學(xué)此處求解的主要的困難來(lái)源于轉(zhuǎn)移率和電位的耦合關(guān)系.為了解決這個(gè)問(wèn)題,定義

從而可以得到具有初始條件v(0)=v0和G(0)=0 的微分方程組

只需要提供時(shí)間步長(zhǎng)h和兩個(gè)均勻分布隨機(jī)變量Y和r,當(dāng)下列條件滿足時(shí)

跳躍時(shí)間和電位的估計(jì)就可由τ≈kh和v(τ)≈v(kh)給出.為了確定所發(fā)生的轉(zhuǎn)移是4 種可能中的哪一種,變量a定義為滿足下列不等式的最小整數(shù)

其中a=1,2,3,4 分別對(duì)應(yīng)于以下4 種轉(zhuǎn)移行為:n→n+1,n→n-1,m→m+1,m→m-1.至此可以看出,通過(guò)將函數(shù)F的積分定義為一個(gè)新的變量,不需要像前人一樣解析計(jì)算出v的表達(dá)式,而是將其轉(zhuǎn)化為一個(gè)微分方程組來(lái)同時(shí)進(jìn)行數(shù)值積分,因此這一方法是通用的.

3 不同閾值下的自發(fā)性動(dòng)作電位

如前所述,自發(fā)性動(dòng)作電位可分為兩個(gè)階段:初始階段和興奮階段.由于噪聲強(qiáng)度較弱,系統(tǒng)將在吸引子附近游蕩很長(zhǎng)時(shí)間,隨后,它將離開(kāi)靜息態(tài),并在某個(gè)隨機(jī)時(shí)間穿過(guò)閾值.一旦到達(dá)閾值,它將沿著確定性軌跡返回到吸引子.此時(shí)便意味著一個(gè)動(dòng)作電位出現(xiàn)了.事實(shí)上,初始階段的本質(zhì)是:在給定閾值的條件下,由噪聲誘導(dǎo)的從吸引域中的逃逸事件(離出問(wèn)題).

迄今為止,我們應(yīng)該能夠意識(shí)到閾值的重要性！因?yàn)樗鼪Q定了什么樣的動(dòng)力學(xué)行為意味著一次放電.如果相空間中存在鞍點(diǎn),那么其穩(wěn)定流形自然起到閾值的作用,最優(yōu)路徑就是連接吸引子和鞍點(diǎn)的路徑,即異宿軌道.穿過(guò)鞍點(diǎn)后,系統(tǒng)會(huì)沿著鞍點(diǎn)的不穩(wěn)定流形運(yùn)動(dòng).然而遺憾的是,對(duì)于不具有鞍點(diǎn)的系統(tǒng)(如本文中的模型),很難定義這樣一個(gè)閾值.盡管如此,本文的系統(tǒng)與具有鞍點(diǎn)的系統(tǒng)有很多相同的動(dòng)力學(xué)特征.為了了解這一事實(shí),我們需要首先選擇一個(gè)具有規(guī)范定義的閾值來(lái)觀察初始逃逸事件.考慮到第2 節(jié)中關(guān)于模型的確定性動(dòng)力學(xué)的分析,虛擬分界線在自發(fā)性動(dòng)作電位中起著重要作用.如果系統(tǒng)受噪聲的驅(qū)動(dòng)下穿過(guò)分界線,那么系統(tǒng)之后將幾乎沿著確定性軌跡繞一個(gè)大圈回到不動(dòng)點(diǎn),形成一個(gè)電位的尖峰.從這一點(diǎn)來(lái)看,把此分界線看成閾值是合理的.

在考察離出路徑之前,我們先來(lái)考察系統(tǒng)的擬勢(shì).依據(jù)OUM 可計(jì)算出相關(guān)區(qū)域的擬勢(shì),其結(jié)果由圖2 所示的等高線和三維圖給出.值得注意的是,該數(shù)值方法是通過(guò)從不動(dòng)點(diǎn)開(kāi)始逐次更新相鄰節(jié)點(diǎn)的值來(lái)計(jì)算規(guī)則網(wǎng)格上的擬勢(shì),其結(jié)果必然是擬勢(shì)的全局最小值.對(duì)于此系統(tǒng),用800×800 的網(wǎng)格覆蓋區(qū)域[?0.5,1.8]×[0,0.6]足夠精確.圖2(b)所示,SN 附近存在一個(gè)很深的勢(shì)阱,在勢(shì)阱右上方有一個(gè)很大的平坦區(qū)域,擬勢(shì)幾乎保持不變,隨后在遠(yuǎn)處迅速增長(zhǎng).我們注意到,圖2(a)中在w=0 處有一個(gè)擬勢(shì)的截痕,這是因?yàn)閣表示鉀離子的開(kāi)通道數(shù)的比例,不可能是負(fù)的.

圖2 利用OUM 在800×800 的網(wǎng)格上計(jì)算的擬勢(shì)Fig.2 Quasi-potential calculated by OUM on a 800×800 grid

由輔助Hamilton 系統(tǒng)(9)的解所刻畫的軌跡,在t0→?∞時(shí)起始于SN 的條件下,在四維擴(kuò)展相空間中張成二維不穩(wěn)定拉格朗日流形.這些路徑在坐標(biāo)空間上的投影決定了極值路徑的模式.在擴(kuò)展的相空間中,由于增加了兩個(gè)不穩(wěn)定維度,原來(lái)的穩(wěn)定不動(dòng)點(diǎn)轉(zhuǎn)化為鞍點(diǎn).為了描述極限路徑的模式,我們使用作用量圖(action plot)[23]方法求解輔助Hamilton系統(tǒng)的初值問(wèn)題.在二維拉格朗日流形上選取均勻分布于以SN 為中心半徑為0.001 的圓周上的點(diǎn)為初始點(diǎn),并且每條極值路徑可由初始點(diǎn)的角位置參數(shù)化,與擬勢(shì)的特征線方程(10)一起積分,直到閾值后終止.其中幾條具有代表性的極值路徑如圖3(a)所示,它們起始時(shí)極為接近,在到達(dá)分界線前散開(kāi).當(dāng)極值路徑觸碰到分界線時(shí),該點(diǎn)的擬勢(shì)值便被記錄,因此我們能夠得到分界線上所有點(diǎn)的擬勢(shì),其與w的函數(shù)關(guān)系由圖3(b)給出.此外,利用插值原理,前面的OUM 方法的結(jié)果也可以用于計(jì)算閾值處的擬勢(shì),其結(jié)果如圖3(b)中的黑線所示.可以看出,兩種方法的結(jié)果吻合得相當(dāng)好.

圖3 選取分界線為閾值Fig.3 The separatrix chosen as threshold

GMAM 是另一種計(jì)算擬勢(shì)的方法.其主要思想是在兩個(gè)端點(diǎn)給定的條件下計(jì)算連接此兩點(diǎn)的最優(yōu)路徑及其上的作用泛函.我們將一個(gè)端點(diǎn)選在SN,另一個(gè)位于分界線,通過(guò)不斷改變分界線上的點(diǎn)以計(jì)算不同位置的擬勢(shì)值.為了提高計(jì)算效率,這一方法的誤差會(huì)比前兩個(gè)結(jié)果稍大一些,但仍在可接受范圍內(nèi).從圖3(b)中可以看出,分界線上的擬勢(shì)在w≈0.023 (對(duì)應(yīng)于左圖中的S0點(diǎn))附近存在明顯的最小值.對(duì)于雙穩(wěn)態(tài)系統(tǒng),鞍點(diǎn)一般是邊界上擬勢(shì)的最小值點(diǎn),因此離出行為總發(fā)生在鞍點(diǎn).對(duì)于本文的模型,盡管不存在鞍點(diǎn),但是S0有著類似于鞍點(diǎn)的作用.

為了驗(yàn)證我們之前的理論結(jié)果,現(xiàn)采用歷程概率分布方法模擬漲落路徑的分布.這一概念由Dykman[37]首先提出,歷程概率分布Ph(x,t|xf,tf)被定義為三時(shí)刻轉(zhuǎn)移概率Ph(xf,tf;x,t|xi,ti)與兩時(shí)刻轉(zhuǎn)移概率Ph(xf,tf|xi,ti)之比,初始時(shí)刻ti被設(shè)為?∞,從而ti和xi均可從表達(dá)式中去掉.這一分布包含了系統(tǒng)到達(dá)xf之前關(guān)于系統(tǒng)演化的所有信息.根據(jù)大偏差理論,系統(tǒng)會(huì)以壓倒性的概率通過(guò)最優(yōu)路徑的鄰域,因此歷程概率分布將在最優(yōu)路徑附近達(dá)到峰值.為了觀察系統(tǒng)的動(dòng)力學(xué)行為,我們固定終點(diǎn)時(shí)間為tf=0 并向前追溯.基于第2 節(jié)中改進(jìn)的Monte Carlo算法,可以得到眾多從SN 出發(fā)直至閾值的隨機(jī)路徑.根據(jù)核密度估計(jì),最終得到歷程概率分布的結(jié)果如圖4(a)所示.自這些軌跡的樣本集合中,提取每一時(shí)刻歷程概率分布的最大值位置,如圖4(b)所示,其結(jié)果與GMAM 計(jì)算的最優(yōu)路徑的理論值吻合得相當(dāng)好.此外,這一結(jié)果證實(shí)了我們之前的觀點(diǎn),即離出位置集中在S0上.

前面的研究表明,將自SN 出發(fā)到分界線的最優(yōu)路徑作為動(dòng)作電位的初始階段是合理的.然而,為了說(shuō)明以分界線作為閾值的合理性,還應(yīng)該觀察系統(tǒng)在越過(guò)分界線后是否存在大振幅運(yùn)動(dòng),并且其結(jié)果是否與模擬結(jié)果一致.為了觀察系統(tǒng)的放電特性,將離出路徑進(jìn)一步延伸,允許最優(yōu)路徑觸及v零傾線的右分支,看看本文的分析是否成立?

圖4 理論計(jì)算與數(shù)值結(jié)果的最優(yōu)路徑對(duì)比Fig.4 The comparison for the optimal paths between the theoretical and numerical results

首先,通過(guò)對(duì)隨機(jī)軌道樣本進(jìn)行積分直到v零傾線右支,經(jīng)統(tǒng)計(jì)可得其歷程概率分布如圖5(a)所示.提取該分布的峰頂線,所得結(jié)果如圖5(b)所示,曲線在穿過(guò)分界線的時(shí)候非常接近S0.另一方面,由GMAM 給出的極值路徑的理論解自SN 開(kāi)始,至S0結(jié)束,進(jìn)一步可延伸至v零傾線的N0點(diǎn).從圖中可以看出,理論解與數(shù)值模擬結(jié)果非常吻合,這說(shuō)明將分界線作為閾值不僅對(duì)初始階段是合理的,對(duì)興奮階段同樣合理.圖5(b)中所示的小箭頭表示確定性向量場(chǎng)的方向,可以看出,在穿過(guò)分界線后,最優(yōu)路徑與確定性流的方向很接近.我們換個(gè)角度觀察動(dòng)量隨時(shí)間變化情況,即圖5(c).可以看出最優(yōu)路徑穿越閾值后,動(dòng)量值迅速衰減.如前所述,動(dòng)量本質(zhì)上是對(duì)漲落路徑偏離確定性流的程度的度量,動(dòng)量的衰減這一現(xiàn)象正說(shuō)明離出路徑接近于確定性軌線這一事實(shí).此外,擬勢(shì)的值源于動(dòng)量的積分,而擬勢(shì)在分界線右邊變化緩慢,這也與OUM 的結(jié)果一致.因此,當(dāng)系統(tǒng)到達(dá)N0后,幾乎沿確定性軌跡返回到SN,發(fā)生一次完整的放電現(xiàn)象.

圖5 初始階段和興奮階段的最優(yōu)路徑及其動(dòng)量、作用量Fig.5 The optimal path as well as the momentum,the action for initial and excitation phase

圖6 另一種閾值選擇Fig.6 Another choice for threshold

根據(jù)圖2 中擬勢(shì)的形狀可知,若系統(tǒng)能夠躍出勢(shì)阱,則其后續(xù)運(yùn)動(dòng)在平坦區(qū)域的能量消耗很小.因此,Newby 選擇勢(shì)阱與平坦區(qū)域的接合處作為閾值,將從勢(shì)阱中的逃逸視為起始階段[11].其具體操作如下:首先,如圖6 所示,相空間中存在一個(gè)尖點(diǎn)(cusp point,CP),延伸出兩條焦散線(caustics)和一條切換線(switching line,SL).這一現(xiàn)象的出現(xiàn)是由于拉格朗日流形在擴(kuò)展相空間中的折疊結(jié)構(gòu),使得其在坐標(biāo)空間上的投影具有多值性[21].而擬勢(shì)取全局最小值,由此會(huì)導(dǎo)致擬勢(shì)在切換線上的不可微性,出現(xiàn)一條折痕.因此,到達(dá)SL 兩側(cè)兩點(diǎn)的最優(yōu)路徑在拓?fù)渖鲜峭耆煌?于是,可以選擇SL 作為閾值的一部分,另一部分由尖點(diǎn)的擬勢(shì)等高線(contour line,CL)提供,分別用粗體的棕色和綠色曲線表示.可以看出最終進(jìn)入平坦區(qū)域的極值路徑在穿過(guò)CL 時(shí)幾乎都位于某個(gè)S 點(diǎn)附近.

在結(jié)束這一節(jié)之前,我們比較一下這兩種閾值選擇的差別.首先,以分界線為閾值更多的是關(guān)注確定性系統(tǒng)的行為,而第二種主要考慮的是擬勢(shì)的結(jié)構(gòu).另外,前者只承認(rèn)大振幅運(yùn)動(dòng)為一次放電現(xiàn)象,這實(shí)際上更符合動(dòng)力學(xué)的觀點(diǎn).更重要的是,由于鞍點(diǎn)的擬勢(shì)在整個(gè)邊界上取最小值,使得在弱噪聲極限情況下,最優(yōu)路徑必然會(huì)通過(guò)鞍點(diǎn),因此鞍點(diǎn)狀態(tài)在逃逸問(wèn)題中起著至關(guān)重要的作用.但是,由于選擇擬勢(shì)等高線作為閾值,CL 上的每個(gè)點(diǎn)具有相同的概率,這削弱了S 的特殊意義.綜合來(lái)看,我們認(rèn)為將分界線作為閾值更為合理.最后,需要強(qiáng)調(diào)的是,這兩種選擇的平均首次離出時(shí)間幾乎是相同的,因?yàn)樗鼈兌荚跀M勢(shì)的平坦區(qū)域,具有幾乎相等的擬勢(shì).

4 不同噪聲組合下的影響

迄今為止,我們研究了隨機(jī)混合Morris-Lecar 系統(tǒng)對(duì)于不同閾值在M=N=40 的情況下自發(fā)性放電現(xiàn)象.本節(jié),我們將分析鈉離子和鉀離子噪聲的各種不同組合,即不同的M和N的值將如何影響離出方案和擬勢(shì)形狀.

首先為探索M和N量級(jí)的影響,令M=N.可以發(fā)現(xiàn),當(dāng)具有不同M和N值和相同初始點(diǎn)時(shí),極值路徑完全一致.然而,當(dāng)使用GMAM 計(jì)算分界線上SN 到S0的最優(yōu)路徑時(shí),擬勢(shì)卻發(fā)生了明顯的變化,如圖7(a)所示,擬勢(shì)與M,N的量級(jí)呈嚴(yán)格的正比關(guān)系.換句話說(shuō),M,N的量級(jí)線性地決定了勢(shì)阱的深度.這一事實(shí)很好理解,因?yàn)槿绻鸐和N的值很大,那么由m和n的加減所帶來(lái)的變化對(duì)v的導(dǎo)數(shù)影響很小,于是便導(dǎo)致了更大的逃逸難度.

圖7 M 和N 不同組合下擬勢(shì)的對(duì)比Fig.7 Comparison for quasi-potential with various combinations of M and N

另一方面,通過(guò)研究M和N值的不同組合,以確定何種漲落對(duì)整個(gè)逃逸過(guò)程的貢獻(xiàn)更大.首先觀察分界線上的擬勢(shì)與M,N的不同組合之間的關(guān)系.將N分別取為1,10,40,400,擬勢(shì)與M=1,10,40,400之間的4 條折線關(guān)系如圖7(b)所示.可以看出,如果固定N,則擬勢(shì)將隨著M的增大而顯著增大.當(dāng)然,也存在一個(gè)N=1 時(shí)的特殊情況,即擬勢(shì)自M=40到M=400 時(shí)增大得不明顯.不過(guò),這一現(xiàn)象也是可以理解的,若M與N的比值過(guò)大,鉀離子的漲落便會(huì)足夠小.

另一個(gè)角度,通過(guò)固定M我們可以分析鈉離子通道噪聲對(duì)擬勢(shì)的影響.若M=1,則當(dāng)N取不同值時(shí),擬勢(shì)值幾乎無(wú)法區(qū)分.此外,當(dāng)M的值較大時(shí),在N

事實(shí)上,根據(jù)表1 中的數(shù)據(jù),我們從另一個(gè)方面也可以發(fā)現(xiàn)這些特性.對(duì)角線顯示了前文中的線性關(guān)系.此外,當(dāng)我們觀察另一條對(duì)角線就會(huì)發(fā)現(xiàn)另一個(gè)更有趣的性質(zhì)—M和N的乘積保持不變而擬勢(shì)存在明顯的最大值,即,存在一個(gè)能使?jié)q落強(qiáng)度最小的M與N的最優(yōu)比例.我們可以利用這一性質(zhì)控制噪聲強(qiáng)度的大小.另外可以發(fā)現(xiàn),在這個(gè)最優(yōu)值的條件下M是N的數(shù)倍,這也從另一方面證實(shí)了鉀離子通道噪聲的主導(dǎo)作用.

表1 M 和N 不同組合下分界線上S0的擬勢(shì)Table 1 Quasi-potential of S0on the separatrix versus various combinations of M and N

接下來(lái),通過(guò)調(diào)整M和N的相對(duì)大小,我們觀察鈉離子和鉀離子通道噪聲如何分別影響離出路徑.對(duì)不同的M,將8 種情況分為N=1 和N=40 兩組,它們的最優(yōu)路徑分別如圖8 所示.可以看出,最優(yōu)路徑的動(dòng)力學(xué)行為由M和N的比值決定,與量級(jí)無(wú)關(guān).如果M遠(yuǎn)大于N,則最優(yōu)路徑接近水平線.當(dāng)M大于或者幾乎等于N時(shí),最優(yōu)路徑位于v零傾線以下.若初始時(shí)M=N,則增大M會(huì)顯著改變最優(yōu)路徑的形狀,而N的增大所引起的改變很小.類似于擬勢(shì)與通道噪聲強(qiáng)度比值的關(guān)系,最優(yōu)路徑也主要取決于鉀離子通道噪聲.最后我們指出,在具有快慢尺度的擴(kuò)散過(guò)程里存在類似的現(xiàn)象,即調(diào)整快慢方向的噪聲強(qiáng)度比例可以控制最優(yōu)路徑水平或彎曲.實(shí)際上,快變方向的噪聲較強(qiáng)時(shí),系統(tǒng)能夠直接克服水平勢(shì)壘從而沿水平方向離出,而快變方向噪聲較弱時(shí)無(wú)法克服,因此只能沿著快變方向零傾線離出,表現(xiàn)為彎曲的最優(yōu)路徑.在文獻(xiàn)[38]中對(duì)一相對(duì)簡(jiǎn)單的快慢尺度擴(kuò)散過(guò)程給出了理論解釋,本文研究的模型雖為混合系統(tǒng),但仍具有快慢特性,因此出現(xiàn)了類似的現(xiàn)象.

圖8 M 和N 不同組合下最優(yōu)路徑的對(duì)比Fig.8 Comparison for the optimal paths with various combinations of M and N

5 結(jié)論

本文研究了因?yàn)橄到y(tǒng)缺乏鞍點(diǎn)狀態(tài)而沒(méi)有明確閾值的隨機(jī)混合Morris-Lecar 系統(tǒng)的自發(fā)性動(dòng)作電位.通過(guò)定義虛擬分界線為擬閾值,可以將一次放電行為分解為一個(gè)長(zhǎng)時(shí)間的初始階段和一個(gè)暫態(tài)的興奮階段.前者可以看成是噪聲誘導(dǎo)的逃逸事件,而后者則幾乎會(huì)沿著確定性的路徑在大范圍運(yùn)動(dòng)后回到靜息態(tài).通過(guò)三種不同的數(shù)值方法,即GMAM、OUM和作用量圖方法,計(jì)算得出閾值上的擬勢(shì),發(fā)現(xiàn)結(jié)果存在明顯的最小值.由此可以肯定離出位置幾乎集中于這個(gè)點(diǎn)附近,因此可以斷定該點(diǎn)發(fā)揮著類似鞍點(diǎn)的作用.通過(guò)改進(jìn)的蒙特卡羅方法我們對(duì)歷程概率分布進(jìn)行了模擬,驗(yàn)證了我們前面使用理論方法得到的最優(yōu)路徑和離出位置的結(jié)果.另外,針對(duì)前人所給出的另一種閾值取法,即切換線和尖點(diǎn)擬勢(shì)等高線的組合,我們比較了兩種選擇的優(yōu)劣性,并認(rèn)為,從動(dòng)力學(xué)的角度來(lái)講,前者更具有實(shí)際物理意義.

最后,研究表明,M,N的量級(jí)嚴(yán)重影響勢(shì)阱的深度,并且它們之間呈現(xiàn)嚴(yán)格的線性關(guān)系.此外,M和N的不同組合對(duì)最優(yōu)路徑的模式和擬勢(shì)有著顯著影響,鉀離子通道噪聲在自發(fā)性動(dòng)作電位過(guò)程中起主導(dǎo)作用.研究還發(fā)現(xiàn),存在一個(gè)M和N的最優(yōu)比值,能使?jié)q落強(qiáng)度達(dá)到最小.

附錄

參數(shù)值

vNa=3.7,gNa=0.22,vK=?0.9,gK=0.4,vleak=?0.36,gleak=0.1,βK=0.04,Iapp=0.06,γNa=1.22,κNa=?1.188,γK=0.8,κK=?0.8,M=40,N=40.