Hindmarsh-Rose神經元陣列自發產生螺旋波的研究?

汪芃 李倩昀 唐國寧

(廣西師范大學物理科學與技術學院,桂林 541004)

1 引 言

自從Winfree于1972年首次在化學系統中觀察到螺旋波以來[1],人們在各種系統中,無論在可激發介質還是在一般介質上,都觀察到了螺旋波[2?12].眾所周知,心臟出現螺旋波電信號會導致危險的心律失常,如心動過速和心室纖維性顫動[7],從而危及生命.心臟中的螺旋波一般不會自發消失,除非螺旋波漫游出邊界消失[7],否則會持續存在,需要人為干預才能消除,所以螺旋波的控制受到人們極大的關注[13,14].研究發現[8,15]:無論在藥物引起的振蕩中,還是睡眠狀態下,在哺乳動物大腦皮層中都會自發出現螺旋波,只是螺旋波的壽命很短;癲癇發作時,在大腦皮層中也會循環出現螺旋波.螺旋波的存在可以在介觀尺度上組織和調節皮層的群體活動,但是它們在哺乳動物大腦皮層中產生的機制和潛在功能仍然不清楚,需要進一步研究.

為了了解螺旋波在神經元系統中產生的機制,人們從理論和實驗上進行了對噪聲誘發螺旋波的研究.因為在宏觀世界中被人們稱為“噪聲”的隨機漲落是普遍存在的,雖然噪聲給人們的第一印象是有害的,但是噪聲和非線性系統在一定條件下會發生隨機共振現象,促進非線性系統響應的增強.1998年,Jung等[16]首次在培養的神經膠質細胞網絡中觀察到由網絡噪聲誘發的螺旋化學波,這種螺旋波在噪聲作用下也會自發消失;1999年,García-Ojalvo等[17]在可激發介質中也觀察到噪聲誘發的螺旋波.這些發現引起人們對噪聲誘發螺旋波的極大興趣[18?23].因為實際神經系統出現噪聲是不可避免的,例如在大腦中,一個神經元通常與許多神經元有耦合,各神經元并不是處于同步狀態,而是處于非同步狀態,因此這種眾多的非近鄰耦合的總效果相當于在只考慮近鄰耦合的神經元上施加一種噪聲擾動.噪聲可以誘發螺旋波,也可以使同質介質系統(由相同神經元組成的系統)處于混沌態,那么系統從隨機的初相位分布態演化是否能自發形成螺旋波?迄今為止這個問題仍缺乏研究,對這一問題進行研究會有助于人們了解大腦中螺旋波產生的機制.

本文采用Hindmarsh-Rose(HR)神經元模型[24]研究了由同質神經元組成的二維陣列系統從具有隨機相位分布的初態演化.我們發現,當單個HR神經元處于一周期態時,只要系統的耦合強度在一定范圍內,系統都可以自發出現螺旋波,適當選取耦合強度和初態,系統還能自發出現圓形波等其他無螺旋波的態,以及出現同步振蕩和振蕩死亡現象.當單個HR神經元處于二周期態時,系統自發出現螺旋波的能力大為減少,這些結果對了解大腦中螺旋波的形成、癲癇和振蕩死亡產生機制具有積極意義.

2 模 型

本文采用1984年Hindmarsh和Rose提出的三變量的HR神經元模型來構造一個二維神經元陣列系統,該系統的動力學方程如下[24]:

