耦合CPG神經(jīng)網(wǎng)絡(luò)的節(jié)律分析

袁泉, 徐潔瓊, 陳慧瑩

(1.廣西大學(xué) 土木建筑工程學(xué)院, 廣西 南寧 530004；2.廣西大學(xué) 數(shù)學(xué)與信息科學(xué)學(xué)院, 廣西 南寧 530004)

0 引言

在動(dòng)物與環(huán)境的互動(dòng)中,包括許多基本的重復(fù)行為,如呼吸、行走、爬行、游泳、咀嚼、抓撓等,這些重復(fù)行為由不同神經(jīng)元集群的節(jié)律運(yùn)動(dòng)的協(xié)調(diào)作用而產(chǎn)生。生物學(xué)家認(rèn)為動(dòng)物的節(jié)律運(yùn)動(dòng)是神經(jīng)中樞的一種自激產(chǎn)生的振蕩行為,產(chǎn)生節(jié)律活動(dòng)的神經(jīng)環(huán)路被稱為中樞模式發(fā)生器(central pattern generator,CPG)[1]。1911年,Brown[1]提出脊髓運(yùn)動(dòng)CPG的半中心振蕩模型,該模型被廣泛應(yīng)用于各種重要運(yùn)動(dòng)節(jié)律的產(chǎn)生和控制中。例如,Rubin等[2]利用幾何動(dòng)力系統(tǒng)方法分析突觸耦合的神經(jīng)元網(wǎng)絡(luò)集群節(jié)律的產(chǎn)生,表明快速抑制耦合將產(chǎn)生同步節(jié)律。

研究表明,節(jié)律運(yùn)動(dòng)是CPG的一種自激產(chǎn)生的振蕩行為,但隨著外界驅(qū)動(dòng)的改變,運(yùn)動(dòng)節(jié)律模式會(huì)相適應(yīng)地改變。不少學(xué)者逐漸研究外界刺激對(duì)CPG運(yùn)動(dòng)節(jié)律模式的變化。Whittington等[3]在實(shí)驗(yàn)中觀察到抑制耦合的2個(gè)神經(jīng)元網(wǎng)絡(luò)的同步振蕩是由谷氨酸受體激活引起的。Traub等[4]通過實(shí)驗(yàn)與理論分析再次驗(yàn)證谷氨基酸將影響神經(jīng)網(wǎng)絡(luò)的同步,并發(fā)現(xiàn)神經(jīng)元的放電頻率呈單調(diào)變化,只有當(dāng)突觸電流在合理范圍時(shí),突觸傳遞將為同步振蕩提供一個(gè)合適的機(jī)制。Xu等[5]基于相關(guān)系數(shù)討論了同步振蕩與同步程度之間區(qū)別與聯(lián)系,發(fā)現(xiàn)在不完全同步振蕩時(shí)可能會(huì)發(fā)生弱相位同步。

而谷氨酸一般可以刺激2種突觸受體:一種是alpha-amino-3-hydroxy-5-methylisoxazole-4-propionic acid(APMA)；另一種是N-methyl-D-aspartic(NMDA)。Talpalar等[6]通過對(duì)轉(zhuǎn)基因小鼠脊髓神經(jīng)元加入NMDA,發(fā)現(xiàn)神經(jīng)元的興奮程度對(duì)NMDA有很強(qiáng)的依賴性,還比較了在增加與刪除V0類中間神經(jīng)元時(shí)對(duì)CPG運(yùn)動(dòng)節(jié)律的差異,最終表明神經(jīng)元振蕩頻率由NMDA濃度調(diào)節(jié)。楊永霞等[7]通過控制興奮性耦合神經(jīng)網(wǎng)絡(luò)的發(fā)放模式和呼吸節(jié)律,分析了CPG的動(dòng)力學(xué)機(jī)制。孫燕等[8]討論了中樞神經(jīng)系統(tǒng)中一類神經(jīng)元的表達(dá)模式與釋放特性。Molkov等[9]通過將中間神經(jīng)元轉(zhuǎn)化為相應(yīng)的連接路徑,觀察興奮性刺激藥物在刪除不同連接路徑時(shí)對(duì)肢體動(dòng)物的運(yùn)動(dòng)步態(tài)的影響,得知運(yùn)動(dòng)頻率由神經(jīng)元興奮和連合通路的參數(shù)所控制,刪除不同類型的路徑時(shí),神經(jīng)元網(wǎng)絡(luò)也表現(xiàn)出不同的同步與交替現(xiàn)象,這一結(jié)果也與文獻(xiàn)[10]相符。

