基于突觸離子通道動力學神經(jīng)元網(wǎng)絡的高效并行仿真算法*

彭 霞，王直杰，韓 芳，顧曉春

(東華大學信息科學與技術學院，上海 201620)

1 引言

生物神經(jīng)網(wǎng)絡是大量的各類神經(jīng)元通過各種突觸相互連接，并相互作用，實現(xiàn)各類神經(jīng)信號的復雜傳遞及處理過程的復雜動力學系統(tǒng)。近年來，為了研究各類神經(jīng)信號復雜傳遞過程，探索和掌握生物神經(jīng)系統(tǒng)的各種功能機理，全球計算神經(jīng)科學專家建立了大量的生物神經(jīng)網(wǎng)絡數(shù)學模型，并通過計算機軟、硬件技術，研發(fā)了大量系統(tǒng)仿真軟件平臺，如GENESIS(the GEneral NEural SImulation System)[1]、NEURON[2,3]、NEST[4]、藍色大腦項目[5]等對神經(jīng)網(wǎng)絡系統(tǒng)進行模擬仿真[6,7]。

生物神經(jīng)元和突觸的動力學特性非常復雜，生物神經(jīng)網(wǎng)絡系統(tǒng)仿真模擬的計算量非常大，譬如，對一個神經(jīng)元個數(shù)N=104的全連接采用遞質-受體(Neurotransmitter-receptor)離子通道動力學[8,9]的復雜突觸模型的生物神經(jīng)網(wǎng)絡，每一突觸后末梢有103個受體，若采用傳統(tǒng)時鐘同步算法[6,10,11]進行該生物系網(wǎng)絡系統(tǒng)的模擬仿真，在每一時間步內，包含N(104)次神經(jīng)元膜電位計算和約N2(108)次突觸電流計算，并且每一突觸電流涉及的受體濃度計算包括受體個數(shù)(103)次綁定和解綁轉化過程的計算。這意味著，每一時間步的每一條突觸的輸入電流計算涉及到大量的復雜離子通道動力學模型的計算，整個系統(tǒng)模擬的計算量尤其是突觸計算量是相當驚人的，單處理器架構下大規(guī)模生物神經(jīng)網(wǎng)絡模擬仿真將變得非常困難。因此，當神經(jīng)網(wǎng)絡的神經(jīng)元及突觸規(guī)模達到一定數(shù)量時，生物神經(jīng)網(wǎng)絡的仿真需要大規(guī)模并行計算技術。因圖形處理器GPU(Graphics Processing Units)[12,13]是強大的計算平臺，且相對于采用大規(guī)模并行計算機而言，它具有價格低、使用方便等特點，本文采用GPU硬件平臺和英偉達NVIDIA公司的統(tǒng)一計算設備架構CUDA(Compute Unified Device Architecture)開發(fā)工具進行生物神經(jīng)網(wǎng)絡并行算法的設計。

為了提高大規(guī)模神經(jīng)元網(wǎng)絡系統(tǒng)在GPU上的并行仿真效率，通過分析整個網(wǎng)絡的計算過程發(fā)現(xiàn)，仿真計算量主要包括每個神經(jīng)元膜電位的計算和更新以及每條突觸電流的計算和更新兩部分；某時刻，某條突觸電流的計算不但與該突觸后神經(jīng)元當前時刻的狀態(tài)相關，與突觸的神經(jīng)遞質傳遞信號時延以及具體突觸特性相關，還跟突觸前神經(jīng)元的所有歷史放電狀態(tài)有關?？梢钥闯觯挥|電流計算中，突觸前神經(jīng)元和突觸后神經(jīng)元的狀態(tài)共同決定突觸的狀態(tài)，突觸前后神經(jīng)元的耦合度非常高。因此，降低突觸電流計算時突觸前后神經(jīng)元之間的耦合度是設計計算模塊獨立性強、模塊間耦合度低的高效生物神經(jīng)網(wǎng)絡并行算法的關鍵。

