頻率編碼和調(diào)制原理在Hodgkin-Huxley系統(tǒng)上的應(yīng)用

唐偉偉，　檀結(jié)慶，　唐元華

(1. 合肥工業(yè)大學(xué)數(shù)學(xué)學(xué)院，合肥230009；

2.ResearchBioinformatics，CompleteGenomics，Inc.，MountainView，California，USA.)

?

頻率編碼和調(diào)制原理在Hodgkin-Huxley系統(tǒng)上的應(yīng)用

唐偉偉1，檀結(jié)慶1，唐元華2

(1. 合肥工業(yè)大學(xué)數(shù)學(xué)學(xué)院，合肥230009；

2.ResearchBioinformatics，CompleteGenomics，Inc.，MountainView，California，USA.)

[摘要]頻率編碼指感覺(jué)神經(jīng)元在受到刺激后產(chǎn)生膜電位的周期變化，變化的頻率往往與刺激的強(qiáng)度成正比.本文在Hodgkin-Huxley系統(tǒng)的基礎(chǔ)上，對(duì)其微分方程進(jìn)行簡(jiǎn)化；對(duì)于頻率編碼，給出去耦時(shí)間尺度上的微分方程系統(tǒng)；研究在單一參數(shù)系統(tǒng)下，頻率編碼行為的變化規(guī)律.

[關(guān)鍵詞]神經(jīng)刺激； 頻率編碼；Hodgkin-Huxley系統(tǒng)； 奇異攝動(dòng)系統(tǒng)； 調(diào)頻機(jī)制

1引言

感覺(jué)神經(jīng)元可根據(jù)刺激強(qiáng)度的變化作出相應(yīng)的頻率調(diào)制.振蕩頻率通常是正比于刺激的強(qiáng)度，理解為輸入的模擬信號(hào)被頻率編碼.當(dāng)刺激強(qiáng)度越大，神經(jīng)元被激發(fā)的速率會(huì)越快；當(dāng)刺激強(qiáng)度越小，神經(jīng)元被激發(fā)的速率會(huì)越慢.當(dāng)然，這個(gè)調(diào)頻機(jī)制具有一定的線性范圍.如果一個(gè)信號(hào)太弱，可能無(wú)法偵測(cè)到它；反之，如果信號(hào)過(guò)于強(qiáng)大，神經(jīng)系統(tǒng)不堪重負(fù)，甚至可能會(huì)徹底關(guān)閉對(duì)刺激的響應(yīng).

自從首次刊登Hodgkin-Huxley具有里程碑意義的成果后[1]，可以用微分方程系統(tǒng)來(lái)呈現(xiàn)一個(gè)神經(jīng)元的調(diào)頻機(jī)制，這個(gè)方程系統(tǒng)同樣具有頻率編碼的能力.

假設(shè)對(duì)一個(gè)可興奮系統(tǒng)輸入一個(gè)模擬信號(hào)S(t)，那么輸出信號(hào)F(t)就是一個(gè)頻率編碼.隨著模擬信號(hào)的強(qiáng)度變化，輸出頻率也作出相應(yīng)的變化.因此，通過(guò)觀察輸出的頻率，我們可以猜測(cè)模擬輸入的量.模擬信號(hào)以某種方式被轉(zhuǎn)化為具有周期的頻率信號(hào).

一個(gè)神經(jīng)元系統(tǒng)可以通過(guò)一組微分方程來(lái)近似模擬，這組微分方程用來(lái)呈現(xiàn)由細(xì)胞膜表面離子通道活動(dòng)而引起的膜電位的變化[2].從數(shù)學(xué)角度，先介紹一種最簡(jiǎn)單的系統(tǒng)，它用來(lái)描述可激發(fā)介質(zhì)中的行為，這個(gè)系統(tǒng)由兩個(gè)微分方程組成[3]：

圖1　典型的可激發(fā)系統(tǒng)

當(dāng)交點(diǎn)x*在x1和x2之間時(shí)，系統(tǒng)不穩(wěn)定，并且是振蕩的；當(dāng)交點(diǎn)x*在x1左邊時(shí)，系統(tǒng)是穩(wěn)定的，并且是可激發(fā)的；當(dāng)交點(diǎn)x*在x2右邊時(shí)，系統(tǒng)是穩(wěn)定的，并且不可激發(fā).