本文通過對(duì)2個(gè)耦合的CPG神經(jīng)網(wǎng)絡(luò)引入刺激藥物NMDA,這一藥物具體表現(xiàn)為應(yīng)用離子通道中的參數(shù)反轉(zhuǎn)電位,通過控制NMDA的濃度觀察肢體動(dòng)物在不同連接路徑的步態(tài)模式。數(shù)值仿真軟件為XPPAUT,積分步長0.01,使用的四階龍格庫塔算法。

1 數(shù)學(xué)模型

本文研究的肢體動(dòng)物左右肢協(xié)調(diào)步態(tài)運(yùn)動(dòng)的神經(jīng)網(wǎng)絡(luò)模型由文獻(xiàn)[9]給出, 該模型由2個(gè)CPG構(gòu)成,分別控制左、右肢的節(jié)律運(yùn)動(dòng)。每個(gè)CPG由相互抑制耦合的屈肌中心和伸肌中心組成。左、右CPG中心之間主要考慮3條路徑:興奮性V3中間神經(jīng)元調(diào)節(jié)左屈肌、右屈肌中心的相互興奮；抑制性V0D中間神經(jīng)元調(diào)節(jié)2個(gè)中心相互抑制；興奮性V0V中間神經(jīng)元也參與這些中心的相互抑制,同時(shí)促進(jìn)相互交替。該模型示意圖如圖1所示。

圖1 簡化模型示意圖

外界驅(qū)動(dòng)改變CPG神經(jīng)網(wǎng)絡(luò)的運(yùn)動(dòng)節(jié)律模式。外界驅(qū)動(dòng)不僅對(duì)興奮路徑產(chǎn)生影響,同時(shí)對(duì)神經(jīng)元的動(dòng)力學(xué)行為產(chǎn)生影響，因此,本文在文獻(xiàn)[9]的基礎(chǔ)上對(duì)神經(jīng)網(wǎng)絡(luò)模型進(jìn)行改進(jìn),得到的模型為

(1)

(2)

其中:

Inapi=gnapim∞(vi)hnap(vi-vna),Ili=gl(vi-vl),Iappi=gapp(vi-vapp),

式中:變量v、h、m分別為電壓、鈉、鉀離子門通道打開的概率；i=1,2,3,4分別為各中心對(duì)應(yīng)的指標(biāo)；Cm,為膜電容,Inap、Iapp、Il、IsynE、IsynI分別為持續(xù)鈉電流、應(yīng)用電流、漏電流、興奮性突觸電流和抑制性突觸電流；g代表相應(yīng)的電導(dǎo)。突觸電流如圖1 所示,這里僅寫出細(xì)胞1的突觸電流,其他的類似。

IsynI1=[adf(v3)+b2f(v2)]gsynI(v1-vsynI),IsynE1=[a3f(v3)+avf(v4)]gsynE(v1-vsynE),

式中:vmin、vmax分別為達(dá)到閾值和飽和狀態(tài)的電壓；gsynE、gsynI分別為興奮性和抑制性突觸通道的最大電導(dǎo)；vsynE、vsynI分別是興奮性和抑制性突觸的反轉(zhuǎn)電位；ad、av、a3、b1、b2分別為連接路徑V0D、V0V、V3、Ini1、Ini2的突觸權(quán)重；其中f(v)表示規(guī)范化輸出電位,關(guān)于其他變量的描述在文獻(xiàn)[9]中可找到。為體現(xiàn)細(xì)胞在外界藥物刺激下運(yùn)動(dòng)節(jié)律模式發(fā)生的變化情況,在模型中引入?yún)?shù)α,并定義

vapp=vapp0(1-α)。

(3)

