任玉強 柴新禹 任秋實 李麗明*

1(上海交通大學(xué)生物醫(yī)學(xué)工程學(xué)院，上海 200240)

2(北京大學(xué)工學(xué)院生物醫(yī)學(xué)工程系，北京 100871)

引言

在現(xiàn)代神經(jīng)科學(xué)研究中，探索神經(jīng)編碼、信息傳遞及處理的生理學(xué)機理越來越成為科學(xué)工作者關(guān)注的焦點。通過數(shù)學(xué)建模及計算機仿真對神經(jīng)元及神經(jīng)元網(wǎng)絡(luò)進行模擬并揭示其信息編碼及傳遞機理已經(jīng)成為神經(jīng)科學(xué)家進行科學(xué)研究的重要手段。一種具體的實驗手段往往只能研究單一層次的問題，而要把許多層次上得到的大量實驗數(shù)據(jù)聯(lián)系和組織起來并得出定量的規(guī)律，對神經(jīng)系統(tǒng)進行建模及仿真的研究手段發(fā)揮了重要的作用。目前，對神經(jīng)元及神經(jīng)元網(wǎng)絡(luò)進行建模的方法有多種，其中，人工神經(jīng)網(wǎng)絡(luò)是一種應(yīng)用類似于大腦神經(jīng)元突觸聯(lián)接的結(jié)構(gòu)進行信息處理的數(shù)學(xué)模型，其目的不在于闡明神經(jīng)信息處理的機制，而在于構(gòu)建出具有某種特定功能的神經(jīng)網(wǎng)絡(luò)模型，以便解決工程上的問題，該內(nèi)容不在本文的綜述范圍之內(nèi)。本文主要對基于神經(jīng)元生物物理參數(shù)進行建模的電生理模型(又稱為生物物理模型)的研究方法及進展進行介紹。

神經(jīng)元電生理模型是指從神經(jīng)元膜的生物物理性質(zhì)出發(fā)，將神經(jīng)元分成不同的部分，不同部分彼此耦合起來構(gòu)成整個神經(jīng)元及神經(jīng)元網(wǎng)絡(luò)的建模方法，其目的是用來闡明神經(jīng)元及神經(jīng)元網(wǎng)絡(luò)的電生理學(xué)特性以及信息編碼機理。神經(jīng)元電生理模型可以在多個層次上進行建立，以實驗研究結(jié)果為基礎(chǔ)，采用精細(xì)的生物物理參數(shù)對單個神經(jīng)元進行建模，可以研究單個離子通道及受體的功能信息以及單個神經(jīng)元的興奮特性及動力學(xué)機理，而在多神經(jīng)元網(wǎng)絡(luò)的水平上建立的模型，能對神經(jīng)元間的功能聯(lián)系及神經(jīng)信息在神經(jīng)網(wǎng)絡(luò)中的傳遞機制進行研究。例如，不同神經(jīng)元通常包含了各種不同的離子通道和膜受體，不同種類的離子通道在激活電壓、失活電壓、電導(dǎo)率等電生理學(xué)特性上往往具有顯著差異，這種差異性表明不同的離子通道在神經(jīng)信息的傳遞中發(fā)揮不同的作用，通過對神經(jīng)元進行電生理建模，就能研究不同種類離子通道在興奮的產(chǎn)生及傳導(dǎo)中所起的作用。

文中將主要從單神經(jīng)元到局部神經(jīng)元網(wǎng)絡(luò)的建模方法、不同模型在神經(jīng)科學(xué)中的應(yīng)用研究進展以及各種仿真工具等方面進行介紹。

1 建模方法

1.1 霍奇金-赫胥黎模型

現(xiàn)代神經(jīng)元電生理模型的核心數(shù)學(xué)框架源自于霍奇金和赫胥黎所建立的模型[1]。通過一系列設(shè)計精妙的實驗，霍奇金和赫胥黎發(fā)現(xiàn)烏賊巨大軸突中離子電流的產(chǎn)生可以用軸突膜中鈉離子和鉀離子通道電導(dǎo)率的變化來解釋，并基于數(shù)理科學(xué)的知識和方法，建立了鈉離子和鉀離子通道電導(dǎo)率隨膜電位和時間變化的數(shù)學(xué)模型，進而推導(dǎo)出了導(dǎo)致動作電位產(chǎn)生的一系列微分方程，這就是霍奇金-赫胥黎模型(Hodgkin-Huxley模型、H-H模型)。該模型與實驗結(jié)果具有高度的一致性，為闡明動作電位是如何產(chǎn)生和傳導(dǎo)的這一神經(jīng)科學(xué)基本問題奠定了理論基礎(chǔ)，霍奇金和赫胥黎也因此獲得了1963年的諾貝爾生理學(xué)或醫(yī)學(xué)獎。

1.1.1 等效電路

在基于生物物理學(xué)的神經(jīng)元建模中，神經(jīng)元的電學(xué)特性是通過等效電路來模擬的。用電容來模擬細(xì)胞膜的電荷存儲能力，用電阻來模擬膜中不同的離子通道，而膜內(nèi)外離子濃度差異產(chǎn)生的電化學(xué)電位用電動勢來描述。一小段烏賊巨大軸突的等效電路示意圖如圖1所示，在等效電路中跨膜電流分為兩部分，一部分是通過膜電容的電流，一部分是通過離子通道的電流，在H-H模型中離子電流分為三種：鈉電流INa、鉀電流IK、主要由氯離子產(chǎn)生的漏電流IL。通過基爾霍夫電流定律可以將圖1中的等效電路表示為微分方程：