式中x代表神經元的細胞膜電位;y代表與內電流相關的恢復變量;z代表與鈣離子激活的鉀離子電流相關的慢變調節電流;a,b,c,d,r,s和x0為HR神經元模型參數;Iext為外部刺激電流;g為耦合強度;角標i,j=1,2,···,N代表N ×N個耦合的HR神經元.為了控制神經元系統的動力學行為,參考文獻[20]中相關模型參數的取值,本文固定取a=1.0,b=3.0,c=1.0,d=5.0,r=0.006,s=4.0,x0= ?1.6,N=200,Iext和g為可調參數,通過適當選擇Iext的值,可使得單個神經元(在情況下方程(1)—(3)描述的系統)處于一周期、二周期等周期態,或者使單個神經元處于混沌態,本文只考慮單個神經元處于一周期和二周期的情況,數值解方程組采用歐拉法,時間步長取Δt=0.02,使用無流邊界條件,每次數值模擬時間長度為12000單位時間.

系統變化狀態為混沌時,對應的方差值變化也就較大.而當系統出現有序狀態時,部分神經元振蕩趨于一致,此時方差值也趨于穩定.因此,方差可以較好地描述系統的狀態變化.為了描述二維神經元陣列系統狀態變化情況,引入系統方差:

當σ無規變化時,系統處于混沌態;如果σ有規律變化,系統處于有序態;如果σ趨于0,系統處于同步態.

3 數值模擬結果

系統的演化結果與系統初態和耦合強度有關,由于系統初態很多,為了比較不同初態演化效果,規定系統初態中各神經元膜電位x＞0數量與神經元總數之比為ρ.本文通過以下方式產生系統演化的初態:在耦合強度g=0和給定Iext的情況下,給系統的每一個神經元賦隨機初值,讓系統演化2000個時間單位后,找ρ滿足給定值的系統狀態作為g/=0和給定Iext情況下系統的初態,該初態具有混沌初相位分布,一般ρ值在2%—8%之間.在下面的數值模擬中,我們都會選擇不同ρ值的不同初態進行數值模擬研究,在沒有特別指出的情況下,給出的斑圖都是系統演化12000單位時間后的斑圖.考慮到耦合強度很小時,系統不會自發出現單個螺旋波,耦合強度取值大于0.1.

3.1 一周期態的結果

固定Iext=1.315,這時單個神經元處于一周期態,研究在相同ρ值下系統從不同初態演化是否自發出現螺旋波.考慮到取不同ρ值都得到類似結果,這些結果與三種類型的初態有關.下面取ρ=5%和三個典型的初態,分別記為初態1、初態2和初態3.為了衡量初態的差異,將單個神經元狀態的周期變化分成5個階段,如圖1(a)所示,圖中相鄰兩個黑點之間的狀態變化為一個階段,圖中數字為階段的編號.在初態1下,處于階段1至階段5的神經元數與總數之比γ分別為γ1=82.75%,γ2=11.03%,γ3=0.938%,γ4=5.28%,γ5=0%;同理得到初態2的各階段占比分別為γ1=0.045%,γ2=4.425%,γ3=4.323%,γ4=91.207%,γ5=0%;初態3的各階段占比分別為γ1=71.195%,γ2=14.058%,γ3=2.13%,γ4=10.788%,γ5=1.825%.這三個典型態初態的特點是:初態1的γ1=82.75%,處于階段1的神經元占絕大多數;初態2的γ4=91.207%,處于階段4的神經元占絕大多數;初態3下的γ5=1.825%,處于階段5的神經元占比不為0,處于比較均勻的分布;圖1(b)—(d)分別給出三個初態的x變量斑圖,可以看出,系統初態中各神經元具有無規的初相位分布.

圖1(a)單個HR神經元的吸引子在z-x平面上的投影和(b)初態1、(c)初態2、(d)初態3對應的x變量斑圖Fig.1.(a)The attractor of the single HR neuron in(z,x)projection and the patterns of variable x for(b)the initial state 1,(c)the initial state 2,and(d)the initial state 3.

下面先使用初態1,研究系統的演化.圖2給出了不同耦合強度下方差隨時間的變化,圖3給出了與圖2對應參數下x變量的斑圖,圖2和圖3結果對比如下.

1)在無耦合下單個神經元周期振蕩,導致系統方差規則變化.

