Ghostburster模型鞍結分岔點附近放電模式的轉遷控制

陳穎源，王 江，曾啟明，陳偉樂，魏熙樂，鄧 斌

(1. 天津大學電氣與自動化工程學院，天津 300072；2. 香港理工大學機電工程學院，香港 999077)

在神經系統中，神經元的放電模式蘊含著大量的信息．聲、光、電以及機械等外界刺激通過影響神經元的生理參數引發不同的放電模式，如峰放電、簇放電等．峰放電在動力學上對應著穩定的極限環[1-2]，而簇放電的動力學特性更加復雜[2-9]．受快慢離子通道的調節，神經元可表現出混沌簇放電[2-4]、反彈簇放電[5-6]等多種簇放電形式．研究認為簇放電在從受體神經元到腦皮層的神經信息處理中發揮了重要的作用[7-8]，而且還可以加強腦皮層各區域之間的聯系[9]．同時，不正常的簇放電則會破壞神經系統的信息傳遞．比如帕金森患者發病時在丘腦中繼細胞上觀察到的大量的反彈簇放電[5-6]．

神經元隨著參數的變化會產生分岔，從而引起不同放電模式的轉換．所以，通過調整分岔點的位置，可以使神經元在某些參數下達到期望的放電特性[1,10]．文獻[11]利用卡爾曼濾波器控制小腦皮層模型的頻率和放電模式．文獻[12]利用狀態反饋控制來抑制心臟中不正常的動作電位．文獻[13]利用washout濾波器控制Hopf分岔點來實現HH模型靜息狀態和周期放電狀態之間的轉換．washout濾波器不改變原系統的平衡點，其在神經系統 Hopf分岔控制[13-16]和混沌控制[16-17]中應用較廣．但是，對神經系統進行鞍結分岔控制的研究鮮見文獻報道．筆者從Ghostburster模型[18-19]極限環上的鞍結分岔點出發，通過washout濾波器調整鞍結分岔點的位置，以實現鞍結分岔點兩端放電模式的轉遷．

Ghostburster模型是來自于弱電魚電敏感側線葉(electrosensory lateral line lobe，ELL)椎體細胞的二室模型，有豐富的放電模式[18-19]．2000年，Lemon和Turner[20]分析了 ELL椎體細胞中頻率的簇放電機制．2002年，Doiron等[18-19]研究了 ELL椎體細胞在外電場作用下的非線性特性，把多室結構降維為二室并命名為 ghostburster模型．他們提出該模型的簇放電是由快子系統固定點詭異(ghost)的鞍結分岔產生的，這種鞍結分岔有別于傳統簇放電發生時快子系統從靜息態到極限環的分岔類型[1]，靜息態并不存在，簇放電的時間間隔主要由無限周期分岔決定，這種分岔形式稱為“ghost”．這種特殊的簇放電為了解放電模式提供了新的見解．

筆者在 washout濾波器產生的控制效應的基礎上，通過快慢系統分解，分析了 washout濾波器作用到鞍結分岔點上的內在機制，發現快子系統固定點的鞍結分岔和快子系統從周期一極限環轉換到周期二極限環的臨界點之間的相互作用是該模型周期放電和簇放電產生的關鍵因素．

1 Ghostburster模型的動力學分析

Ghostburster模型[18-19]是弱電魚 ELL椎體細胞的二室模型．它描述了胞體和樹突之間信息的傳導，從整體上刻畫了柱狀細胞的動力學行為．Ghostburster模型包含了6個非線性微分方程．

1)胞體部分

2)樹突部分

其中，式(1)和(3)分別為胞體和樹突膜電壓的微分方程，Iu為控制電流項．胞體和樹突之間通過簡單的耗散電流耦合．各個離子通道電流(INa,s，IDr,s，INa,d，IDr,d)受最大電導 gmax(單位為 m s/cm2)，S型激活和失活概率曲線w = n ,h,p)及通道時間常數τ(單位為ms)影響，其參數選擇參見文獻[18]．其他參數選擇如下：gc=1，κ= 0 .4，VNa= 4 0mV，VK=?8 8.5mV ，Vleak=?7 0mV，gleak= 0 .18，Cm= 1 μF/cm2．本文用去極化電流 Is作為分岔參數來研究 Ghostburster模型的動力學特性．圖1是放電頻率f( f = 1 /TISI，ISI為峰峰間期)和最大Lyapunov指數λ隨著分岔參數 Is變化的圖形．