式中，Iext是外部施加的電流，Vm是跨膜電位(Vin-Vout)，該式表示了跨膜電位改變和離子電流之間的關(guān)系。

圖1 一小段烏賊巨大軸突的等效電路示意圖[1]Fig.1 Electrical equivalent circuit for a short segment of squid giant axon[1]

圖1中的離子電流代表大量離子通道形成的宏觀電流，在模型中宏觀電流遵循歐姆定律 I=GV，G代表電導(dǎo)。將此公式代入式(1)中，V就是跨膜電位Vm和平衡電位Ei(i代表不同的離子通道)的差值，整個離子電流就是不同離子通道電流的代數(shù)和：

式中，GNa、GK、GL分別代表鈉通道、鉀通道、漏電流通道的電導(dǎo)，ENa、EK、EL分別代表鈉通道、鉀通道、漏電流通道的平衡電位。通常，電導(dǎo)都不是常數(shù)，而是依賴于膜電位并隨時間發(fā)生變化。

1.1.2 門控通道

H-H模型中的宏觀電導(dǎo)可以看作是鑲嵌在膜上眾多離子通道共同作用的結(jié)果，每個離子通道可以理解為由一個或幾個電壓門所調(diào)控，每個電壓門都有兩種狀態(tài)：自由態(tài)和非自由態(tài)。當(dāng)一個離子通道的所有電壓門都處于自由態(tài)時，離子通道打開，離子就能通過離子通道形成離子電流，如果任意一個電壓門處于非自由態(tài)，離子通道就關(guān)閉。

考慮單獨的一個電壓門 i，定義pi為這個門處于自由態(tài)的概率(0～1)。如果考慮大量的電壓門，那么pi可以看作是所有的類型為i的電壓門中處于自由態(tài)的比例，pi(t)就代表在時刻t時處于自由態(tài)的電壓門比例，1-pi(t)代表處于非自由態(tài)的電壓門比例。

電壓門從非自由態(tài)轉(zhuǎn)化到自由態(tài)的速率為αi(V)，從自由態(tài)轉(zhuǎn)化到非自由態(tài)的速率為 βi(V)，兩種速率都依賴于跨膜電壓。在H-H模型中，假定電壓門在兩種狀態(tài)間的改變遵循一階動力學(xué)方程：

當(dāng)一個通道開放時，電導(dǎo)會有一定程度的上升，因此大量離子通道形成的宏觀電導(dǎo)正比于開放通道的數(shù)量，也就正比于相關(guān)電壓門處于自由態(tài)的概率。因此離子通道j的宏觀電導(dǎo)Gj正比于組成通道j的所有電壓門處于自由態(tài)比例的乘積：

式(8)～式(11)以及式(1)合在一起就能夠完整描述H-H模型中跨膜電位隨時間的變化。為了完整的求解這個模型，必須求出式(9)～式(11)中的6個速率常數(shù)?；羝娼鸷秃振憷柰ㄟ^電壓膜片鉗實驗，把軸突膜鉗制到不同的電壓得到電導(dǎo)隨時間的變化曲線，通過曲線擬合求得各速率常數(shù)的表達(dá)式[1]，在此不做進一步介紹。

1.2 鈣依賴通道模型

H-H模型中的鈉通道和鉀通道都是電壓門控通道，但某些離子通道同時受跨膜電壓和鈣離子濃度的影響[2-3]，因此將鈣離子濃度整合到 H-H模型中能夠更加完整地描述某些神經(jīng)元的特性。Traub建立的海馬區(qū)神經(jīng)元模型中引入了一個代表胞內(nèi)鈣離子濃度的狀態(tài)變量[4]，這個模型中包含了一種鈣依賴的慢鉀通道Gs，它的速率常數(shù)公式同時依賴于跨膜電壓V和胞內(nèi)的鈣離子濃度χ。慢鉀通道的電導(dǎo)公式為，其中是慢鉀通道的最大電導(dǎo)率，q是一個狀態(tài)變量，它的一階動力學(xué)方程為

式中，αq和βq是速率常數(shù)，其同時依賴于跨膜電壓V和胞內(nèi)鈣離子濃度χ：

為了理解鈣依賴離子通道中鈣離子的作用，一般認(rèn)為在緊靠神經(jīng)元膜內(nèi)側(cè)處有一個扁平空間，其中含有大量的鈣離子，其濃度要明顯高于細(xì)胞內(nèi)其他區(qū)域的濃度。鈣離子主要通過膜上的鈣通道進入這個扁平區(qū)域，然后依靠濃度差擴散到胞內(nèi)其他區(qū)域。可以用一個微分方程來模擬胞內(nèi)鈣離子的動力學(xué)特性[4]：

式中，ICa代表鈣離子通道形成的鈣離子電流，A是一個常數(shù)，代表庫倫轉(zhuǎn)到離子摩爾量的轉(zhuǎn)換系數(shù)，B是速率常數(shù)，代表鈣離子的擴散速率。

在Traub的研究中建立的包含鈣依賴鉀通道的CA3海馬神經(jīng)元模型能重現(xiàn)實驗中在胞體和頂樹突記錄到的動作電位，同時也能重現(xiàn)連續(xù)刺激下的重復(fù)放電現(xiàn)象，從而能夠?qū)A3海馬神經(jīng)元進行更加精確的模擬。