本文通過分析化學突觸信息傳遞機理及遞質-受體離子通道動力學特征，提出了突觸前神經(jīng)元釋放的遞質濃度計算與突觸后神經(jīng)元已綁定受體濃度計算分離的思想，即遞質-受體計算分離思想，大大增強了這兩個計算過程的獨立性，有效降低突觸電流計算中突觸前神經(jīng)元和突觸后神經(jīng)元的耦合度，并基于此思想設計一種計算模塊獨立性強、模塊間耦合度低的并行算法，提高了生物神經(jīng)網(wǎng)絡并行仿真效率。

2 具有復雜突觸動力學的生物神經(jīng)網(wǎng)絡模型

本文以全連接網(wǎng)絡結構、IAF (Integrate-And-Fire)神經(jīng)元模型[14,15]及遞質-受體離子通道動力學突觸模型的生物神經(jīng)網(wǎng)絡模型為例設計并行仿真算法(其它類型的網(wǎng)絡模型的仿真算法可類似設計，如神經(jīng)元為H-H(Hodgkin-Huxley)方程[16]描述的模型)。

2.1 神經(jīng)元計算模型

基于IAF神經(jīng)元模型的全連接神經(jīng)網(wǎng)絡可用動態(tài)方程(1)表示。

(1)

通過Euler數(shù)值計算方法對公式(1)進行差分計算，設步長為h=Δt，則神經(jīng)元在時間步t+1的膜電位Vi(t+1)可由公式(2)來計算。

Vi(t+1)=Vi(t)+

(2)

2.2 遞質-受體離子通道動力學突觸計算模型

(3)

其中，th為時間常數(shù)，h(t)在t=th時取峰值h(th)=e-1/th，qmax為遞質分子濃度最大值，t0為突觸前神經(jīng)元的放電時刻，時間常數(shù)tq決定了q(t)的衰減快慢。

根據(jù)化學突觸的動力學機理[8]，突觸傳遞遞質分子T綁定突觸后末梢的各類受體R(Receptor)的動力學過程可用動力學過程(4)表示。

(4)

其中，R表示突觸后末梢上未與遞質分子綁定的受體;TR*表示突觸后末梢已經(jīng)與遞質分子綁定的受體;α是受體從未綁定到與遞質分子綁定的狀態(tài)轉化率常數(shù);β是受體從綁定狀態(tài)轉化到未綁定狀態(tài)的轉化率常數(shù)，轉化率常數(shù)由具體突觸本身特性決定，其值越大，表明在同一時間段內受體狀態(tài)轉換越快。公式(4)表示在每一單位時間Δt，任何一個綁定的受體，有β的可能轉化為未綁定，對于任意未被綁定的受體有αq(t)的可能轉化成已綁定狀態(tài)。

神經(jīng)元j到神經(jīng)元i之間的突觸記為：j→i，t時刻，突觸電導記為Gij(t)，設突觸后神經(jīng)元的末梢的受體總數(shù)為NR，分別分布在突觸后末梢的第k(k=1,…,NR)個位置，受體跟突觸間隙的遞質分子結合并綁定，從而使得突觸后末梢該受體位置的離子通道門打開，gk(t)為k位置的單離子通道門開度，突觸單離子通道門的最大開度即最大電導為gc，已綁定遞質的受體個數(shù)(記為nc(t))直接決定實際整個突觸j→i此時刻的離子通道門的開度。因此，對于每條突觸，在任意t時刻,突觸電導Gij(t)可以用以下公式表示：

(5)

(6)

gc·nc(t)·(Esyn-Vi)

(7)

下面本文將根據(jù)突觸信息傳遞機理及突觸動力學特性來分析突觸電流計算特點，為降低突觸電流計算中突觸前后神經(jīng)元的耦合度，提出遞質-受體計算分離的思想，從而設計一種基于突觸離子通道動力學特征的神經(jīng)元網(wǎng)絡并行算法，實現(xiàn)生物神經(jīng)網(wǎng)絡的高效并行仿真。

3 遞質-受體計算分離思想