圖1 放電頻率f和Lyapunov指數λ隨去極化電流Is變化的分岔圖Fig1 Bifurcatim diagram of f and λ as Is varied

(1) 當Is＜Is1時，神經元處于靜息狀態，系統存在一個穩定的平衡點和一個不穩定的平衡點，這時f=0，λ=0．

(2) 當Is= Is1( ≈ 5 .8μ A /cm2)時產生了不變環上固定點的鞍結分岔，2個平衡點碰撞消失，產生一對穩定的極限環和不穩定的極限環．

(3) 當 Is1＜Is＜Is2時，穩定的極限環和不穩定的極限環共存．神經元呈現周期放電狀態，這時峰峰間期(ISI)為常數，f也為常數，λ=0．在該區間內的典型響應如圖2(a)和(b)所示．

(4) 當Is= Is2( ≈ 8 .5μ A /cm2)時產生極限環的鞍結分岔，這是一種接近余維2(near codimension-2)的鞍結分岔[1]．穩定的極限環和不穩定的極限環碰撞消失，從而產生混沌現象．

(5) 當 Is2＜ Is＜ Is3時，神經元不存在穩定的極限環．在大部分參數下其 ISI是變化的，故f也是變化的．λ＞0表明神經元處于混沌狀態，而λ=0表明系統進入了周期窗口(比如 Is在 [ 13.13μA/cm2，13.73μA /cm2]區間內存在一個六周期解)．在該區間內的典型響應如圖2(c)和(d)所示．

(6) 當 Is= Is3≈ 1 7.65μA /cm2時，系統在該參數附近產生倍周期分岔，產生了二周期解．

(7) Is＞Is3時神經元呈二周期簇放電狀態，所以λ= 0 ．在該區間內的典型響應如圖2(e)和(f)所示．

圖2 Ghostburster模型在3類典型電流輸入下的響應Fig.2 Time series of somatic membrane potential Vs and dendritic membrane potential Vd of ghostburster model with 3 different current stimulations

2 Ghostburster模型的 washout濾波器反饋控制

2.1 Washout濾波器

Washout濾波器實際上是一個高通濾波器，其穩態輸入被濾除，而暫態輸入則得到保留．一維washout濾波器[13,17,21]可設計為

式中：x為 washout濾波器的輸入；d為 washout濾波器時間常數的倒數；z為 washout濾波器的狀態；y為輸出函數． ()gy是待設計的控制函數．濾波器有2個條件需要滿足，一個是 d ＞ 0 以保證濾波器的穩定性．另一個是 g ( 0)= 0 以保留初始平衡點．加入washout濾波器后，新系統(式(1)～式(7))將變成七階系統．因為 Vs在反饋系統中易于測量，令 x = Vs作為唯一的控制器輸入．為了簡化控制參數的選擇，控制函數只包含一個線性項，表示為

式中K為反饋控制器增益．由圖 1可知 Is=Is2(≈8.5μA/cm2)2端存在著2種不同的狀態：周期放電和混沌簇放電．故本文把控制目標設定為實現這2種狀態之間的相互轉換．根據 Izhikevich的描述，Is2附近發生的是極限環上接近余維 2的鞍結分岔[1]，要從數值上精確確定分岔點的位置非常困難．所以，只能從經驗上通過調節參數K和d來達到期望的特性．

2.2 Washout濾波器控制結果

2.2.1 從混沌簇放電到周期放電的轉遷