Ekeberg等建立的七鰓鰻運動神經(jīng)元網(wǎng)絡(luò)模型中也考慮了鈣依賴的鉀離子通道[5]，研究發(fā)現(xiàn)，鈣依賴鉀離子電流對動作電位的后超極化時相起主要作用，從而對發(fā)放率起到重要的調(diào)制作用。Yamada等人對鈣依賴通道的機理及建模方法也進行了深入的研究[6]，在此不做詳細(xì)介紹。

1.3 單離子通道馬爾科夫模型

H-H模型對于描述神經(jīng)元的宏觀電流非常成功，然而如果要對單個離子通道電流進行模擬就必須采用其他的方法。在微觀層面上，單個離子通道的開放與關(guān)閉是一個隨機過程，電壓門在自由態(tài)和非自由態(tài)之間的轉(zhuǎn)變依賴于離子通道在不同構(gòu)型間的變化，特定的構(gòu)型能使離子通過而其他的構(gòu)型則不能。通過對單個離子通道進行膜片鉗記錄，就能觀測到離子通道在開放與關(guān)閉狀態(tài)間的隨機改變[7-8]。

Destexhe等提出的馬爾科夫模型為描述單個離子通道的微觀電流提供了理論框架[9]。馬爾科夫模型的基本假設(shè)是離子通道的關(guān)閉與開放，可以被描述為離子通道在一系列不同構(gòu)型狀態(tài)間的改變，不同狀態(tài)間的改變有相應(yīng)的轉(zhuǎn)移概率。假設(shè)離子通道處于狀態(tài)Si時的概率為Pi(t)，pij為處于狀態(tài)i的離子通道轉(zhuǎn)移到狀態(tài)j的條件概率，那么Pi(t)可以表示為：

式中右邊的第1項代表離子通道由其他狀態(tài)轉(zhuǎn)移到Si狀態(tài)時概率的增加量，第2項代表離子通道從狀態(tài)Si轉(zhuǎn)移到其他狀態(tài)時概率的減小量。如果有大量相同類型的離子通道，那么Pi(t)可以被解釋為離子通道處于狀態(tài)Si的比率，而pij就可以被解釋為不同狀態(tài)間轉(zhuǎn)變的速率常數(shù)。因此通過馬爾科夫模型就能夠?qū)⑽⒂^離子通道的特性與H-H模型描述的宏觀電流聯(lián)系起來。

圖2表示了一個包含五個狀態(tài)的馬爾科夫模型，這個模型對應(yīng)了H-H模型中的鉀離子通道。該模型中n0狀態(tài)表示離子通道的四個電壓門都處于非自由態(tài)，n1狀態(tài)表示離子通道有一個電壓門處于自由態(tài)，依次類推，n4狀態(tài)表示離子通道的四個電壓門都處于自由態(tài)，此時離子可以通過離子通道。不同狀態(tài)間的轉(zhuǎn)移概率可以通過H-H模型中的速率常數(shù)計算得到，比如從n0狀態(tài)轉(zhuǎn)移到 n1狀態(tài)共有4種途徑，因此相應(yīng)的轉(zhuǎn)移概率為4αn。最終單個離子通道所處狀態(tài)隨時間變化的規(guī)律可以通過蒙特卡洛模擬得到[9]。

圖2 代表H-H模型中鉀電導(dǎo)的馬爾科夫模型，該模型中的n0到n3為關(guān)閉狀態(tài)，n4為開放狀態(tài)Fig.2 A Markov model of the HH K+conductance.The Markov model has four closed states n0～ n3 and one open state n4

對于鈉離子通道，Destexhe等也通過相同的方法研究[9]。在H-H模型中，假設(shè)鈉通道的電壓門之間都是相互獨立的，在馬爾科夫模型中也基于該假設(shè)進行了研究。但更精細(xì)的研究發(fā)現(xiàn)鈉通道電壓門的激活和失活過程并不是完全獨立的，因此Patlak提出了一個包含電壓門間相互作用的鈉離子馬爾科夫模型[10]，此模型能更好地對通過膜片鉗技術(shù)獲得的實驗數(shù)據(jù)進行描述。

1.4 突觸模型

文中前述內(nèi)容主要對電壓門控離子通道的建模方法進行了介紹，而配體門控的化學(xué)突觸在神經(jīng)信號的傳遞中同樣具有重要的作用。當(dāng)動作電位到達(dá)突觸前末梢時，神經(jīng)遞質(zhì)被釋放到突觸間隙，然后遞質(zhì)與突觸后膜的配體門控受體結(jié)合導(dǎo)致離子電流流入膜內(nèi)。

為了模擬突觸的活動，一般把突觸遞質(zhì)的釋放、擴散、與受體結(jié)合等一系列細(xì)節(jié)抽象為突觸后的一個時間依賴的配體門控通道，Rall描述了配體門控通道的電導(dǎo)方程[11]。假設(shè)在 tspike時刻一個動作電位到達(dá)突觸，則突觸電導(dǎo)方程可表示為：

式中，τsyn為時間常數(shù)，突觸電導(dǎo)在 t=tspike+τsyn時刻達(dá)到最大值gpeak，突觸產(chǎn)生的突觸電流為Isyn(t)=Gsyn(t)(V-Esyn)，其中Esyn為突觸的反轉(zhuǎn)電位。