2)在0.1＜g＜2.9的情況下,系統經過暫態后一般自發出現多個螺旋波和螺旋波對,出現螺旋波與初態無關.耦合強度不同,螺旋波和螺旋波對的位置和大小不相同(即斑圖不同),螺旋波的數量與耦合強度有關,這表明螺旋波的出現具有偶然性,因此不同初態和不同耦合強度時,斑圖一般都不一樣.一般有這樣的規律:當耦合強度較小時,神經元簇狀放電的峰個數少,螺旋波的波臂較細,因此螺旋波尺寸小,數量就多,系統不存在單個螺旋波;增大耦合強度,神經元簇狀放電的峰個數增多,螺旋波波臂變粗,系統中螺旋波尺寸就大,數量相應就少,所以適當選擇耦合強度可以得到單螺旋波.通常多個螺旋波和螺旋波對相互作用導致規則或無規則斑圖在系統中重復出現,方差接近規則變化,如圖2(b)—(e)所示.這些結果表明系統出現單螺旋波的概率比較小,尺寸小的螺旋波出現的概率比較大,這些結果與實驗觀察到的結果一致.實驗觀察到[16]:在背景網絡噪聲(由紅藻氨酸濃度調節)作用下膠質細胞網絡中自發出現螺旋波,小尺寸螺旋波出現的概率大,大尺寸螺旋波出現的概率小.

3)在g≥3.0的情況下,系統一般不會出現螺旋波,一般出現圓形波或環形波,它們相互作用導致系統出現平面波、不規則的曲線波等斑圖.圖3(e)所顯示的亮斑就會形成圓形波,但是這種圓形波不穩定,會隨時間變化,在該狀態下,存在所有神經元滿足x＜0,因此系統方差表現為無規律大幅度振蕩變化.在數值模擬中,我們在g=2.84,4.3,4.5,4.8情況下分別觀察到單個螺旋波,如圖3(d)和圖3(f)所示.

4)從圖3可以看出,螺旋波的波臂存在明顯的無規律明暗分布,這是因為在螺旋波態下,神經元簇狀放電有多個峰,而不是一周期態的單個峰,而且峰的數量及高度和間隔都在無規則變化,這是導致在螺旋波態下系統方差隨時間變化不是很有規律的原因.

圖2 在使用初態1和不同耦合強度下方差隨時間的變化(a)g=0;(b)g=0.2;(c)g=1.2;(d)g=2.2;(e)g=2.84;(f)g=3.0Fig.2.The time evolution of the variance for different values of coupling strength when the initial state 1 is applied:(a)g=0;(b)g=0.2;(c)g=1.2;(d)g=2.2;(e)g=2.84;(f)g=3.0.

圖3 在使用初態1和不同耦合強度下x變量的斑圖(a)g=0.2;(b)g=1.2;(c)g=2.2;(d)g=2.84;(e)g=3.0;(f)g=4.3Fig.3.Pattern of the variable x for different values of coupling strength when the initial state 1 is applied:(a)g=0.2;(b)g=1.2;(c)g=2.2;(d)g=2.84;(e)g=3.0;(f)g=4.3.

圖4 在使用初態2和不同耦合強度下方差隨時間的變化(a)g=0.2;(b)g=0.5;(c)g=1.0;(d)g=1.4;(e)g=1.5;(f)g=2.5Fig.4.The time evolution of the variance for different values of coupling strength when the initial state 2 is applied:(a)g=0.2;(b)g=0.5;(c)g=1.0;(d)g=1.4;(e)g=1.5;(f)g=2.5.