圖3為washout濾波器參數的變化對最大Lyapunov指數λ的影響．由圖 3可知在同一d下隨著K增加λ呈遞減趨勢，但是并不單調．而且，d不同時λ作為K的函數變化并無太大區別．K增加到一定值后λ為負值，說明神經元呈現周期放電狀態．為了分析方便，把濾波器參數設置為 K = 0 .3、d = 0 .4.圖 4表示的是加上控制器前后胞體和樹突膜電壓的時間序列．控制量在 t = 2 00ms 時加入．在經過一個短的暫態過程后，神經元從混沌簇放電轉換為周期放電．圖 5是加入控制器以后，λ和f隨著 Is變化．因為washout濾波器并不改變平衡點的特性，第1分岔點仍然在 Is1= 5 .8μA /cm2處．第 2分岔點 Is2從8.5μA/cm2移至11.7μA/cm2，使I=9.0在區間s[Is1, Is2]內，模型呈周期放電狀態．第 3分岔點 Is3從17.65μA/cm2延后至20.8μA/cm2．混沌區域[I ,I ]的s2s3大小幾乎保持不變，只是做了一段平移．這是因為washout濾波器只是改變了全系統的極限環的鞍結分岔位置，而快子系統固定點的鞍結分岔幾乎沒有任何改變造成的．

圖3 在不同d下λ隨K變化的情況(Is=9.0 μA/cm2)Fig.3 Maximum Lyapunov exponents λ as a function of K while parameter d varied(Is=9.0 μA/cm2)

圖4 Is =9.0 μA /cm2，K =0.3和 d =0.4時神經元的放電時間序列(控制項在t=200 ms時加入)Fig.4 Time series of ghostburster model with I s =9.0 μA /cm2，K=0.3 and d=0.4 (t=200 ms)

圖 5 加入控制器(參數為 0.3K= ， 0.4d= )后 ghostburster模型狀態隨著去極化電流 sI變化的分岔圖Fig.5 Bifurcation diagram of fand λafter adding controller with 0.3=K ，0.4=d

2.2.2 從周期放電到混沌簇放電或其他放電模式的轉遷

在極限環的鞍結分岔點左方取 Is= 8 .0μA /cm2．為分析方便仍然取 d = 0 .4．圖 6給出了在 K =?0 .2、K=?1 .0和 K =?2 .0時胞體的動作電位．控制項在t= 1 00ms時加入．隨著K減小，神經元先是轉變成周期簇放電，然后再轉變成高頻峰放電．在混沌區域內，λ劇烈波動．通過λ的變化確定系統混沌時K的取值范圍．當d=0.4時，K的取值在[?0.8,0]之間．隨著d的增加，K的取值范圍也在增加．圖7(b)圖．因為 pd對快子系統的直接作用對象是樹突膜電壓 Vd，為分析方便選擇 Vd作為快子系統的代表變量進行研究．分岔圖頂部帶三角形標志實線為在不同pd下 Vd在周期軌道上的極大值． pd= pd1是快子系統將 -Kd平面( 0K＜ )劃分成2部分，神經元在該平面上半部分呈混沌簇放電狀態，在下半部分呈周期簇放電狀態．所以，K為負時可以達到將周期放電轉換為混沌簇放電或高頻放電的目的．

圖 6 不同的 washout控制參數下胞體的放電時間序列(Is＝8.0 μA/cm2，d＝0.4，控制項在 t＝100 ms時加入)Fig.6 Membrane potential of soma before and after control with different parameters (Is＝8.0,μA/cm2，d＝0.4，and t＝100 ms)

圖7 不同控制參數下放電模式的區域劃分Fig.7 Regions of firing patterns divided by different control paramters

2.3 Washout濾波器對 ghostburster模型的作用機制分析

為了分析washout濾波器對ghostburster模型分岔點的作用機制，利用奇異攝動理論[22-23]，將ghostburster模型分成一快一慢2個子系統，表示為

式中x為快變量的向量．在原始系統中，該向量包含了 Vs、Vd、ns、hd和 nd5個變量．因為 washout濾波器的常數大約為2.5 ms(1d)，約為τpd的一半，所以washout濾波器變量z也應為快變量．在快子系統變化過程中 pd可以近似為常數，所以其為慢變量． f ( x , pd)在原始模型中表示的是式(1)～式(5)，而在控制模型中表示的是式(1)～式(5)和式(7)．