當(dāng)突觸被一系列的動作電位激活時，將每個動作電位引起的電導(dǎo)改變進行代數(shù)相加就能得到電導(dǎo)的整體變化，然而這種方法會造成計算復(fù)雜性的上升因而很少被用于較復(fù)雜的模擬中。有許多有效的方法可以簡化計算，比如用二階動力學(xué)方程而不是一階動力學(xué)方程來表示電導(dǎo)的變化[12];或者是重組計算順序，Sriinaivasan等基于突觸配體門控通道模型提出了一種算法，對每個突觸只存儲相繼動作電位產(chǎn)生的電導(dǎo)變化量之和，而不是將所有動作電位產(chǎn)生的電導(dǎo)變化都進行存儲再求和[13]，這種算法能提高計算速度，可以用來研究大尺寸神經(jīng)元網(wǎng)絡(luò)中短時程突觸抑制的集合效應(yīng)。

另一種對突觸電導(dǎo)進行建模的方法是Destexhe等人提出的馬爾科夫模型[14]。此模型中的配體門控通道僅僅有開放和關(guān)閉兩種狀態(tài)，可以表示為：

式中，C是關(guān)閉狀態(tài)，O是開放狀態(tài)，T代表神經(jīng)遞質(zhì)，α和β是速率常數(shù)并且獨立于跨膜電位。突觸電流可表示為：

式中，神經(jīng)遞質(zhì)的濃度為[T]。模擬神經(jīng)遞質(zhì)濃度的一種簡化方法是假定動作電位到達(dá)突觸前末梢時，遞質(zhì)釋放量隨時間的變化為幅值不變的脈沖，通過數(shù)值計算就能求得r隨時間變化的方程[15]。

以上討論的模型中遞質(zhì)與受體結(jié)合會直接導(dǎo)致通道的開放形成離子電流，然而有許多神經(jīng)元胞體和樹突含有代謝型受體，這些受體與遞質(zhì)結(jié)合不會導(dǎo)致通道的開放而是激活細(xì)胞內(nèi)的第二信使，經(jīng)過一系列的胞內(nèi)生物化學(xué)過程最終導(dǎo)致離子通道的開放。對這種級聯(lián)化學(xué)反應(yīng)可以用馬爾科夫模型對每一步電化學(xué)反應(yīng)進行模擬，最后將所有的模型進行整合來實現(xiàn)[15]。

包含突觸的神經(jīng)元模型已經(jīng)被廣泛地應(yīng)用于神經(jīng)科學(xué)的研究中，Rattay對中樞系統(tǒng)神經(jīng)元進行了建模，模型中整合了突觸模型，用于研究人工電刺激對中樞神經(jīng)元產(chǎn)生的影響[16]。McIntyre等建立的丘腦皮層中繼神經(jīng)元的模型中包含了興奮性突觸和抑制性突觸，此模型被用于深部腦刺激的機理研究中[17]。

1.5 多室模型

圖1的等效電路可以用來模擬神經(jīng)元細(xì)胞膜的一個局部區(qū)域，然而神經(jīng)元細(xì)胞膜的特性并非各處都是一樣的，不同的膜區(qū)域可能有不同的幾何形狀、不同的離子通道種類及數(shù)量、不同的受體種類及數(shù)量，在神經(jīng)元的軸突上也可能包裹著髓鞘。多室模型是將神經(jīng)元在空間上分割成一系列的小室，每個小室具有相同的生物物理學(xué)特性，并且可以通過如圖1中的等效電路來表示，不同小室等效電路的電學(xué)參數(shù)根據(jù)室與室間生物物理學(xué)參數(shù)的不同而有所不同，室與室之間通過胞內(nèi)的軸向電流聯(lián)系，多室模型的等效電路如圖3所示。室i的跨膜電位Vi與相鄰室的跨膜電位Vi-1和Vi+1的關(guān)系為：

式中，ri±1，i代表相鄰室間的胞內(nèi)軸向電阻，(Vi±1-Vi)/ri±1，i代表胞內(nèi)軸向電流[18-19]。利用多室模型能夠表示任意復(fù)雜形狀的細(xì)胞結(jié)構(gòu)，包括有髓鞘及無髓鞘的神經(jīng)軸突、復(fù)雜的樹突、帶有分叉的軸突及樹突等[20-21]。

圖3 單個神經(jīng)元多室模型的等效電路Fig.3 Electrical equivalent circuit of the compartmental approach for single-neuron modeling