下面使用初態2,研究系統的演化.圖4給出了在不同耦合強度下方差隨時間的變化,圖5給出了與圖4對應參數下x變量的斑圖.從圖4和圖5可以看出:當0.1＜g≤0.5時,系統很容易出現螺旋波和螺旋波對,只是隨耦合強度增加,系統出現小螺旋波和螺旋波對所用時間不斷增加;當0.5＜g≤1.4時,系統不容易出現螺旋波,只是在適當的耦合強度下才出現螺旋波,參見圖5(c),這時系統較容易出現不規則的平面波,參見圖5(c)和圖5(d),方差無規律小幅變化;當g≥1.5時,方差大部分時間接近0,當神經元狀態處于階段2和3時,方差才比較大,如圖5(e)和圖5(f)所示,這表明系統整體出現間歇式全局同步振蕩,而且隨著耦合強度的增加,系統整體同步程度逐漸增強.神經元容易同步的原因是,絕大部分神經元處于階段4,各神經元膜電位的值差別不大,通過膜電位耦合就容易實現全局同步.系統出現間歇同步的原因是,當神經元狀態處于階段4和5時,盡管方差接近0,但是神經元并沒有達到精確全局同步,x變量的差的最大值在千分之二左右,而且系統的x變量的弱非均勻分布呈中心在對角的靶波狀,這樣當神經元狀態處于階段2時,神經元出現依次被激發,導致系統方差又增大,當神經元處于階段4時,神經元又處于接近精確全局同步,方差又接近0.

圖5 在使用初態2和不同耦合強度下x變量的斑圖(a)g=0.2;(b)g=0.5;(c)g=1.0;(d)g=1.4;(e)g=1.5;(f)g=2.5Fig.5.Pattern of the variable x for different values of coupling strength when the initial state 2 is applied:(a)g=0.2;(b)g=0.5;(c)g=1.0:(d)g=1.4;(e)g=1.5;(f)g=2.5.

最后使用初態3研究系統的演化,同樣得到:當耦合強度在[0.2,0.9]范圍時,系統都很容易產生螺旋波,且隨著耦合強度的增加,螺旋波或螺旋波對的數量逐漸減少,與前兩種初態演化結果一致,當g=0.9時,系統出現單個螺旋波;當g≥1.0時,系統快速進入振蕩死亡,所有神經元振蕩方式演化到不動點x=?1.31742,y=?7.687,z=1.13032,如圖6所示.產生振蕩死亡的原因是,不動點靠近階段5和階段1,當初態中的神經元處于階段5的比率足夠大時,就可以出現振蕩死亡.為了證明是處于階段5的神經引起振蕩死亡,通過取不同ρ值(在2%—8%范圍內)的初態,發現只要,不論初態是γ1最大或還是γ4最大都會出現振蕩死亡現象,而且γ5越大,達到振蕩死亡需要的耦合強度就越小.

圖6 在使用初態3和g=1.0情況下系統中某個神經元在z-x平面上的相圖Fig.6.Phase diagram of a neuron in the system on the plan(z,x).The initial state 3 and g=1.0 are applied.

以上是取ρ=5%時得到的結果,當ρ取其他值時,同樣得到:當耦合強度在小于臨界值gc的一個范圍內取值時,系統容易出現螺旋波和螺旋波對,增大耦合強度,螺旋波和螺旋波對的數量將減少,偶爾出現單個螺旋波;當耦合強度大于臨界值gc時,系統出現三種不同的動力學行為,分別對應三類初態.系統從第一類初態演化,一般不容易出現螺旋波,而是出現圓形波,只有適當選擇耦合強度才能出現螺旋波;不同耦合強度會得到不同斑圖,即使出現單螺旋波斑圖,其波頭位置也可能不同;系統從第二、三初態演化,在大的耦合強度下同樣分別觀察到了同步振蕩和振蕩死亡現象.

3.2 二周期結果