2.3.1 從混沌簇放電到周期放電轉遷的情況

圖 8 對于延后鞍結分岔點的快慢系統分析( Is= 9 .0μA /cm2，K= 0 .3， d = 0 .4)Fig.8 Analysis of fast-slow subsystem in the case of delaying bifurcation point( I s = 9 .0 μ A /cm2,K = 0 .3,d = 0 .4)

圖 8(a)是當 Is= 9 .0μA /cm2時快子系統的分岔周期一極限環和周期二極限環的分界點．在 pd＞pd1時 Vd只有 1個極大值，而在 pd＜pd1時有2個極大值．分岔圖中部帶方塊標志實線為 2個連續放電之間 Vd的平均值 Vd，表示為式中 ti是第i個放電的放電時間．樹突膜電壓 Vd在pd＜pd1時呈二周期簇放電狀態，其第2個放電的峰值減小，和下一個放電之間的時間間隔拉長．所以Vd在pd= pd1處也出現了一個明顯的拐點．分岔圖底部帶空心圓點的長虛線 d pdd t = 0 為 pd零值線，其和 Vd曲線交于 pd= pd2＜ pd1，在該點處由正變負．在 pd=(= 2 .207749＞ pd1)處快子系統出現了固定點的鞍結分岔．因為 pd是鉀離子通道的失活變量，增加 pd相當于減小去極化電流 Is．

把全系統的解投射到 Vd-pd相平面上(圖 8(b)中帶箭頭的實線)來分析簇放電的產生過程．假設解軌跡從位于 pd零值線下方的 A 點出發，因為dpddt ＞ 0 ，軌跡向右邊緩慢移動．在穿過 pd零值線后，由于 d pdd t ＜ 0 ，軌跡向左移動并突然產生一次放電上跳至 Vd最大值曲線的某一點．之后，受復極化變量 hd的作用，Vd下落跌回 pd零值線的下方．Vd-pd軌跡盤旋左移遠離固定點的鞍結分岔點 pd=，放電頻率也逐漸增加．在進入 pd(t)＜pd1區域以后，軌跡從周期一峰放電轉換成周期二峰放電，這是一次高頻二次放電，第 2次放電的幅值會減小從而使 Vd減小．當，p (t)重新增加逐漸d進入“詭異的”鞍結分岔點． pd越接近，移動也就越緩慢，所以簇放電之間的間隔拉長．這種簇放電是混沌的，具有不可重復性，其產生取決于解軌跡能否從周期一區域進入周期二區域從而產生不完全的放電[18]．

圖8(c)是加上washout濾波器反饋控制(參數為d= 0 .4， K = 0 .3)以后快子系統的分岔圖．當dpddt = 0 時，pd由 Vd唯一確定，所以 pd零值線和圖8(a)相比沒有變化．通過對比還可以發現 pd1左移而pd2右移，它們相互靠近使 Vd-pd軌跡進入周期二區域越來越難． pd1和 pd2越靠近，產生不完全的放電就越困難．解軌跡向右移動，而固定點的鞍結分岔并沒有變化． pd在進入 pd(t)＜pd1區域之前返回高值區域，在這種情況下產生簇放電的條件并不滿足．所以系統只能產生周期放電(圖 4(a)和(b)，t = 2 00ms 后部分)．

2.3.2 從周期放電到混沌簇放電或其他放電模式的情況

圖 9(a)是當 Is= 8 .0μA /cm2時，原系統的快子系統分岔圖及全系統解在 Vd-pd平面上的投影疊加到分岔圖上的情形，pd1和 pd2的距離和圖 8(a)相比縮短了，而固定點的鞍結分岔= 1 .472756則更接近pd1．在 pd1，pd2和相互作用下，解軌跡被夾在 pd1和之間，系統周期放電．