多室模型現(xiàn)在已成為神經(jīng)科學(xué)中電生理模型的重要建模手段，并被廣泛地用于神經(jīng)科學(xué)的研究之中。自從Fitzhugh在1962年首次提出神經(jīng)元有髓鞘軸突的多室模型后[18]，后人對多室模型進行了不斷的擴展研究。Durand等在1992年用多室模型研究得出朗飛氏節(jié)去極化原因一是來自于朗飛氏節(jié)上的激活函數(shù)，二是來自于相鄰朗飛氏節(jié)的電流分布[22]。Foulmeister等在1997年建立了蠑螈視網(wǎng)膜神經(jīng)節(jié)細(xì)胞的多室模型，模型中包含了鈣離子通道及不同種類的鉀離子通道，他們用這個模型研究了神經(jīng)節(jié)細(xì)胞的編碼機理[23]。Parini等在2000年通過對比不同的多室模型，建立了適合于人類視神經(jīng)纖維的多室模型，通過與實驗數(shù)據(jù)對比，證實了模型的可行性[24]。McIntyre等在2002年提出了一種雙層電纜多室模型，并對哺乳動物的運動神經(jīng)纖維進行建模[25-26]，此模型對朗飛氏節(jié)間包裹著髓鞘的部分進行了精細(xì)模擬，模型的特性與實驗結(jié)果有很好的吻合，他們用此模型研究了哺乳動物運動神經(jīng)元興奮時。產(chǎn)生去極化后動作電位和超極化后動作電位的機理，證實了一種慢鉀通道在其中起了重要作用，隨后在2004年，他們又用這種模型模擬了丘腦中繼神經(jīng)元，并對深部腦刺激的作用機理進行了研究[17]。Rattay等在2004年采用多室模型對視網(wǎng)膜中的雙極細(xì)胞和神經(jīng)節(jié)細(xì)胞進行了建模[27]，模擬結(jié)果表明在外部電流下神經(jīng)節(jié)細(xì)胞的軸突是首先興奮的位置。Oozeer等在2005和2006年通過對比前人的多種模型，提出了符合人類視神經(jīng)纖維的多室模型，通過對外部電極刺激產(chǎn)生的外部電場進行三維仿真，研究了在外部電場作用下視神經(jīng)纖維的興奮特性[28-29]。Schiefer等在2006年對單個視網(wǎng)膜神經(jīng)節(jié)細(xì)胞進行了建模[30]，研究表明在外部電刺激下神經(jīng)節(jié)細(xì)胞首先興奮的位置是在神經(jīng)節(jié)軸突的彎曲部位。Sotiropoulous等在2007年對McIntyre提出的雙層電纜多室模型進行了改進[31]，研究了深部腦刺激中電場對神經(jīng)元的直接影響。Bellinger等在2008年建立了包含軸突外鉀離子濃度變化的有髓鞘神經(jīng)纖維計算模型[32]，用來研究深部腦刺激的直接作用，仿真結(jié)果表明軸突外鉀離子的累積會對軸突的興奮產(chǎn)生功能性的阻斷。Jiang等在2009年建立了大鼠丘腦錐體神經(jīng)元的完整模型[33]，用來研究興奮性突觸和抑制性突觸對神經(jīng)元興奮共同作用的內(nèi)在機理，通過模型和實驗的研究他們得出了谷氨酸能興奮性輸入和GABA能抑制性輸入空間加和作用的數(shù)學(xué)規(guī)律。Morse等在2010年用多室模型建立了丘腦CA1神經(jīng)元的完整模型[34]，用來研究遠(yuǎn)端樹突引起的超興奮特性，結(jié)果表明動作電位在樹突的反向傳播導(dǎo)致了超興奮的產(chǎn)生。

1.6 神經(jīng)元網(wǎng)絡(luò)模型

前面幾部分介紹了離子通道、突觸、單個神經(jīng)元的建模方法，利用這些方法可以直接建立神經(jīng)元網(wǎng)絡(luò)的模型。在神經(jīng)元網(wǎng)絡(luò)中神經(jīng)沖動從一個神經(jīng)元通過突觸連接傳遞到下一個神經(jīng)元，通過多室模型對神經(jīng)元軸突、樹突、胞體進行建模，然后通過突觸模型模擬突觸連接就可以建立整個神經(jīng)元網(wǎng)絡(luò)的模型。在建立網(wǎng)絡(luò)級的模型時，主要考慮兩方面的問題，一是要有一種合適的數(shù)學(xué)方法來描述動作電位在神經(jīng)元之間的傳遞，其次是如何確定網(wǎng)絡(luò)中的突觸連接。

在描述動作電位在神經(jīng)元之間的傳遞時有兩種數(shù)學(xué)算法，一種稱為同步算法(時間驅(qū)動算法)[35]，通過這種算法，所有的神經(jīng)元在每一個時鐘時刻都進行一次模擬計算，求得各處的膜電位，這種算法非常適合如基于向量的數(shù)值計算工具如Matlab、Scilab等。另一種算法稱為異步算法(事件驅(qū)動算法)，在這種算法中所有的神經(jīng)元只在接收到一個動作電位或發(fā)放一個動作電位時進行一次模擬計算。異步算法跟同步算法相比更加復(fù)雜，它的優(yōu)勢在于當(dāng)神經(jīng)沖動沒有到達(dá)時不需要對每個神經(jīng)元都進行模擬計算，因此計算速度更快[36-38]。目前，異步算法已經(jīng)被用于很多復(fù)雜的模型中[39-44]。

在建立神經(jīng)元網(wǎng)絡(luò)模型時需要考慮的第二個問題是明確神經(jīng)元之間的連接，小的網(wǎng)絡(luò)模型可以精確地確定每一個突觸連接的細(xì)節(jié)，然而對大的網(wǎng)絡(luò)模型就必須采用一定的策略來簡化模型。比如無脊椎動物海兔的小的神經(jīng)網(wǎng)絡(luò)中僅包含10個神經(jīng)元，每個神經(jīng)元有5個突觸，總共約有50個突觸，因此可以精確地確定每個突觸的受體類型、反轉(zhuǎn)電位、最大電導(dǎo)和傳導(dǎo)延遲等電生理學(xué)特性[45];然而，一個哺乳動物大腦皮層中的局部神經(jīng)元網(wǎng)絡(luò)就可能包含一萬級別的神經(jīng)元數(shù)量，如果每個神經(jīng)元有100個突觸，總共就會有一百萬個突觸連接，要精確地確定每個突觸連接的細(xì)節(jié)是非常困難的，必須采用一定的方法來簡化模型[46]。一種方法是對于特定種類神經(jīng)元(比如中間神經(jīng)元)和一定半徑內(nèi)的另一種神經(jīng)元(比如皮層錐體神經(jīng)元)之間的突觸連接都采用同一種突觸(比如GABA能受體突觸)，這樣就能極大地簡化網(wǎng)絡(luò)模型，并且能較好地符合生理學(xué)結(jié)構(gòu)。