另外,為了更好地模擬刺激性藥物對(duì)興奮路徑的影響,需要對(duì)這些連接進(jìn)行突觸加權(quán),以增加參數(shù)α所代表的NMDA濃度,因此記a3=a30(1+β3α)，av=av0(1+βvα)。

2 動(dòng)力學(xué)分析

2.1 單個(gè)神經(jīng)元的動(dòng)力學(xué)行為

此時(shí)只考慮單個(gè)細(xì)胞的運(yùn)動(dòng)節(jié)律,以左屈為例,在方程(1)中不需要再考慮突觸電流,因此,這一情形的膜電位是由Inap、Il、Iapp三者所決定。單細(xì)胞模型在(v,h)相平面的周期軌和零值線如圖2所示。圖中v零值線(三次型)與h零值線(單調(diào)遞減)的交點(diǎn)位于v零值線的中支,這也是系統(tǒng)的平衡點(diǎn),并且是不穩(wěn)定的,閉軌線是模型的一個(gè)穩(wěn)定極限環(huán),關(guān)于這類研究在文獻(xiàn)[11-12]已有介紹。右側(cè)為周期軌對(duì)應(yīng)的膜電位波形圖,周期T=1 283 ms,單位周期中的4個(gè)部分分別對(duì)應(yīng)細(xì)胞在基準(zhǔn)驅(qū)動(dòng)上的靜息狀態(tài)、向上跳躍狀態(tài)、活躍狀態(tài)以及向下跳躍狀態(tài)。根據(jù)神經(jīng)元的運(yùn)動(dòng)狀態(tài),周期軌可分為4個(gè)部分:①靜息狀態(tài)；②向上跳躍狀態(tài)；③激活狀態(tài)；④向下跳躍狀態(tài)。

圖2 單細(xì)胞模型在(v,h)相平面的周期軌與零值線

在靜息和活躍狀態(tài),解曲線分別靠近v零值線的左、右分支；而向上、向下跳躍狀態(tài)出現(xiàn)時(shí)在軌跡上到達(dá)左右拐點(diǎn)。由幾何奇異攝動(dòng)理論,向上、向下跳躍狀態(tài)滿足方程的快子系統(tǒng),而靜息和活躍狀態(tài)對(duì)應(yīng)慢子系統(tǒng)的解。細(xì)胞從①區(qū)出發(fā),經(jīng)過Ts時(shí)間后,到達(dá)左拐點(diǎn)處附近；在②區(qū)神經(jīng)元向上跳躍,為瞬時(shí)活動(dòng)；在③區(qū),細(xì)胞逐漸活躍,經(jīng)過Ta時(shí)間后到達(dá)④區(qū)；④區(qū)為向下跳躍狀態(tài)也為瞬時(shí)活動(dòng),最終完成一個(gè)周期運(yùn)動(dòng)。

Iapp隨時(shí)間的變化情況如圖3所示。從圖中可見,Iapp具有周期性,而構(gòu)成這一項(xiàng)的主要參數(shù)是vapp與gapp,vapp的分岔圖如圖4所示。從圖中可見,隨著vapp的增加,平衡點(diǎn)從最初的穩(wěn)定變?yōu)椴环€(wěn)定,在HB處穩(wěn)定性發(fā)生改變,再次回到穩(wěn)定狀態(tài),并從HB產(chǎn)生穩(wěn)定的極限環(huán),說明有一個(gè)Hopf 分岔出現(xiàn),而且是超臨界Hopf分岔。極限環(huán)由穩(wěn)定變?yōu)椴环€(wěn)定,SN為周期軌上的鞍結(jié)點(diǎn)。本文主要研究的是神經(jīng)元在周期軌的運(yùn)動(dòng)模式,故取vapp∈(-5.541,91.8)。由前所述,每個(gè)單一的神經(jīng)元可處于以下3種狀態(tài):靜息、振蕩、強(qiáng)直。在該范圍內(nèi),神經(jīng)元是振蕩的。再根據(jù)α與vapp的關(guān)系(圖5),得出本文所考慮的α的范圍為(0,2.8]。

圖3 Iapp隨時(shí)間的變化情況

圖4 單神經(jīng)元模型關(guān)于參數(shù)vapp的分岔圖

圖5 外界參數(shù)α與vapp的關(guān)系