本文采用時鐘同步算法進行生物神經(jīng)網(wǎng)絡系統(tǒng)的模擬，傳統(tǒng)意義上的時鐘同步算法框架如圖1所示，在每一時鐘步Δt，針對每一個神經(jīng)元進行神經(jīng)元狀態(tài)的計算，包括兩部分：對神經(jīng)元連接的每條突觸進行輸入突觸電流計算(公式(3)、公式(4)和公式(7))，并累加獲得該神經(jīng)元總的突觸輸入電流；計算每一神經(jīng)元的膜電位(公式(2))。時鐘不斷向前推進，實現(xiàn)整個網(wǎng)絡的仿真。

Figure 1 Flow of neuronal network serial simulation圖1 神經(jīng)元網(wǎng)絡串行仿真流程圖

結合公式(2)～公式(4)和公式(7)，對于由N個神經(jīng)元組成的全連接網(wǎng)絡，在每一時間步，主要包括四個主要計算部分：N(N-1)次突觸前神經(jīng)元產(chǎn)生的遞質分子濃度q(t)計算；對每條突觸要通過NR個隨機過程來統(tǒng)計已綁定的該突觸后末梢受體個數(shù)，并且由文獻[8]可知，NR越大，離子通道動力學突觸模型中突觸電導Gij(t)函數(shù)曲線將越平滑，因此共有NR·N(N-1)次突觸后受體nc(t)計算；N(N-1)次突觸電流計算；N次神經(jīng)元的膜電位Vi(t)計算?？梢钥闯?，當神經(jīng)網(wǎng)絡規(guī)模增大、模型的復雜度增加時，整個網(wǎng)絡的計算量急劇增大，而且絕大多數(shù)的計算都耗費在突觸電流的計算上。因此，大規(guī)模神經(jīng)元網(wǎng)絡的模擬仿真需要并行計算，而且并行算法設計中，突觸電流計算的并行策略設計尤其重要。

然而，從圖1可以看出，從神經(jīng)元j到神經(jīng)元i的突觸電流計算包括三部分：與突觸前神經(jīng)元放電狀態(tài)相關的遞質分子濃度計算；突觸后神經(jīng)元末梢已綁定受體的統(tǒng)計計算；突觸后神經(jīng)元的輸入電流計算。每條突觸電流計算過程中，突觸后神經(jīng)元上已綁定受體個數(shù)的統(tǒng)計計算必須要獲得當前時刻的遞質分子濃度q(t)，而此時刻遞質分子濃度跟突觸前神經(jīng)元歷史放電信息以及突觸信息傳遞時延等有關，這意味著，每條突觸電流計算過程中，突觸前神經(jīng)元和突觸后神經(jīng)元狀態(tài)信息之間耦合度很高。如何減少突觸電流計算中突觸前后神經(jīng)元之間的耦合度，是設計高效神經(jīng)元網(wǎng)絡并行算法的關鍵問題。

Figure 2 Core framework of the algorithm based on neurotransmitter-receptor computing separation圖2 基于遞質-受體計算分離的核心算法框架

從基于遞質-受體計算分離思想的神經(jīng)元網(wǎng)絡核心算法可以看出，通過直接從循環(huán)數(shù)組Q中獲取突觸前神經(jīng)元的信息，突觸電流計算主要由突觸后神經(jīng)元的狀態(tài)決定；而由突觸前神經(jīng)元狀態(tài)決定的遞質分子濃度q(t)，從突觸電流計算中分離出來，通過直接獲取突觸前神經(jīng)元膜電位等信息直接計算并保存到數(shù)組Q中；顯然，每條突觸電流計算時，突觸前神經(jīng)元和突觸后神經(jīng)元狀態(tài)之間的耦合度通過引入循環(huán)數(shù)組Q而得到顯著降低。因此，整個網(wǎng)絡計算量最大、耦合度最高的突觸電流計算模塊的獨立性得到了很大的提高，有利于提高大規(guī)模生物神經(jīng)網(wǎng)絡并行仿真效率。

4 基于遞質-受體計算分離的神經(jīng)元網(wǎng)絡并行算法

Figure 3 Parallel strategy of the algorithm based on neurotransmitter-receptor computing separation圖3 基于遞質-受體計算分離算法的并行策略

根據(jù)圖3所示策略，本文設計了基于遞質-受體計算分離的神經(jīng)元網(wǎng)絡仿真并行算法如下：