以上介紹的神經(jīng)元網(wǎng)絡(luò)的建模方法在需要精確描述動作電位在神經(jīng)元間的傳遞時比較有效[47]，然而如果用這種建模方法模擬一個較大的神經(jīng)元網(wǎng)絡(luò)時會有巨大的計算量。如果不需要知道神經(jīng)元的內(nèi)部機制，而只關(guān)心神經(jīng)元的輸入輸出關(guān)系，則可采用神經(jīng)元的簡化模型構(gòu)建，比如人工神經(jīng)網(wǎng)絡(luò)。這種模型的目的不在于闡明神經(jīng)回路的機制，而在于構(gòu)造出具有特定功能的神經(jīng)網(wǎng)絡(luò)模型，以解決工程上的實際問題。

2 建模工具

目前已經(jīng)有復(fù)雜的仿真軟件工具包來幫助我們實現(xiàn)依賴于生物物理學(xué)細(xì)節(jié)的神經(jīng)元及網(wǎng)絡(luò)建模。兩種最廣泛應(yīng)用的神經(jīng)元建模工具為GENESIS(GEneral NEural SImulation System)[48-49]和NEURON[50]。這兩種建模環(huán)境可以用來建立神經(jīng)元多室模型，都提供了從分子級別到網(wǎng)絡(luò)級別的建模工具，并且有預(yù)定義的各種功能模塊及可視化圖形用戶界面。功能模塊是指依賴于神經(jīng)元結(jié)構(gòu)特性的各個神經(jīng)片段、電壓門控通道、配體門控通道、化學(xué)突觸等，通過編程語言將相應(yīng)的功能模塊進行連接就能建立神經(jīng)元的模型，同時建模工具還提供了模擬膜片鉗實驗的方法，可以進行各種的模擬實驗。目前很多的神經(jīng)建模工作都是通過GENESIS或NEURON來進行的[51-54]，應(yīng)用范圍包括胞內(nèi)記錄、樹突信號傳遞、神經(jīng)沖動傳導(dǎo)、信息編碼以及學(xué)習(xí)記憶等。

其他較常用的仿真工具還有NEST、NCS、CSIM、XPPAUT、SPLIT和 Mvaspike。NEST主要用于建立大的神經(jīng)元網(wǎng)絡(luò)，這個軟件可以非常容易地建立105數(shù)量級的網(wǎng)絡(luò)，并且可以在每個神經(jīng)元加上符合生理學(xué)數(shù)量的突觸[55];NCS適合于對哺乳動物新皮層的神經(jīng)元結(jié)構(gòu)進行建模;CSIM可以用來模擬多達(dá)數(shù)千個到十萬數(shù)量級突觸的神經(jīng)元網(wǎng)絡(luò);XPPAUT是一種模擬和分析動態(tài)系統(tǒng)的數(shù)值計算工具[56];SPLIT是一種專門對基于H-H方程的大規(guī)模多室模型進行建模的工具[42];Mvaspike是一種利用事件驅(qū)動算法來進行神經(jīng)元網(wǎng)絡(luò)建模的工具。概括來說，可以將各種仿真工具依據(jù)它們最適宜的應(yīng)用范圍分為四類：(1)單室模型：CSIM、NEST和 NCS;(2)多室模型：NEURON、GENESIS、SPLIT;(3)事件驅(qū)動模擬：Mvaspike;(4)動態(tài)系統(tǒng)分析：XPPAUT。

3 總結(jié)

本文對基于生物物理學(xué)參數(shù)的神經(jīng)元電生理模型的建模方法進行了綜述，并介紹了采用相關(guān)建模方法的研究進展以及常用的建模工具。本文可以為從事神經(jīng)建模研究的科研人員提供參考。

神經(jīng)建模作為一種技術(shù)手段已經(jīng)成為神經(jīng)科學(xué)中非常重要的研究方向，已經(jīng)并將繼續(xù)為神經(jīng)科學(xué)的研究做出重要貢獻。然而，神經(jīng)建模也有它的局限性，即使建模過程依據(jù)了神經(jīng)元膜的生物物理學(xué)特性，這仍是對復(fù)雜神經(jīng)元結(jié)構(gòu)的高度抽象。模型建立的準(zhǔn)確性依賴于對重要細(xì)節(jié)的描述，然而在一種神經(jīng)元中重要的細(xì)節(jié)可能在另一種神經(jīng)元中并不那么重要，因此不能認(rèn)為神經(jīng)元模型就是對神經(jīng)元的完整描述。神經(jīng)元的建模過程就是一個從實驗觀察到提出假設(shè)(建立模型)，到模型預(yù)測(模型仿真)，再到模型檢測(同實驗數(shù)據(jù)比較)的不斷迭代的過程。