2簡(jiǎn)化Hodgkin-Huxley微分方程系統(tǒng)

2.1　Hodgkin-Huxley方程

Hodgkin-Huxley系統(tǒng)將烏賊軸突神經(jīng)元行為賦予了數(shù)學(xué)形式.以下是Hodgkin-Huxley微分方程組系統(tǒng)的原始形式[1]：

(1)

Hodgkin-Huxley系統(tǒng)是由一個(gè)表示膜電位的方程及附加的三個(gè)表示通道門控的方程組成的.Na+離子通道是由兩個(gè)門控變量m和h來(lái)支配的，m表示激活狀態(tài)門控系數(shù)，h表示失活狀態(tài)門控系數(shù).K+離子通道僅由一個(gè)門控變量n來(lái)支配，表示激活狀態(tài)門控系數(shù).方程中沒(méi)有給出K+離子通道的失活狀態(tài)門控系數(shù).盡管如此，K+離子通道，作為一個(gè)延遲外流控制器，當(dāng)?shù)竭_(dá)一定程度時(shí)也會(huì)有失活狀態(tài).但是，在正常情況下，它對(duì)于神經(jīng)元行為來(lái)說(shuō)是多余的.

2.2　簡(jiǎn)化Hodgkin-Huxley方程

下面介紹如何將Hodgkin-Huxley系統(tǒng)簡(jiǎn)化到只含兩個(gè)微分方程的系統(tǒng).首先，Na+離子的激活過(guò)程非常迅速，通常在1ms內(nèi)完成.假設(shè)通道的激活狀態(tài)總是在預(yù)定的時(shí)間尺度上完成.可以將原系統(tǒng)的快動(dòng)作微分方程改為代數(shù)方程，即

0=am(1-m)-bmm.

其次，Na+離子失活過(guò)程和K+激活過(guò)程幾乎是在同一時(shí)間尺度上進(jìn)行的.它們的變化速率幾乎是相同的，因此用同一個(gè)門控變量來(lái)表示就足夠了.在這里，用0.5-n4來(lái)代替h，替換后，可以保證h>0，并且h總是與n步調(diào)相反.

于是得到兩個(gè)方程的Hodgkin-Huxley系統(tǒng):

(2)

方程(2)相對(duì)于原始方程(1)來(lái)說(shuō)，使Hodgkin-Huxley系統(tǒng)變得簡(jiǎn)單明了，因此更容易被理解和分析.它可以更有效的幫助我們理解神經(jīng)元細(xì)胞的復(fù)雜生理行為.

這并不是大家第一次簡(jiǎn)化Hodgkin-Huxley系統(tǒng).在這之前，F(xiàn)itzhuh-Nagoma方程[4]將原始系統(tǒng)簡(jiǎn)化成兩個(gè)微分方程系統(tǒng)，即是FitzHugh-Nagumo模型[5]：

圖2展示了簡(jiǎn)化Hodgkin-Huxley系統(tǒng)后所得方程(方程2)的nullclines.令dV/dt=0，得到一條N形曲線(圖2中的黑線)；令dn/dt=0，得到一條單調(diào)遞增曲線(圖2中的紅線).N形曲線(dV/dt=0)有一個(gè)局部極大值點(diǎn)和一個(gè)局部極小值點(diǎn).從圖2中兩條曲線的交點(diǎn)可以得出，系統(tǒng)是振蕩的.與原始Hodgkin-Huxley系統(tǒng)相比，簡(jiǎn)化后的方程顯示相同的性質(zhì).

圖2　簡(jiǎn)化后Hodgkin-Huxley系統(tǒng)(方程2)的Null clines

3頻率編碼理論

3.1　奇異攝動(dòng)系統(tǒng)

在本節(jié)中，我們研究微分方程組系統(tǒng)的頻率編碼，并給出一些初步研究結(jié)果.著重點(diǎn)是，解出下面這個(gè)微分方程組系統(tǒng)(方程3)的周期解[6].

(3)

其中x∈m，y∈n.p∈k是一個(gè)參變量，ε是一正數(shù)，并且足夠小.f和g是定義在n+m上的函數(shù).