t=t0；

whilet

/*GPU上執(zhí)行，每條線程對應一條突觸電流計算，并行實現(xiàn)每個神經(jīng)元的總輸入電流計算*/

(1)計算突觸后末梢的被綁定受體數(shù)nc1：

nc-new=0；//新綁定的受體個數(shù)

nc-left=nc0；//當前時刻已綁定的受體個數(shù)

從j的長度為Dmaxdelay的循環(huán)數(shù)組Q的(j,(t-dij)modDmaxdelay)位置獲取突觸前神經(jīng)元j產(chǎn)生的遞質分子濃度q(t)；

For每個未綁定的受體R(共NR-nc0個)：產(chǎn)生隨機數(shù)，若大于α·q(t)，則該未綁定的受體轉變成已綁定，即：nc-new=nc-new+1；

For每個已綁定的受體TR*(共nc0個)：產(chǎn)生隨機數(shù)，若大于受體的解綁率常數(shù)β，則該綁定受體將解綁成未綁定狀態(tài)，即nc-left=nc0-1；nc1=nc-new+nc-left；

/*GPU上執(zhí)行，每一條線程對應每個神經(jīng)元的膜電位及遞質濃度計算*/

步驟3根據(jù)公式(2)計算神經(jīng)元i在t時刻的膜電位Vi(t)；

步驟4計算突觸前神經(jīng)元狀態(tài)及其產(chǎn)生的遞質濃度：

if(Vi(t)>Vth)，重置神經(jīng)元i的膜電位：Vi(t)=0；根據(jù)公式(3)計算神經(jīng)元i的虛擬突觸在t時刻的遞質濃度q(t-t0)并存儲到循環(huán)數(shù)組Q的(i,tmodDmaxdelay)位置；

else

根據(jù)公式(3)計算神經(jīng)元i的虛擬突觸在t時刻的遞質濃度q(t)并存儲到循環(huán)數(shù)組Q的(i,tmodDmaxdelay)位置；

end if

//CPU上執(zhí)行，時鐘推進

t=t0+Δt；

end while

5 仿真結果及性能分析

本文在單機上分別針對基于遞質-受體計算分離的神經(jīng)元網(wǎng)絡時鐘同步串行算法、傳統(tǒng)時鐘同步算法進行了神經(jīng)元網(wǎng)絡系統(tǒng)模擬仿真實驗，并在單NVIDIA GPU模式下使用CUDA對基于遞質-受體計算分離的神經(jīng)元網(wǎng)絡并行算法進行并行模擬實驗，實驗測試軟件平臺為NVIDIA GeForce GTX 960，4 GB顯存,Windows 8.0,Visual Studio 2015，C++語言，機器主機頻率4.2 GHz，內存8 GB。生物神經(jīng)網(wǎng)絡采用包含興奮性、抑制性兩類IAF神經(jīng)元及遞質-受體離子通道動力學突觸模型(受體個數(shù)NR=1000)的全連接網(wǎng)絡。模擬50 ms(5 000步)生物時間，步長為0.01 ms，網(wǎng)絡中最大延遲步長為2 ms(200步)。

對基于遞質-受體計算分離的時鐘同步串行算法及傳統(tǒng)時鐘同步串行算法分別進行了上述神經(jīng)元網(wǎng)絡系統(tǒng)的模擬仿真實驗，其結果如圖4所示，當網(wǎng)絡的神經(jīng)元數(shù)量大于3 000時，新算法的串行執(zhí)行時間與傳統(tǒng)時鐘同步算法相比，明顯減少，并且隨著網(wǎng)絡規(guī)模的增大，減少的幅度更加明顯。這是因為，對于每一時間步，傳統(tǒng)時鐘同步算法中每條突觸電流的計算包含遞質濃度的計算和受體濃度的計算，遞質濃度q(t)的計算次數(shù)為突觸數(shù)量次，即N·(N-1)，而將遞質濃度計算從突觸電流計算中分離出來，只要在進行每個突觸前神經(jīng)元的膜電位計算時進行該神經(jīng)元引起的遞質濃度計算，因此每一時間步遞質濃度計算次數(shù)減少到神經(jīng)元數(shù)量次，即N。當神經(jīng)元數(shù)量較少時，突觸數(shù)量相對不大，突觸電流計算中的遞質濃度計算量雖然得到減少，但對整個算法的性能改善不明顯；而當網(wǎng)絡規(guī)模擴大時，基于遞質-受體計算分離算法的串行執(zhí)行效率較傳統(tǒng)時鐘同步算法串行執(zhí)行效率優(yōu)勢逐漸明顯。因此可以看出，隨著網(wǎng)絡規(guī)模的擴大，基于遞質-受體計算分離的新時鐘同步算法提高了傳統(tǒng)時鐘同步算法的性能。