[1]Hodgkin AL， Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve[J].The Journal of Physiology，1952，117(4)：500-544.

[2]Brown AM，Westenbroek RE，Catterall WA，et al.Axonal L-type Ca2+channels and anoxic injury in rat CNS white matter[J].Journal of Neurophysiology，2001，85(2)：900-911.

[3]Stys P，Lopachin R.Mechanisms of calcium and sodium fluxes in anoxic myelinated central nervous system axons [J].Neuroscience，1997，82(1)：21-32.

[4]Traub R.Simulation of intrinsic bursting in CA3 hippocampal neurons[J].Neuroscience，1982，7(5)：1233-1242.

[5]Ekeberg ?，Wallen P，Lansner A，et al.A computer based model for realistic simulations of neural networks[J].Biological Cybernetics，1991，65(2)：81-90.

[6]Yamada WM，Koch C，Adams PR.Multiple channels and calcium dynamics[M]//Koch C，Segev I.Methods in Neuronal Modeling：From Synapses to Networks.(2nd Edition).Cambridge：MIT Press，1998：137-170.

[7]Neher E，Sakmann B.Single-channel currents recorded from membrane of denervated frog muscle fibres[J].Nature，1976，260(5554)：779-802.

[8]Colquhoun D，Hawkes A.On the stochastic properties of single ion channels[J].Proceedings of the Royal Society of London.Series B，Biological Sciences，1981，211(1183)：205-235.

[9]Destexhe A，Huguenard J.Which formalism to use for modeling voltage-dependent conductances? [M]// Schutter ED.Computational Neuroscience： Realistic Modeling for Experimentalists.Boca Raton：CRC Press，2001：129-157.

[10]Patlak J.Molecular kinetics of voltage-dependent Na+channels[J].Physiological Reviews，1991，71(4)：1047-1080.

[11]Rall W.Distinguishing theoretical synaptic potentials computed for different soma-dendritic distributions of synaptic input[J].Journal of Neurophysiology，1967，30(5)：1138-1168.

[12]Wilson MA，Bower JM.The simulation of large-scale neural networks[M]//Koch C，Segev I.Methods in Neuronal Modeling：From Synapses to Networks.Cambridge：MIT Press，1989：291-334.

[13]Srinivasan R， Chiel HJ. Fast calculation of synaptic conductances[J].Neural Computation，1993，5(2)：200-204.

[14]Destexhe A，Mainen ZF，Sejnowski TJ.Kinetic models of synaptic transmission[M]//Koch C，Segev I.Methods in Neuronal Modeling：From Synapses to Networks(2nd edition).Cambridge：MIT Press，1998：1-26.

[15]Destexhe A，Mainen ZF，Sejnowski TJ.Synthesis of models for excitable membranes，synaptic transmission and neuromodulation using a common kinetic formalism[J].Journal of Computational Neuroscience，1994，1(3)：195-230.

[16]Rattay F.Analysis of the electrical excitation of CNS neurons[J].IEEE Transactions on Biomedical Engineering，1998，45(6)：766-772.

[17]McIntyre C，Grill W，Sherman D，et al.Cellular effects of deep brain stimulation：model-based analysis of activation and inhibition[J].Journal of Neurophysiology，2004，91(4)：1457-1469.

[18]Fitzhugh R.Computation of impulse initiation and saltatory conduction in a myelinated nerve fiber[J]. Biophysical Journal，1962，2(1)：11-21.

[19]McNeal D.Analysis of a model for excitation of myelinated nerve[J].IEEE Transactions on Biomedical Engineering，1976，23(4)：329-337.

[20]Segev I.Cable and compartmental models of dendritic trees[M]//Bower JM，Beeman D.The Book of GENESIS：Exploring Realistic Neural Models with the General Neural Simulation System.New York：Springer-Verlag，1998：51-77.

[21]Zhou XL， Chiu XS. Computer model for action potential propagation through branch point in myelinated nerves[J].Journal of Neurophysiology，2001，85(1)：197-210.

[22]Durand D.Modeling the effects of electric fields on nerve fibers：determination of excitation thresholds [J].IEEE Transactions on Biomedical Engineering，1992，39(12)：1244-1254.

[23]Fohlmeister J， Miller R. Impulse encoding mechanisms of ganglion cells in the tiger salamander retina[J].Journal of Neurophysiology，1997，78(4)：1935-1947.

[24]Parrini S，Delbeke J，Legat V，et al.Modelling analysis of human optic nerve fibre excitation based on experimental data[J].Medical and Biological Engineering and Computing，2000，38(4)：454-464.

[25]McIntyre C，Richardson A，Grill W.Modeling the excitability of mammalian nerve fibers：influence of afterpotentials on the recovery cycle[J].Journal of Neurophysiology，2002，87(2)：995-1006.

[26]McIntyre CC，Grill WM.Extracellular stimulation of central neurons：influence of stimulus waveform and frequency on neuronal output[J].Journal of Neurophysiology，2002，88(4)：1592-1604.

[27]Rattay F， Resatz S. Effective electrode configuration for selective stimulation with inner eye prostheses[J]. IEEE Transactions on Biomedical Engineering，2004，51(9)：1659-1664.

[28]Oozeer M，Veraart C，Legat V，et al.Simulation of intra-orbital optic nerve electrical stimulation[J].Medical and Biological Engineering and Computing，2005，43(5)：608-617.