2.2 完整神經(jīng)網(wǎng)絡(luò)模型的動(dòng)力學(xué)行為

現(xiàn)在討論完整的神經(jīng)元網(wǎng)絡(luò)模型,本文的完整模型是由2個(gè)CPG耦合的神經(jīng)網(wǎng)絡(luò)構(gòu)成,經(jīng)過突觸連接,構(gòu)成一個(gè)8維的非線性微分方程組。首先通過相平面分析,分析4個(gè)神經(jīng)元的動(dòng)力學(xué)行為。

4個(gè)中心的零值線、周期軌與動(dòng)作電位,如圖6所示。由圖中可見,左右中心的周期軌和動(dòng)作電位具備相同的特征。CPG在受外界刺激情況下,盡管左右中心的電位振幅有所不同,但相位差保持一致,這就是相位同步,表現(xiàn)為圖7中的相位差恒為0。而同側(cè)中心運(yùn)動(dòng)模式較為復(fù)雜,左側(cè)中心在受外界刺激時(shí)運(yùn)動(dòng)初期不太穩(wěn)定,在靜息狀態(tài)存在小范圍的擾動(dòng),但擾動(dòng)會(huì)隨著時(shí)間而消失,后面將說明這種擾動(dòng)也會(huì)隨著NMDA濃度的增加而消失。由于2個(gè)CPG具有相同的特性,因此只需討論其中一個(gè)在外界刺激下的運(yùn)動(dòng),現(xiàn)著重比較左側(cè)中心的變化。

圖6 4個(gè)中心的零值線,周期軌與動(dòng)作電位

圖7 中心相位差

左側(cè)中心在α=0.1,α=1.0,α=2.8時(shí)的零值線、周期軌與動(dòng)作電位,如圖8所示。隨著α增加,左屈肌與右屈肌始終是完全同步的,而左側(cè)中心在一定時(shí)刻以后對(duì)不同的參數(shù)都始終趨向一個(gè)非零數(shù),說明在完整模型中左側(cè)中心在不同外界刺激下受α的影響比較大,其相位差的幅度逐漸減少,而且由對(duì)稱逐漸轉(zhuǎn)化為非對(duì)稱。屈肌與伸肌中心的周期軌和動(dòng)作電位具備相同的特征,即細(xì)胞1與細(xì)胞3,細(xì)胞2與細(xì)胞4兩組中心的具有對(duì)稱性。從圖1中就能看出,這也是文獻(xiàn)[1、13]所研究的基于持續(xù)性鈉電流的半中心振蕩。

同側(cè)中心有許多共同點(diǎn),h是與通道門有關(guān)的狀態(tài)變量,外界藥物不會(huì)影響各離子通道打開的概率,另外,在與h有關(guān)的方程中不顯含α,無論α取何值,都不會(huì)影響到h零值線。相比之下,2個(gè)中心的v零值線受α的影響比較大,都是隨著α的增加,其左膝值的h軸坐標(biāo)下降,平衡點(diǎn)與左膝值的變化相同,不管α在所給定的區(qū)間上如何變動(dòng),各中心的零值線始終保持三部分。并由此確定了α的區(qū)間,在圖8中,細(xì)胞1與細(xì)胞2對(duì)α取不同值都始終存在閉軌,這與圖4中參數(shù)vapp的分岔的情況相符。其次,隨著α增加,周期變小,意味著在外界藥物的刺激下,細(xì)胞活動(dòng)變得更加興奮,這與右欄周期放電的頻率增加相符,且尖峰電位也隨之下降。

同側(cè)的電位活動(dòng)也存在著較大差異,前者在α=0.1,α=1.0處的平衡點(diǎn)位于零值線的左支,在α=2.8處位于中支,說明細(xì)胞1在外界藥物刺激的作用下,顯示處于興奮狀態(tài),在藥物達(dá)到一定劑量后開始變?yōu)檎袷?而細(xì)胞2始終是處于振蕩狀態(tài),而不受外界藥物的影響,此外,細(xì)胞1在剛接觸到外界刺激時(shí)靜息狀態(tài)表現(xiàn)出很大差異,在這活動(dòng)期間還出現(xiàn)了一個(gè)短暫的向上跳躍狀態(tài),稱為閾下放電,為了更好表現(xiàn)二者在周期軌與尖峰電位受α影響程度的差異,將區(qū)間細(xì)分,取步長0.3,同側(cè)中心的規(guī)范化周期與尖峰電位如圖9所示。