在方程(3)中，如果f(x，y，p)是一階無(wú)窮小o(1)，相對(duì)于變量y，變量x隨時(shí)間t變化的速度更快.這表示變量x相對(duì)于變量y是快變量.總的來(lái)說(shuō)，方程組(3)中，快動(dòng)作系統(tǒng)是由第二個(gè)微分方程決定的；慢動(dòng)作系統(tǒng)是由第一個(gè)微分方程決定的.

在物理學(xué)和生物學(xué)當(dāng)中，有許多振蕩系統(tǒng)或可興奮系統(tǒng)都可以寫(xiě)成方程(3)這樣的形式.例如，Van der Pol’s 方程[7].

方程(3)中，令ε=0，得到其所對(duì)應(yīng)的退化系統(tǒng)，即方程(4)

f(x，y，p)=0，

(4)

其中x和y是標(biāo)量變量，參數(shù)p固定不變.(3)中第二個(gè)方程定義了相平面2上一曲線Γ.假設(shè)在曲線Γ上，

(fx(x，y，p))2+(fy(x，y，p))2>0，

曲線Γ上，使fx(x，y，p)=0的點(diǎn)為奇異點(diǎn)，否則為非奇異點(diǎn).假設(shè)Γ上奇異點(diǎn)是孤立的，并且fxx(x，y，p)≠0.

根據(jù)Implicit Function Theorem[8]，在非奇異點(diǎn)上，可用y=F(x，p)來(lái)表示曲線Γ.滿足fx(x，y，p)<0點(diǎn)的集合為曲線Γ的穩(wěn)定部分(stable)；滿足fx(x，y，p)>0點(diǎn)的集合為曲線Γ的不穩(wěn)定部分(unstable)；這兩部分被奇異點(diǎn)fx(x，y，p)=0隔開(kāi).以下給出方程(3)的不連續(xù)解.此方程的解由有限次慢動(dòng)作與快動(dòng)作交替組成.

慢動(dòng)作代表曲線Γ上的穩(wěn)定部分；快動(dòng)作代表曲線Γ上的不穩(wěn)定部分，它幾乎是瞬間完成的.在節(jié)點(diǎn)fx(x，y，p)=0處，慢動(dòng)作向快動(dòng)作轉(zhuǎn)化. 若在(x，y)平面上，不連續(xù)解是一極限環(huán)Z0，我們稱Z0是方程(4)的不連續(xù)周期解.不連續(xù)解的周期為

(5)

在方程(4)中，T0僅由慢動(dòng)作部分所用時(shí)間組成.快動(dòng)作部分所用時(shí)間為零.

定理1假如方程(4)有一極限環(huán)Γ0，并且滿足上述條件，那么對(duì)于一足夠小ε，其非退化系統(tǒng)，即方程(3)有唯一且穩(wěn)定極限環(huán)Zε，當(dāng)ε→0時(shí)，Zε一致收斂于Z0.另外，Zε的周期為

(6)

其中，α為一常量，o(ε)是ε的高階無(wú)窮小，當(dāng)ε→0.

參考文獻(xiàn)[8]中給出定理1的證明過(guò)程.關(guān)于Γ0周期解的實(shí)值近似，慢動(dòng)作部分是ε的無(wú)窮小，快動(dòng)作部分是εa的無(wú)窮小，其中a≥2/3.在節(jié)點(diǎn)附近，近似值為εa的無(wú)窮小，其中a≥0.然而，在節(jié)點(diǎn)附近沒(méi)有給出a的下界.這種近似方法我們稱之為零級(jí)近似[6].

定理2令εM=max(ε1，ε2)，ε1和ε2足夠小.若退化系統(tǒng)，即方程(4)有極限環(huán)Γ0，并且滿足定理1的條件，那么有

(i) 對(duì)εi，方程(3)有周期解φ(t，εi)，i=1，2；

(iii) Tε1-Tε2=o((εM)2/3).

定理1可以被擴(kuò)展到x和y都是矢量變量的情形，并且結(jié)論非常相似.但是在擴(kuò)展的過(guò)程中，要引入一些附加的條件.

3.2　單一縮放系數(shù)的慢動(dòng)作系統(tǒng)

假設(shè)m=n=1，即x∈，y∈，進(jìn)一步假設(shè)能夠從g(x，y，p)中將h提取出來(lái)，并作為一縮放系數(shù)放在前面，即g(x，y，p)).其他參數(shù)不變，并且h不出現(xiàn)在f(x，y，p)中.基于以上假設(shè)，可將方程(3)寫(xiě)成下面這種形式：