圖 9(b)為當 washout濾波器參數為 d = 0 .4，K=? 0 .2時快子系統的分岔圖．和圖 9(a)相比，固定點的鞍結分岔保持不變，但是 pd1的位置右移而pd2左移，兩者距離拉長使解軌跡易于進入周期二區域，所以，整個軌跡向左移動，系統呈混沌放電狀態(圖6(a)，t = 1 00ms后部分)．繼續減小washout濾波器的參數K(圖 9(c)，K =?1 .0)，pd1和解軌跡都繼續左移，而 Vd曲線上升，和 pd零值線不再具有交點．全系統只能產生一個完全放電和一個不完全放電，呈周期二簇放電狀態(圖 6(b)，t = 1 00ms 后部分)．進一步減小K( K=?2 .0)時(圖 9(d)～(f))，pd1和之間的距離繼續減小，周期二區域右移．在圖 9(d)中 Vd,max左端，快子系統從右到左經歷了一次從周期二極限環到周期一極限環的轉換過程，其臨界點記為 pd2．為簡單起見，把加入控制器以前，加入控制器以后的過渡階段和穩定階段的解軌跡分別疊加到分岔圖上．從圖 9(d)可知軌跡是周期的，但是在washout濾波器加入以后，解軌跡馬上遠離固定點的鞍結分岔．周期二區域同時右移把解軌跡包含在內．因為解軌跡不能馬上變化，它只能盤旋左移產生不完全的放電并穿越 pd2圖9(e)．在 pd＜pd2時，解軌跡逐漸收斂成一個極限環(圖 9(f))．所以，樹突膜電壓呈高頻小幅震蕩狀態，而在胞體則體現為高頻峰放電．

綜上所述，外電流的加入改變了 Ghostburster模型快子系統從周期一極限環轉換到周期二極限環的臨界點 pd1，樹突平均模電壓的正負分界點 pd2和快子系統固定點的鞍結分岔的位置．而 washout濾波器只改變 pd1和 pd2的位置．考慮到 pd2并不總是存在(如圖9(c) 和(d)的情況)，可以得出washout濾波器通過影響 pd1產生樹突的不完全放電，從而影響整個ghostburster神經元是否產生簇放電的結論．

圖9 對于提前鞍結分岔點的快慢系統分析( Is=8.0μA /cm2，d=0.4)Fig.9 Analysis of fast-slow subsystem in case of advancing bifurcation point( I s=8.0μA /cm2，d=0.4)

3 結 語

將 washout濾波器作為反饋控制器加入到ghostburster模型中以實現極限環上的鞍結分岔點附近的放電模式轉換．然后，通過把系統分成快慢 2個子系統，研究了 washout濾波器作用于 ghostburster模型的機制．因為washout濾波器并不改變系統平衡點(包括全系統和快子系統)的位置，所以其通過改變從周期一極限環轉換到周期二極限環的臨界點的位置來實現放電模式的轉遷．而外加電流除了改變上述分岔外，還可以改變快子系統固定點的鞍結分岔．本文的研究明確了 ghostburster模型的簇放電形成機制．在以后的研究中會繼續探討鞍結分岔控制對神經系統的作用機制，以期為神經系統疾病診斷和治療提供參考．

［1］ Izhikevich E M. Dynamical Systems in Neuroscience：The Geometry of Excitability and Bursting[M]. Cambridge：The MIT Press，2007.

［2］ 楊卓琴，陸啟韶. 神經元 Chay模型中不同類型的簇放電模式[J]. 中國科學，2007，37(4)：440-450.Yang Zhuoqin，Lu Qishao. Different bursting patterns in chay neuron model[J]. Science in China，2007，37(4)：440-450(in Chinese).

［3］ 楊卓琴. 神經元電活動不同節律模式的幾種變化過程[J]. 物理學報，2010，59(8)：5319-5324.Yang Zhuoqin. Transitions of different rhythm patterns of neuronal firing activity[J]. Acta Phys Sin，2010，59(8)：5319-5324(in Chinese).