Figure 4 Serial simulation results of traditional algorithm and the new algorithm with different neural network scales圖4 不同網(wǎng)絡規(guī)模下的傳統(tǒng)時鐘同步算法 與新算法的串行仿真結果

采用基于遞質-受體計算分離的并行算法在不同規(guī)模下對上述神經(jīng)元網(wǎng)絡模擬實驗的加速比及GPU運行時間如圖5a和圖5b所示(注：分別采用尺寸為640、1 024的一維線程塊)。從圖5a可以看出，本文提出的并行算法對神經(jīng)元網(wǎng)絡并行仿真時，其加速比整體上隨著網(wǎng)絡規(guī)模的增大而快速增加。這是因為隨著整個網(wǎng)絡規(guī)模N的擴大，突觸數(shù)量N·(N-1)迅速增大，每一時間步的突觸電流計算次數(shù)N·(N-1)也將迅速增長，并且每一時間步內突觸電流計算中，突觸受體濃度計算量(每一時間步有NR·N·(N-1)次)急劇增大，使得基于遞質-受體計算分離的串行算法的執(zhí)行時間急劇增大。而基于遞質-受體計算分離的時鐘同步并行算法的映射策略是采用每條突觸電流計算(包括突觸受體濃度計算)對應每條GPU線程，這些線程大量并行執(zhí)行。因此，相對于串行算法執(zhí)行時間增加幅度而言，并行算法在GPU上執(zhí)行時間的增加比較平緩(圖5b)，所以并行算法的加速比增長速度較快(圖5a)。然而當網(wǎng)絡規(guī)模達到6 000或10 000附近(塊尺寸分別為640和1024)的兩條加速比曲線都出現(xiàn)了小幅度的下降(圖5a)。這是GPU硬件結構特點導致的，因為隨著網(wǎng)絡規(guī)模的擴大，神經(jīng)元數(shù)量以及突觸數(shù)量每達到硬件一次上下文切換的線程條數(shù)時，都將要進行上下文切換，這一特點影響了突觸電流并行和神經(jīng)元狀態(tài)計算并行兩模塊的執(zhí)行時間，從而表現(xiàn)出整個算法在某些網(wǎng)絡規(guī)模(6 000、10 000)下出現(xiàn)GPU并行執(zhí)行時間加速上升(圖5b)，導致圖5a中兩條加速比曲線在相應的規(guī)模下出現(xiàn)小幅度的下降，但隨后兩條加速比曲線又隨著網(wǎng)絡規(guī)模的增大而增加。因此，本文提出的并行算法隨著網(wǎng)絡規(guī)模的增大整體上表現(xiàn)出很好的并行性能。

Figure 5 Performance of the parallel algorithm圖5 并行算法的性能

上述仿真實驗結果表明，與傳統(tǒng)時鐘同步算法相比，本文提出的遞質-受體計算分離算法既減少了傳統(tǒng)時鐘同步算法中遞質濃度計算量，提高了傳統(tǒng)時鐘同步算法的性能，還通過將遞質濃度計算從突觸電流計算中分離出來單獨計算，降低了突觸電流計算中突觸前神經(jīng)元和突觸后神經(jīng)元之間狀態(tài)耦合度，更有利于占整個網(wǎng)絡最大計算量的突觸電流計算部分的并行化實現(xiàn)。因此，本文提出的基于遞質-受體計算分離的并行算法對大規(guī)模的含復雜突觸模型的神經(jīng)元網(wǎng)絡的模擬仿真是高效的。