(7)

這里x和y是標(biāo)量變量，h 是參變量，并且h>0，ε1是一很小的參數(shù).

令h=h0，有

(8)

方程(8)對(duì)應(yīng)的退化系統(tǒng)為

f(x，y)=0.

(9)

假設(shè)，方程(9)有一周期解φ(t，h0，0)，根據(jù)上一節(jié)定義，其周期為T(h0，0).由定理2，知道對(duì)于一很小的值ε1，方程(8)也有一周期解φ(t，h0，ε1).此外，方程(7)的周期解φ(t，h，ε1)在軌跡和周期上非常接近方程(9)的周期解φ(t，h0，0).

(10)

定理3假設(shè)方程(9)滿足定理1 的條件，并且，r=h/h0不是太大.那么有以下結(jié)論：

(i) 對(duì)應(yīng)方程(9)的周期解φ(t，h0，0)，方程(7)有周期解φ(t，h)；

(ii) 方程(7)解φ(t，h)的周期，即T(h)與方程(8)解φ(t，h0)的周期，即T(h0)滿足

(iii)方程(7)周期解φ(t，h)與方程(8)周期解φ(t，h0)滿足

定理3說(shuō)明，隨著h的變化，方程(3.6)依然可以頻率編碼.例如，若h=2h0，有r=2.根據(jù)定理3，當(dāng)h=2h0時(shí)，解的周期為

圖3說(shuō)明在不同h值下，膜電位振蕩模式.當(dāng)h加倍，振蕩頻率加倍.

(a) h=1　　　　　　　(b)h=0.5圖3　膜電位振蕩模式

3.3　不同縮放系數(shù)的慢動(dòng)作系統(tǒng)

假如f(x，y)=0有一局部極小值x1和一局部極大值x2，并且退化系統(tǒng)周期解由四個(gè)部分組成.考慮

(11)

定理4假設(shè)方程(11)滿足正則性條件.令εM=max(εS1，εS2，εS0).對(duì)于足夠小的εM，εM>0有

(i) 系統(tǒng)(11)有解φ(t，S1，S2，ε)，周期為T(S1，S2，ε)，相對(duì)于方程(11)的周期解φ(t，S0，S0，ε)，周期為T(S0，S0，ε).另外，

(ii)

及

(12)

(iii)

及

(13)

定理3與定理4 的不同之處是，定理4 將慢動(dòng)作部分分為兩個(gè)部分.令S2=S0，通過(guò)改變S1，S1>1，可以在不改變波峰寬度的情況下，改變恢復(fù)時(shí)間，從而達(dá)到改變振蕩頻率的目的；或者令S1=S0，通過(guò)改變S2，S2>1，在不改變恢復(fù)時(shí)間的情況下，改變波峰寬度，同樣達(dá)到改變振蕩頻率的效果.如果令S1=S2改變縮放系數(shù)r.波峰寬度及恢復(fù)時(shí)間同時(shí)改變.

圖4說(shuō)明在不同S1，S2值下，膜電位的振蕩模式.

(a) S1=S2=1　　　(b) S1=0.25, S2=1　　　(c) S=1, S2=0.25圖4　膜電位振蕩模式

最后，解決調(diào)頻機(jī)制的問(wèn)題.如果有一隨時(shí)間變化的信號(hào)，將此信號(hào)引入系統(tǒng)，那么系統(tǒng)的頻率是如何變化的，如圖5所示.反之，如果觀察到系統(tǒng)中頻率發(fā)生改變，能推導(dǎo)相對(duì)應(yīng)輸入的模擬信號(hào)嗎？ 這個(gè)問(wèn)題非常復(fù)雜，并牽扯到多時(shí)間尺度.這里給出一個(gè)最簡(jiǎn)單結(jié)論.假設(shè)輸入的模擬信號(hào)是分段函數(shù)，在這一假設(shè)條件下，可衍生出一公式，這個(gè)公式反應(yīng)輸入信號(hào)與輸出頻率的關(guān)系.

圖5　膜電位振蕩模式S1(t)=1+sin(2π(t/1000)-π/2)

4結(jié)論與認(rèn)識(shí)