當固定Iext=1.60,這時單個神經元處于二周期態.選擇不同ρ值,研究系統從不同初態演化,觀察到兩種演化結果,分別對應兩類初態,一類初態不能使系統自發出現螺旋波,另一類初態能使系統出現螺旋波,但是只能在一定的耦合強度范圍內才能出現螺旋波,出現螺旋波的耦合強度范圍比一周期情況大為縮小.圖7給出了與圖1類似的結果,這里仍取ρ=5%.將神經元二周期振蕩分成4個階段,如圖7(a)所示,圖7(b)和圖7(c)是兩個典型的初態,分別記為初態4和初態5.初態4的γ1=50.515%,γ2=18.895%,γ3=6.89%,γ4=23.7%,其分布特點是不夠均勻,因為γ3偏小;初態5的γ1=34.19%,γ2=26.13%,γ3=10.07%,γ4=29.645%,其分布特點是比較均勻.

圖8給出了系統從初態4和5演化得到的部分結果.系統從初態4演化的研究結果是:當0.2≤g≤0.4時,系統只能產生小波,最終形成向中心傳播方形波,即反靶波,如圖8(a)和圖8(b)所示;當g≥0.5時,系統形成圓形波,如圖8(c)所示.產生圓形波的原因是,處于階段3的神經元較少,相當于振蕩的循環過程,少了一個過程,而螺旋波、靶波形成必須有四個過程,所以系統演化過程不容易形成螺旋波和靶波.

圖7(a)單個HR神經元的吸引子在z-x平面上的投影和(b)初態4、(c)初態5對應的x變量斑圖Fig.7.(a)The attractor of the single HR neuron in(z,x)projection and the patterns of variable x for(b)the initial state 4 and(c)the initial state 5.

圖8 不同初態和不同耦合強度下x變量的斑圖(a)初態4,g=0.3;(b)初態4,g=0.4;(c)初態4,g=0.8;(d)初態5,g=0.2;(e)初態5,g=1.4;(f)初態5,g=1.6Fig.8.The patterns of variable x for different coupling strengths and initial states:(a)The initial state 4 and g=0.3;(b)the initial state 4 and g=0.4;(c)the initial state 4 and g=0.8;(d)the initial state 5 and g=0.2;(e)the initial state 5 and g=1.4;(f)the initial state 5 and g=1.6.

系統從初態5演化的研究結果是:當g=0.2時,系統開始只能產生小波,最終形成了螺旋波;當g=0.3時,系統產生的波與圖8(a)相似;當0.4≤g≤1.3時,系統出現多個螺旋波和螺旋波對;當g=1.4,1.5時系統出現單個螺旋波,如圖8(e)所示;當g≥1.6時,系統只出現圓形波,如圖8(f)所示.可見只有當神經元的初相位分布比較均勻時系統才能產生螺旋波,耦合強度過大和過小都不能產生螺旋波.

從上面的結果可以看出,單個神經元處于一周期態時,系統比較容易出現螺旋波,而且出現螺旋波的耦合強度范圍較寬.單個神經元處于二周期態時,系統只有在神經元的初相位分布比較均勻時才能產生螺旋波,因此神經元處于更高周期態,系統較難自發出現螺旋波,數值模擬也得到這樣的結論.

4 結 論

采用HR神經元模型,研究神經元陣列系統從一個具有隨機初相位分布的初態演化是否能自發形成螺旋波(包括螺旋波對),隨機的初相位分布模擬了神經元的非近鄰耦合行為,發現系統是否能自發形成螺旋波與系統初態、單個神經元的狀態、耦合強度有關,且單個神經元的狀態對系統形成螺旋波影響最大.

當單個神經元處于一周期態時,系統從任意初態演化,在一定的耦合強度范圍內,系統都能自發產生螺旋波和螺旋波對,還觀察到其他有序斑圖(如平面波、反靶波、圓形波);初態相同而耦合強度不同或者是耦合強度相同而初態不同時,系統演化一般會得到不同結果;當耦合強度超過臨界值后,系統不容易自發形成螺旋波,演化結果依賴初態.我們觀察到兩種重要現象:1)當系統初態中處于階段4的神經元占絕大多數時,系統會演化到間歇全局周期振蕩態上,甚至實現精確全局同步振蕩;2)當系統初態中處于階段5的神經元數量足夠多時,系統就會出現振蕩死亡現象.如果系統不出現間歇全局同步振蕩和振蕩死亡,在大的耦合強度下系統一般出現圓形波.