圖9 同側(cè)中心的規(guī)范化周期與尖峰電位

前面已經(jīng)說明左右中心是同步的,而同側(cè)是交替的,從圖8(a)、8(c)可以看出作圖的周期幾乎一致,僅僅是在α=2.8處稍有不同,確切的說是這一點(diǎn)處的細(xì)胞活動(dòng)反常,在此之前,同側(cè)中心的活動(dòng)周期與動(dòng)作電位都是隨著α單調(diào)下降,而在這點(diǎn)卻增加,這是因?yàn)?當(dāng)外界藥物達(dá)到一定劑量之后,細(xì)胞的反應(yīng)不再像之前那樣強(qiáng)烈,而是逐漸到達(dá)一種平衡狀態(tài),其相軌線具體表現(xiàn)為靠近平衡點(diǎn),而在α=2.8時(shí)已經(jīng)逐漸接近這種臨界狀態(tài),此前的2個(gè)周期已經(jīng)變?yōu)?個(gè)次周期,這就是周期2峰放電。換言之,將觀察到的實(shí)際周期除以2依舊滿足之前單調(diào)下降的特征。圖8(c)、8(d)中細(xì)胞在該點(diǎn)處的動(dòng)作電位也符合這一特征。

2.3 刪除V0D

在完整的耦合模型中主要考慮了V3中間神經(jīng)元、抑制性V0D中間神經(jīng)元與興奮性V0V中間神經(jīng)元,而外界刺激藥物主要是影響了細(xì)胞的興奮性程度,路徑V0D的權(quán)重并不與α直接相關(guān),因此,此處考慮刪除路徑V0D以后細(xì)胞隨α的變化情況。構(gòu)成Iapp電流的還有基準(zhǔn)應(yīng)用逆轉(zhuǎn)勢電位,刪除V0D路徑之后,這一參數(shù)也可能會(huì)影響模型的穩(wěn)定性,對(duì)該參數(shù)進(jìn)行分岔分析與圖4基本相同,依舊存在超臨界的Hopf分岔,其中分岔點(diǎn)為vapp0=168。在這一模型中仍取vapp0=-50。在刪除V0D后,周期軌上的鞍結(jié)分岔點(diǎn)消失。

刪除路徑V0D的模型下周期軌及動(dòng)作電位,如圖10所示。各細(xì)胞的活動(dòng)情況與完整模型中的并無很大差異,周期軌與零值線的位置也幾乎一致,放電周期與尖峰電位也沒有太大改變。可以說這一情形下的各細(xì)胞的運(yùn)動(dòng)節(jié)律模式與完整模型幾乎一致,完全符合這條路徑的屬性,路徑V0D是由抑制性V0D中間神經(jīng)元改變得到的,而抑制性神經(jīng)元與外界刺激藥物非直接相關(guān),所引入的參數(shù)α無法影響這條路徑,故關(guān)于這條路徑V0D的周期變化情況與前者類似。屈肌與伸肌中心仍然是完全同步。在完整模型中,同側(cè)中心在一定時(shí)刻后的相位差逐漸趨于一個(gè)固定值,而且隨著參數(shù)的改變這一固定值也會(huì)減小,而在刪除V0D模型下同側(cè)中心的相位差不再穩(wěn)定在某一固定值,而是沿著v=0周期性的波動(dòng),當(dāng)參數(shù)增加到臨界值時(shí),變化幅度變小,而此時(shí)的振幅并無關(guān)聯(lián),說明同側(cè)中心隨著參數(shù)的逐漸增加開始出現(xiàn)同步現(xiàn)象。

圖10 刪除路徑V0D的模型下周期軌及動(dòng)作電位

2.4 刪除V0V

前面已經(jīng)討論了刪除與外界刺激無關(guān)的路徑V0D的情況,剩下的2條興奮性路徑V3、興奮性V0V均與外界藥物刺激直接相關(guān),現(xiàn)討論刪除路徑V0V的情形。