將系統(tǒng)從二維平面推廣到更高維空間.

(14)

Hodgkin-Huxley系統(tǒng)就是一個(gè)案例.在HH系統(tǒng)中，改變某個(gè)離子通道的狀態(tài)只能一次調(diào)整系統(tǒng)振蕩頻率.目前能夠詮釋這種系統(tǒng)的最佳方法是仿真模擬.關(guān)于研究心臟竇房結(jié)(SA節(jié)點(diǎn))機(jī)制，我們已經(jīng)找到一個(gè)非常好的模型[9]，這個(gè)模型反應(yīng)電信號(hào)及振蕩行為.但是， 神經(jīng)遞質(zhì)是如何影響心率調(diào)節(jié)的呢？我們知道在恢復(fù)階段，有多重電流被激活，例如A-current[10].但是在細(xì)胞中， A-current 是如何調(diào)整頻率編碼的[11].這些問(wèn)題有待于研究.

[1]HodgkinAL，HuxleyAF.Aquantitativedescriptionofmembranecurrentanditsapplicationtoconductionandexcitationinnerve[J].TheJournalofphysiology， 1952， 117(4): 500-544.

[2]HodgkinAL，HuxleyAF.CurrentscarriedbysodiumandpotassiumionsthroughthemembraneofthegiantaxonofLoligo[J].TheJournalofphysiology, 1952， 116(4): 449-472.

[3]TangY，OthmerHG.Frequencyencodinginexcitablesystemswithapplicationstocalciumoscillations[J].Proc.NationalAcademySciences, 1995， 92(17): 7869-7873.

[4]FitzHughR.Impulsesandphysiologicalstatesintheoreticalmodelsofnervemembrane[J].BiophysicalJournal, 1961， 1(6): 445-466.

[5]NagumoJ，ArimotoS，YoshizawaS.Anactivepulsetransmissionlinesimulatingnerveaxon[J].ProcIRE, 1962， 50(10):2061-2070.

[6]TangY.TheoreticalStudiesonSecondMessengerDynamics[M].AnnArboor,Michigan:UniversityMicrofilms,ProQuestLIC., 1993.

[7]VanderPolB.Onrelaxation-oscillations[J].TheLondon，EdinburghandDublinPhil.Science1927， 2(7): 978-992.

[8]MishchenkoEF，RosovNK.DifferentialEquationswithSmallParametersandRelaxationOscillations[M].NewYork:PlenumPress, 1980.

[9]NobleD.Modellingtheheartfromgenestocellstothewholeorgan[J].Science, 2002， 295(5560): 1678-1682.

[10]RushME，RinzelJ.ThepotassiumA-current，lowfiringratesandreboundexcitationinHodgkin-Huxleymodels[J].BulletinofMathematicalBiology, 1995， 57(6): 899-929.

[11]ConnerJA，StevensCF.Predictionofrepetitivefiringbehaviourfromvoltageclampdataonanisolatedneuronesoma[J].TheJournalofPhysiology, 1971， 213(1): 31-53.

TheoryofFrequencyEncodingandModulationwith

ApplicationtotheHodgkin-HuxleySystem

TANG Wei-wei1，TAN Jie-qing1，TANGY.Tom2

(1.DepartmentofMathematics，HefeiUniversityofTechnology，Hefei230009，China;

2.ResearchBioinformatics，CompleteGenomics，Inc.，MountainView，California，USA.)

Abstract:Sensoryneuronsrespondtostimulatorysignalswithperiodicfiringsofthemembranepotential.Thefrequencyoffiringisoftenproportionaltotheintensityofstimulus，andonesaysthatthestimulusisfrequencyencoded.HerebasedontheoriginalHodgkin-Huxleysystem，wesimplifyitsdifferentialequations;forfrequencyencoding，developsystemsdescribedbywithdecouplingtimescales;studydifferentmodesoffrequencyencodingariseasasingleparameterofthesystemisvaried.

Keywords:nervestimulation;frequencyencoding;Hodgkin-Huxleysystem;singularlyperturbedsystems;frequencymodulation

[收稿日期]2014-09-05

[中圖分類號(hào)]O193

[文獻(xiàn)標(biāo)識(shí)碼]A

[文章編號(hào)]1672-1454(2015)01-0014-07