[29]Oozeer M， Veraart C， Legat V， et al. A model of the mammalian optic nerve fibre based on experimental data[J].Vision Research，2006，46(16)：2513-2524.

[30]Schiefer M，Grill W.Sites of neuronal excitation by epiretinal electrical stimulation [J].IEEE Transactions on Neural Systems and Rehabilitation Engineering，2006，14(1)：5-13.

[31]Sotiropoulos S，Steinmetz P.Assessing the direct effects of deep brain stimulation using embedded axon models[J].Journal of Neural Engineering，2007，4(2)：107-119.

[32]Bellinger SC，Miyazawa G，Steinmetz PN.Submyelin potassium accumulation may functionally block subsets of local axons during deep brain stimulation：a modeling study[J].Journal of Neural Engineering，2008，5(3)：263-274.

[33]Hao J，Wang XD，Dan Y，et al.An arithmetic rule for spatial summation of excitatory and inhibitory inputs in pyramidal neurons[J].Proc Natl Acad Sci USA，2009，106(51)：21906-21911.

[34]Morse TM， Carnevale NT， Mutalik PG， et al.Abnormal excitability of oblique dendrites implicated in early alzheimer's：a computational study [J].Front Neural Circuits，2010，4(16)：1-11

[35]Rotter S，Diesmann M.Exact digital simulation of time-invariant linear systems with applications to neuronal modeling[J].Biological Cybernetics，1999，81(5)：381-402.

[36]Fujimoto RM.Parallel Simulation： Parallel and Distributed Simulation Systems[M].New York：Wiley，2001：147-157.

[37]Zeigler BP，Praehofer H，Kim TG.Theory of Modeling and Simulation[M].(2 nd Edition).New York：Academic Press，2000.

[38]Sloot P， Kaandorpa JA， Hoekstra G， et al.. Distributed simulation with cellular automata：Architecture and applications[C]//Pavelka J，Tel G，Bartosek M，eds.SOFSEM’99：Theory and Practice of Informatics.Berlin：Springer-Verlag，1999：203-248.

[39]Makino T.A discrete-event neural network simulator for general neuron models[J].Neural Computing ＆ Applications，2003，11(3)：210-223.

[40]Gerstner W，Kistler WM.Mathematical formulations of Hebbian learning[J].Biological Cybernetics，2002，87(5)：404-415.

[41]Djurfeldt M，Johansson C，?rjan Ekeberg R.Massively parallel simulation of brain-scale neuronal network models[R].TRITA-NA-P0513，2005.

[42]Brette R.Exact simulation of integrate-and-fire models with exponential currents[J].Neural Computation，2007，19(10)：2604-2609.

[43]Izhikevich EM.Simple model of spiking neurons [J].IEEE Trans Neural Networks，2003，14(6)：1569-1572.

[44]Brette R，Gerstner W.Adaptive exponential integrate-and-fire model as an effective description of neuronal activity[J].Journal of Neurophysiology，2005，94(5)：3637-3642.

[45]White JA，Ziv I，Cleary LJ，et al.The role of interneurons in controlling the tail-withdrawal reflex in Aplysia：a network model[J].Journal of Neurophysiology，1993，70(5)：1777-1786.

[46]Wang Xiaojing， Buzsaki G.Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model[J].J Neurosci，1996，16(20)：6402-6413.

[47]Manor Y， Gonczarowski J， Segev I.Propagation of action potentials along complex axonal trees.Model and implementation[J].Biophysical Journal，1991，60(6)：1411-1423.

[48]Bower JM，Beeman D，Wylde AM.The book of GENESIS：exploring realistic neural models with the GEneral NEural SImulation System[M].New York：Springer-Verlag，1998.

[49]Bower JM，Beeman D，Hucka M.GENESIS simulation system[M]//The Handbook of Brain Theory and Neural Networks.(2 nd Edition).Cambridge：MIT Press，2002：475-478.

[50]Hines ML， Carnevale NT. The NEURON simulation environment[J].Neural Computation，1997，9(6)：1179-1209.

[51]Botta P，de Souza FMS，Sangrey T，et al.Alcohol excites cerebellar golgi cells by inhibiting the Na+/K+ATPase [J].Neuropsychopharmacology，2010，35(9)：1984-1996.

[52]Hendrickson E，Edgerton J，Jaeger D.The capabilities and limitations of conductance-based compartmental neuron models with reduced branched or unbranched morphologies and active dendrites[J].Journal of Computational Neuroscience，2011，30(2)：301-321.

[53]Edgerton JR，Hanson JE，Gunay C，et al.Dendritic sodium channels regulate network integration in globus pallidus neurons：a modeling study[J].J Neurosci，2010，30(45)：15146-15159.

[54]Uebachs M，Opitz T，Royeck M，et al.Efficacy loss of the anticonvulsant carbamazepine in mice lacking sodium channel beta subunits via paradoxical effects on persistent sodium currents[J].J Neurosci，2010，30(25)：8489-8501.

[55]Morrison A， Mehring C， Geisel T， et al. Advancing the boundaries of high-connectivity network simulation with distributed computing [J].Neural Comput，2005，17(8)：1776-801.

[56]Ermentrout B. Simulating， Analyzing， and Animating Dynamical Systems：A Guide to XPPAUT for Researchers and Students[M].Philadelphia：SIAM，2004.