左側(cè)屈肌與左側(cè)伸肌中心對(duì)所考慮的α是交替的,但交替周期隨著α的增加而減少,使運(yùn)動(dòng)頻率增加。這點(diǎn)在動(dòng)作電位與輸出函數(shù)上表現(xiàn)出較大的相關(guān)性,因?yàn)樗x的輸出函數(shù)就是對(duì)電壓的規(guī)范化。最大的差別表現(xiàn)在向下跳躍狀態(tài),左屈中心受到外界刺激時(shí)在向下跳躍的過程中出現(xiàn)較長的閾下電位,而這一狀態(tài)隨著時(shí)間的增加持續(xù)時(shí)間越來越短,直至消失,當(dāng)外界藥物增加到一定程度后,又重新開始出現(xiàn)新的放電,此時(shí)的放電不同于前者,從輸出函數(shù)圖上表現(xiàn)為周期2峰放電。

屈肌與伸肌中心始終表現(xiàn)為完全同步(圖11)。而左屈與左伸的相位差與其動(dòng)作電位間的變化情況幾乎相同,尤其是在向下跳躍狀態(tài)依舊表現(xiàn)出閾下放電,這一行為在相位差的圖中漸趨于0,說明左屈與左伸在交替中仍有出現(xiàn)相位同步,只不過隨著α的增加,相位同步逐漸消失。

圖11 反相同步與完全同步

3 結(jié)論

不同的外界刺激可以引起放電節(jié)律模式發(fā)生改變,也會(huì)導(dǎo)致神經(jīng)元產(chǎn)生不同的生理效應(yīng),刺激幅度偏離需求值以及周期偏大或偏小都無法對(duì)疾病產(chǎn)生好的效果,甚至?xí)又夭∏閇14]，因此,了解不同的節(jié)律模式與外界刺激之間的關(guān)系將有助于解釋外界刺激導(dǎo)致的疾病所產(chǎn)生的機(jī)理,在臨床治療時(shí),刺激的相關(guān)參數(shù)控制在合理的范圍內(nèi)顯得尤為重要。本文考慮了一個(gè)基于持續(xù)性鈉電流和應(yīng)用電流的2個(gè)CPG耦合的神經(jīng)網(wǎng)絡(luò)模型,研究了外界刺激藥物NMDA對(duì)運(yùn)動(dòng)步態(tài)節(jié)律的影響。首先對(duì)應(yīng)用反轉(zhuǎn)電位vapp作分岔分析,由此確定外界參數(shù)α的范圍,本文章后面部分都是基于這一范圍對(duì)模型進(jìn)行分析。通過研究完整網(wǎng)絡(luò)中各細(xì)胞的放電活動(dòng)發(fā)現(xiàn)左右中心是同步的,而且這一現(xiàn)象不會(huì)受到NMDA濃度變化的影響,同側(cè)中心則表現(xiàn)為交替活動(dòng),但這種交替行為容易受到外界藥物的干擾,在交替與同步中相互轉(zhuǎn)化。接著又分別討論了刪除抑制性路徑V0D和交替興奮性路徑V0V的模型,并相互比較,觀察到左右中心在這2種情形下都表現(xiàn)出與完整模型相同的特征,即神經(jīng)元始終是同步放電,不受外界刺激的影響,周期的變化情況也大致相同,主要的差別是在同側(cè)中心,動(dòng)作電位及相位差都隨著參數(shù)的變化而變化,本文也是將重點(diǎn)放在了這一部分討論,這一結(jié)論符合文獻(xiàn)[15]中關(guān)于復(fù)雜網(wǎng)絡(luò)魯棒性的研究。關(guān)于具體的實(shí)驗(yàn)數(shù)據(jù)見文獻(xiàn)[6],實(shí)驗(yàn)分別在野生型和轉(zhuǎn)基因去雄小鼠脊髓中進(jìn)行。關(guān)于多個(gè)CPG耦合的運(yùn)動(dòng)機(jī)制是如何整合起來產(chǎn)生穩(wěn)定肢體間的相互作用,這在很大程度上仍是一個(gè)未知的領(lǐng)域,仍將是一個(gè)活躍的研究領(lǐng)域。