［4］ 王寶燕，徐 偉，邢真慈. 外界電場激勵下的耦合FitzHugh-Nagumo神經元系統的放電節律研究[J]. 物理學報，2009，58(9)：6590-6595.Wang Baoyan，Xu Wei，Xing Zhenci. Fire patterns of coupled FitzHugh-Nagumo neurons exposed to external electric field[J]. Acta Phys Sin，2009，58(9)：6590-6595(in Chinese).

［5］ Rubin J E，Terman D. High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model[J]. J Comput Neurosci，2004，16(3)：211-235.

［6］ Terman D，Rubin J E，Yew A C，et al. Activity patterns in a model for the subthalamopallidal network of the basal ganglia[J]. J Neurosci，2002，22(7)：2963-2976.

［7］ Singer W，Gray C M. Visual feature integration and the temporal correlation Hypothesis[J]. Annu Rev Neurosci，1995，18：555-586.

［8］ Wehr M，Laurent G. Odour encoding by temporal sequences of firing in oscillating neural assemblies[J]. Nature，1996，384(14)：162-166.

［9］ Singer W. Synchronization of cortical activity and its putative role in information processing and learning[J].Annu Rev Neurosci，1993，55(3)：349-374.

［10］ Chen G R，Moiola J L，Wang H O. Bifurcation control：Theories，methods，and applications[J]. Int J Bif and Chaos，2000，10(3)：511-548.

［11］ Schiff S J，Sauer T. Kalman filter control of a model of spatiotemporal cortical dynamics[J]. J Neural Eng，2008，5(1)：1-8.

［12］ Munoz L M，Stockton J F，Otani N F. Applications of control theory to the dynamics and propagation of cardiac action potentials[J]. Annals of Biomedical Engineering，2010，38(9)：2865-2876.

［13］ Xie Y，Chen L N，Kang Y M，et al. Controlling the onset of Hopf bifurcation in the Hodgkin-Huxley model[J]. Physical Review E，2008，77(6)：061921.

［14］ Wang J，Chen L Q，Fei X Y. Bifurcation control of the Hodgkin-Huxley equations[J]. Chaos，Solitons & Fractals，2007，33(1)：217-224.

［15］ 羅 明，彭建華，吳勇軍. Washout-filter輔助的FitzHugh-Nagumo神經元放電行為的分岔控制[J]. 動力學與控制學報，2009，7(3)：226-229.Luo Ming，Peng Jianhua，Wu Yongjun. Bifurcation control of FitzHugh-Nagumo neuronal firing behabior aided with Washout-filter[J]. Journal of Dynamics and Contrrol，2009，7(3)：226-229(in Chinese).

［16］ Wang H O，Abed E H. Bifurcation control of a chaotic system[J]. Automatica，1995，31(9)：1213-1326.

［17］ Zhou S B，Liao L F，Wu Z F，et al. Washout filterbased chaotic control in a neural network with uncertain time delays [C]//Proceedings of International Conference on Communications，Circuits and Systems. Hong Kong，China，2005：959-963.

［18］ Doiron B，Laing C R，Longtin A. Ghostbursting：A novel neuronal burst mechanism[J]. J Comput Neurosci，2002，12(1)：5-25.

［19］ Laing C R， Doiron B， Longtin A， et al.Ghostbursting：The effects of dendrites on spike patterns[J]. Neurocomp，2002，44(3)：127-132.

［20］ Lemon N，Turner R W. Conditional spike backpropagation generates burst discharge in a sensory neuron[J].J Neurophysiol，2000，84(3)：1519-1530.

［21］ 趙益波，羅曉曙. 基于washout濾波器技術的Colpitts振蕩器混沌控制研究[J]. 物理學報，2007，56(11)：6258-6262.Zhao Yibo，Luo Xiaoshu. Study on chaotic control of Colpitts oscillator based on washout filter technique[J].Acta Phys Sin，2007，56(11)：6258-6262(in Chinese).

［22］ Guckenheimer J. Progress in Nonlinear Differential Equations and Their Applications[M]. Boston：Birkhauser，1996.

［23］ Rinzel J，Lee Y S. Nonlinear Oscillations in Biology and Chemistry[M]. New York：Springer，1986.