6 結束語

本文通過分析化學突觸遞質分子和受體離子通道動力學特征，提出了遞質-受體計算分離的思想，增強了突觸前神經(jīng)元引起的遞質分子濃度計算與突觸后綁定狀態(tài)的受體計算之間的獨立性，降低突觸電流計算中突觸前神經(jīng)元狀態(tài)和突觸后神經(jīng)元狀態(tài)之間的耦合度?；谶f質-受體計算分離的思想，通過不同的映射策略分別實現(xiàn)了突觸前神經(jīng)元的狀態(tài)信息并行計算和突觸電流的并行計算，從而設計了一種大規(guī)模生物神經(jīng)網(wǎng)絡并行算法。仿真結果表明了該算法的高效性。

參考文獻:

[1] Basu I,Kudela P,Korzeniewska A,et al.A study of the dynamics of seizure propagation across micro domains in the vicinity of the seizure onset zone[J].Journal of Neural Engineering,2015,12(4):046016.

[2] Carnevale N T,Hines M L.The NEURON book[M].Cambridge:Cambridge University Press,2006.

[3] Kim D,Paré D,Nair S S.Mechanisms contributing to the induction and storage of Pavlovian fear memories in the lateral amygdala[J].Learning & Memory,2013,20(8):421-430.

[4] Jan M,Tomohiro S,Kenji D.The mechanism of saccade motor pattern generation investigated by a large-scale spiking neuron model of the superior colliculus[J].Plos One,2013,8(2):e57134.

[5] Markram H.The blue brain project[J].Nature Rev Neurosci,2006,7(2):153-160.

[6] Brette R,Rudolph M,Carnevale T,et al.Simulation of networks of spiking neurons:A review of tools and strategies[J].Journal of Computational Neuroscience,2007,23(3):349-398.

[7] Milkowski M. Explanatory completeness and idealization in large brain simulations:A mechanistic perspective[J].Synthese,2016,193(5):1457-1478.

[8] Chapeau-bloudeau F, Chambet N.Synapse models for neural networks:From ion channel kinetics to multiplicative coefficient wij[J].Nueral Computation,1995,7(4):713-734.

[9] Destexhe A,Mainen Z F,Sejnowski T J.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.

[10] Hanuschkin A,Kunkel S,Helias M,et al.A general and efficient method for incorporating precise spike times in globally time-driven simulations[J].Frontiers in Neuroinformatics,2010,4(8):113.

[11] Lin Xiang-hong.Simulation and evolution of large-scale spiking neural networks[D].Harbin:Harbin Institute of Technology,2009.(in Chinese)

[12] Nageswaran J M,Dutt N,Krichmar J L,et al.2009 special issue:A configurable simulation environment for the efficient simulation of large-scale spiking neural networks on graphics processors[J].Neural Networks the Official Journal of the International Neural Network Society,2009,22(5-6):791-800.

[13] Brette R,Goodman D F.Vectorized algorithms for spiking neural network simulation[J].Neural Computation,2011,23(6):1503-1535.

[14] Dayan P,Abbott L F.Theoretical neuroscience[M]. Massachusetts:MIT Press,2001:154-155.

[15] Brette R.Exact simulation of integrate-and-fire models with synaptic conductances[J].Neural Computation,2006,18(8):2004-2027.

[16] Hodgkin A L,Huxley A F.Currents carried by sodium and potassium ions through the membrane of the giant axon of Loligo[J].Journal of Physiology,1952,116(4):449.

[17] Li Lian-jiang.Research on GPU based spiking neural network learning [D].Wuhan:Huazhong University of Science and Technology,2015.(in Chinese)

附中文參考文獻：

[11] 藺想紅.大規(guī)模脈沖神經(jīng)網(wǎng)絡的模擬與進化研究[D].哈爾濱：哈爾濱工業(yè)大學,2009.

[17] 李連江.基于GPU的脈沖神經(jīng)網(wǎng)絡學習研究[D].武漢：華中科技大學,2015.