當單個神經元處于二周期態時,只有當各神經元的初相位分布比較均勻時才可以自發出現螺旋波,而且出現螺旋波的耦合強度范圍相比一周期情況有較大的縮小;當單個神經元處于更高周期態時,系統越不容易自發出現螺旋波.

上述結果有助于了解大腦皮層中螺旋波是如何自發形成的,特別是能幫助人們了解癲癇產生的機制,因為癲癇是大量神經元同步引起的一種功能性腦失調[25].本文結果能夠回答同步振蕩是如何產生的,而且大量神經元的振蕩死亡形成機制也為治療癲癇提供了有用信息.

[1]Winfree A T 1972 Science 175 634

[2]Larionova Y,Egorov O,Cabrera-Granado E,Esteban-Martin A 2005 Phys.Rev.A 72 033825

[3]Plapp B P,Egolf D A,Bodenschatz E,Pesch W 1998 Phys.Rev.Lett.81 5334

[4]Müller S C,Plesser T,Hess B 1985 Science 230 661

[5]Vanag V K,Epstein I R 2001 Science 294 835

[6]Foerster P,Müller S C,Hess B 1990 Development 109 11

[7]Davidenko J M,Pertsov A V,Salomonsz,Baxter W,Jalife J 1992 Nature 355 349

[8]Huang X,Xu W,Liang J,Takagaki K,Gao X,Wu J 2010 Neuron 68 978

[9]Huang X,Troy W C,Yang Q,Ma H,Laing C R,SchiffS T,Wu J Y 2004 J.Neurosci.24 9897

[10]Chen J X,Zhang H,Qiao L Y,Liang H,Sun W G 2018 Commun.Nonlinear Sci.Numer.Simulat.54 202

[11]Chen J X,Guo M M,Ma J 2016 Europhys.Lett.113 38004

[12]Chen J X,Liang P,Zheng Q,Zhao Y H,Ying H P 2014 Chaos 24 033103

[13]Pumir A,Nikolski V,Horning M,Isomura A,Agladze K,Yoshikawa K,Gilmour R,Bodenschatz E,Krinsky V 2007 Phys.Rev.Lett.99 208101

[14]Gan Z N,Cheng X M 2010 Chin.Phys.B 19 050514

[15]Viventi J,Kim D H,Vigeland L,Frechette E S,Blanco J A,Kim Y S 2011 Nat Neurosci.14 1599

[16]Jung P,Cornell-Bell A,Madden K S,Moss F 1998 J.Neurophysiol.79 1098

[17]García-Ojalvo J,Schimansky-Geier L 1999 Europhys.Lett.47 298

[18]Gu H G,Jia B,Li Y Y,Chen G R 2013 Physica A 392 1361

[19]Li Y Y,Zhang H M,Wei C L,Yang M H,Gu H G,Ren W 2009 Chin.Phys.Lett.26 030504

[20]Wang Q Y,Perc M,Duan Z,Chen G 2008 Phys.Lett.A 372 5681

[21]Jung P,Cornell-Bell A,Moss F,Kadar S,Wang J,Showalter K 1998 Chaos 8 567

[22]Ullner E,Zaikin A,García-Ojalvo J,Kurths J 2003 Phys.Rev.Lett.91 180601

[23]Ma J,Wu Y,Ying H P,Jia Y 2011 Chin.Sci.Bull.56 151

[24]Hindmarsh J L,Rose R M 1984 Proc.Roy.Soc.Lond.B 221 87

[25]Jiruska P,De-Curtis M,Jefferys J G R,Schevon C A,SchiffS J,Schindler K 2013 J.Physiol.591 787