神經元是構成生物體神經系統的基本單位。在神經元的電生理特性描述上, Hodgkin-Huxley(HH)模型最接近于生物學實際。硬件實現神經元可為臨床脊髓損傷治療、仿生學、人工智能等方面提供新的研究思路。本文以HH模型神經元為研究對象, 采用DSP Builder技術, 在現場可編程門陣列(FPGA)上完成了單個HH模型神經元的硬件實現。對FPGA實現的神經元施加了不同類型的電流刺激, 分析研究其動作電位響應特性, 并與相同參數刺激下的數值模擬結果進行了相關系數的計算分析。結果表明, FPGA實現的神經元動作電位響應與數值模擬結果高度一致, 為后續通過硬件構建神經元網絡奠定了基礎。
引用本文： 王金龍, 逯邁, 胡延文, 陳小強, 潘強強. 基于現場可編程門陣列的Hodgkin-Huxley模型神經元動作電位的數值模擬功能實現. 生物醫學工程學雜志, 2015, 32(6): 1302-1309. doi: 10.7507/1001-5515.20150231 復制
0 引言
19世紀50年代，Hodgkin與Huxley發表了一系列文章，首次提出了描述神經元動作電位發放與傳導的數學模型，即Hodgkin-Huxley(HH)模型[1-4]。由于這一杰出的貢獻，他們獲得了1963年諾貝爾生理學與醫學獎。此后，HH模型便成了神經元模型的典型代表。神經系統是一個復雜的非線性系統，現階段對神經元本身特性及神經系統的研究是神經信息學研究領域的熱點[5-8]，其中神經元是構成生物體神經系統的基本單位，故對神經元及神經元網絡的研究是神經系統研究的基礎。近期，通過硬件去實現神經元功能是一個研究創新點[9-10]，在硬件實現人工神經網絡方面，學者們已有很多的研究[11-14]，但是在以生物系統神經元為研究對象的計算神經科學領域，此方面的研究才剛剛開始。
在硬件實現神經元方面，早在1989年，Mead用模擬電路設計并制造了神經元網絡模型，也首次提出了神經形態芯片的創意[15]；1991年，Mahowald等[16]設計出一個基于HH模型的硅神經元，該神經元是最早的硅神經元之一；1998年，Toumazou等[17]用模擬電路實現了HH模型神經元；2011年，Folowosele等[18]用CMOS(complementary metal oxide semiconductor)工藝實現了Mihalas-Niebur模型神經元，再現了該神經元具有的豐富的非線性動力學特性。
如前所述，硬件實現神經元的方式均為搭建模擬電路的方式，相對于數字電路實現方式，模擬電路實現神經元的方式耗費大、耗時長，所以現階段用數字電路去實現神經元是主流。而具有并行處理方式的現場可編程門陣列(field programmable gate array, FPGA)的優勢在此時就顯現出來，但FPGA在硬件實現神經元方面的研究還較少。2004年，Graas等[19]在研究HH模型時使用了Xilinx生產的FPGA芯片，該方法可以極大地提高模擬復雜單神經元及其網絡的速度；2007年，Weinstein等[20-21]對FitzHugh-Nagumo神經元進行了硬件實現，并分析了其特性；2011年，張榮華等[22-23]對Morris-Lecar神經元模型進行了FPGA仿真，再現了Morris-Lecar神經元網絡的非線性動力學特性。
本文通過FPGA輔助開發工具箱DSP Builder與FPGA開發平臺QuartusⅡ聯合設計的方法，對HH模型神經元進行了硬件實現，在FPGA上再現了神經元功能，并對FPGA工作頻率進行了約束，對輸出進行了相應的設置。通過對硬件實現的神經元施加不同類型的電流刺激，將其動作電位響應與數值模擬結果進行分析對比，驗證了硬件實現的神經元具有和數值模擬結果一樣幅值、時間尺度的動作電位響應。
1 HH模型
在計算神經科學領域，HH模型是描述眾多神經元模型中最接近真實生物系統神經元的模型。單個神經元的HH模型由下面的“四個一階”微分方程組描述：
| $\left\{ \begin{gathered} C\frac{{{\text{d}}V}}{{{\text{d}}t}}={G_{{\text{Na}}}}{m^3}h\left({{E_{{\text{Na}}}}-V} \right) + {G_{\text{K}}}{n^4}\left({{E_{\text{K}}}-V} \right) + \hfill \\ {G_1} \pm \left({{E_1}-V} \right) + I \hfill \\ \frac{{{\text{d}}n}}{{{\text{d}}t}}={\alpha _n}\left({1-n} \right)-{\beta _n}n \hfill \\ \frac{{{\text{d}}h}}{{{\text{d}}t}}={\alpha _h}\left({1-h} \right)-{\beta _h}h \hfill \\ \frac{{{\text{d}}m}}{{{\text{d}}t}}={\alpha _m}\left({1-m} \right)-{\beta _m}m \hfill \\ \end{gathered} \right.$ |
式中，V是神經元膜電位，C是膜電容，I是外部刺激電流，GNa、GK分別為鈉離子通道、鉀離子通道最大電導，Gl為漏電導，m為鈉通道活化過程參數，h為鈉通道失活過程參數，n為鉀通道活化過程參數。α函數和β函數是與膜電位有關而與時間無關的速率常數，αm、βm、αh、βh、αn、βn的計算表達式如下[5]：
| $\begin{gathered} {\alpha _m}=\frac{{0.1\left({V + 40} \right)}}{{1-\exp \left({-\left({V + 40} \right)/10} \right)}}, \hfill \\ {\beta _m}=4\exp \left({-\left({V + 65} \right)/18} \right), \hfill \\ {\alpha _h}=0.07\exp \left({-\left({V + 65} \right)/20} \right), \hfill \\ {\beta _h}=\frac{1}{{1-\exp \left({-\left({V + 35} \right)/10} \right) + 1}}, \hfill \\ {\alpha _n}=\frac{{0.01\left({V + 55} \right)}}{{1-\exp \left({-\left({V + 55} \right)/10} \right)}}, \hfill \\ {\beta _n}=0.125\exp \left({-\left({V + 65} \right)/80} \right) \hfill \\ \end{gathered} $ |
HH模型各參數取值如表 1所示。
表1
HH模型的各參數取值
Table1.
Parameter values in HH model
V、m、n、h四個變量的初值選擇，V的初值為神經元細胞膜的靜息電位，本文選V的初值為-65 mV，m、n、h的初始值由下式算出：
| $\begin{gathered} {m_o}\left(V \right)=\frac{{{\alpha _m}}}{{{\alpha _m} + {\beta _m}}}, \hfill \\ {n_o}\left(V \right)=\frac{{{\alpha _n}}}{{{\alpha _n} + {\beta _n}}}, \hfill \\ {h_o}\left(V \right)=\frac{{{\alpha _h}}}{{{\alpha _h} + {\beta _h}}} \hfill \\ \end{gathered} $ |
2 FPGA
本文所用FPGA開發板為睿智FPGA開發板，FPGA芯片型號為Altera公司CycloneⅡ系列EP2C8Q208C8。FPGA一般使用硬件描述語言VHDL(Very-High-Speed Integrated Circuit Hardware Description Language)或Verilog HDL進行開發，開發平臺為QuartusⅡ。Altera公司還提供了FPGA輔助設計開發工具箱DSP Builder，本文選擇DSP Builder進行主系統的設計，對一些外部設計用硬件描述語言進行開發。所使用各軟件版本為QuartusⅡ11.0、Altera DSP Builder 11.0以及Matlab2009b。
3 DSP Builder系統建模
3.1 基本原理
HH模型是一個多變量耦合的復雜非線性系統，如式(1)所示，它由一個“四個一階”微分方程組描述。FPGA實現的數字神經元即是對HH方程進行硬件方式解析，以一個典型的微分方程來說明：
| $\left\{ \begin{gathered} y'\left(t \right)=f\left({t, y\left(t \right)} \right) \hfill \\ y\left({{t_0}} \right)={y_0} \hfill \\ \end{gathered} \right.$ |
式中，f函數為微分方程，y0是微分方程中變量y的初值。運用歐拉法對其進行離散求解，則有：y(n+1)=yn+ΔT*f(tn, yn)(n=0, 1, 2…)
在Simulink環境下，可根據上式的原理用DSP Builder庫文件對上式進行模型搭建，原理框圖如圖 1所示，主要是實現數字積分功能。
圖1
數字積分實現框圖
Figure1.
Block diagram of digital integral
3.2 HH模型離散化
首先選擇合適的方法對神經元模型進行離散化，本文選擇歐拉算法對HH模型進行離散。沒有選擇精度相對較高的龍格庫塔法去離散模型，是因為后期在硬件實現時要考慮FPGA的硬件資源消耗問題，如果使用龍格庫塔法離散的模型在FPGA上實現，勢必會造成大量的硬件資源消耗，甚至導致一些硬件資源較少的FPGA芯片無法完成運算。而相對來說歐拉法離散模型則需要較少的硬件消耗，且在步長選擇合適時，一樣能保持計算精度。用歐拉法對式(1)進行離散化處理，可以得到以下的差分方程。
| $\left\{ \begin{gathered} V\left({k + 1} \right)=V\left(k \right) + \Delta T * \left(\begin{gathered} \left({{G_{{\text{Na}}}}{m^3}h} \right)\left({{E_{{\text{Na}}}}-V\left(k \right)} \right) + \hfill \\ {G_{\text{K}}}{n^4}\left({{E_{\text{K}}}-V\left(k \right) + {G_1}\left({{E_1}-V\left(k \right)} \right) + I} \right) + D/C \hfill \\ \end{gathered} \right) \hfill \\ n\left({k + 1} \right)=n\left(k \right) + \Delta T * {\alpha _n}\left({1-n\left(k \right)} \right)-{\beta _n}n\left(k \right) \hfill \\ h\left({k + 1} \right)=h\left(k \right) + \Delta T * {\alpha _h}\left({1-h\left(k \right)} \right)-{\beta _h}h\left(k \right) \hfill \\ m\left({k + 1} \right)=m\left(k \right) + \Delta T * {\alpha _m}\left({1-m\left(k \right)} \right)-{\beta _m}m\left(k \right) \hfill \\ \end{gathered} \right.$ |
式中，V(k)、n(k)、m(k)、h(k)是4個變量的迭代序列，ΔT是離散化采樣時間。DSP Buildere根據以上離散化后的離散化模型完成具體的模型搭建。
3.3 DSP Builder建模
對HH模型神經元的硬件實現，本質上就是將離散化后的HH模型差分方程用數字信號處理的方式在FPGA上予以計算，結合3.1節所描述的基本原理以及圖 1的數字積分框圖，理解把握4個耦合變量的關系，下面以m變量的計算來說明模塊搭建過程，m變量的計算主要涉及三個方程：
| $\begin{gathered} \frac{{{\text{d}}m}}{{{\text{d}}t}}={\alpha _m}\left({1-m} \right)-{\beta _m}m, \hfill \\ {\alpha _m}=\frac{{0.1\left({V + 40} \right)}}{{1-\exp \left({-\left({V + 40} \right)/10} \right)}}, \hfill \\ {\beta _m}=4\exp \left({-\left({V + 65} \right)/18} \right) \hfill \\ \end{gathered} $ |
由離散化后的差分方程[式(5)]可知，在狀態k時刻，m變量的差分方程如下式：
| $m\left({k + 1} \right)=m\left(k \right) + \Delta T * {\alpha _m}\left({1-m\left(k \right)} \right)-{\beta _m}m\left(k \right)$ |
αm、βm的計算式變為：
| $\begin{gathered} {\alpha _m}=\frac{{0.1\left({V\left(k \right) + 40} \right)}}{{1-\exp \left({-\left({V\left(k \right) + 40} \right)/10} \right)}}, \hfill \\ {\beta _m}=4\exp \left({-\left({V\left(k \right) + 65} \right)/18} \right) \hfill \\ \end{gathered} $ |
通過式(7)可以由前一個狀態的值(m(k))算出下一個狀態的值(m(k+1))，逐次循環，可以算出所有離散點的值，根據式(7)，在DSP Builder中搭建計算模型，如圖 2所示。
圖2
m變量離散積分模塊
Figure2.
m variables discrete integral module
結合實現數字積分功能原理框圖(如圖 1所示)，圖 2為一個離散積分模塊。它可以完成輸入In1(dm(k)/dt)到輸出Out1(m(k+1))的計算。中間的虛線框內是一個初值選擇模塊，初值選擇模塊使用了一個比較輸入器(Multiplexer)、一個信號脈沖(Single Pulse)和一個固定值模塊，原理是在計算開始時，信號脈沖輸入值為0，則比較輸入器的0端導通，輸入固定值，此值即為變量的初值，在經過幾個延時之后，信號脈沖輸入值變為1，此時1端導通，完成m當前值的輸入。圖中Memory Delay起寄存器的作用。
m變量計算子系統如圖 3所示。圖中Subsystem的輸入In1為αm(1-m(k))-βmm(k)，此式結果由V(k)算出，圖中的Subsystem為前面搭建的離散積分模塊(圖 2)，前面的模塊就是對In1端αm(1-m(k))-βmm(k)的計算，查找表(Look up tables)LUT和LUT1分別完成膜電位V到αm、βm的映射。
圖3
m變量計算子系統圖
Figure3.
m variable calculation subsystem
以LUT為例(LUT負責計算αm的值)，說明查找表的參數設置，αm的計算式為：
| ${\alpha _m}=\frac{{0.1\left({V\left(k \right) + 40} \right)}}{{1-\exp \left({1-\left({V\left(k \right) + 40} \right)/10} \right)}}$ |
由公式可以看出α函數和β函數只與膜電位有關，α函數和β函數的計算式中含有指數函數以及除法，在DSP Builder中沒有合適的模塊來計算，所以在搭建模型時α函數和β函數的計算使用了查找表，作用就是用輸入V來計算α函數和β函數的值，先將輸入V轉換為地址，再由地址查找出在查找表中相應的值，查找表地址位寬為8位，共256個地址，輸入V值限定-80 mV(地址0)到50 mV(地址255)之間，地址(add)和膜電位V的關系式為：
| ${\text{add=}}\left({V-\left({-80} \right)} \right) * \frac{{{2^8}-1}}{{50-\left({-80} \right)}}=1.965 * \left({V + 80} \right)$ |
如圖 3所示，輸入add即是按式(10)將V處理后的地址值。根據式(10)，查找表LUT的MATLAB array項應該設置為：0.1*((-40+80):-130/2^8:(-40-50))./(exp(((-40+80):-130/2^8:(-40-50))/10)-1)。通過LUT可完成輸入V算出αm，βm的計算可由LUT1同理算出。由m變量子系統的搭建，可同理搭出n、h變量計算的子系統，有了m、n、h三個變量計算的子系統圖，可由HH模型離散化后的差分方程(式(4))綜合搭建出HH模型的總系統圖。
3.4 時序控制以及硬件資源消耗
在使用DSP Builder技術進行相關設計建模時，很重要的一點就是要把握好各變量的時序問題，否則很容易造成時序混亂，導致計算錯誤，為此文獻[20]提出了一種使用圖形化的設計方案，將各變量時序關系用流水線結構圖的形式表示出來，這種方法使各變量之間時序關系清楚明了，很有利于建模的成功。本文按照類似的方法，設計了HH模型建模時的流水線結構圖(圖 4)，并以此作為DSP Builder建模時的時序參考，實踐證明，這種方法是很有效的。
圖4
流水線結構圖
(a)主流水線結構圖；(b) (a)中虛線框內流水線結構圖
Figure4. Pipeline construction(a) the main pipeline structure; (b) pipeline architecture of dashed frame in (a)
芯片資源消耗量是FPGA相關設計工程中必須加以考慮的問題，因為設計方案所消耗的芯片資源消耗量直接關系到實驗可行度、實驗成本各方面的問題，而設計的計算精度和芯片消耗量又是互相制約的，這就需要找出一個平衡點，盡量減少資源消耗的同時保持計算的精度。如表 2所示，是本文設計中的FPGA資源消耗情況，在表中可以看出本次設計芯片內部的乘法器使用率達到了100%，基本邏輯單元與存儲器還有剩余。由于HH模型比較復雜，乘法運算比較多，在乘法器不足而又需要進行大量乘法運算時，可以用減少乘法器數據位寬的方式來實現。復雜的乘法運算可以通過查找表來完成，查找表在綜合時消耗的資源為Memory Bits與M4ks。要減少查找表綜合消耗的資源，可以通過減少查找表地址位寬以及減少查找表數據位寬的方式來實現。
表2
FPGA資源消耗量
Table2.
FPGA resource usage
4 結果與分析
首先對HH模型進行數值模擬仿真，對其施加不同的方波和正弦波兩種電流刺激，得到其數值模擬結果的動作電位波形。然后對相同參數設置下的硬件實現的神經元施加與數值模擬一致的電流刺激，在示波器上觀察其動作電位響應。將硬件實現的神經元動作電位響應與數值模擬結果進行對比，并取樣計算結果之間的相關系數。實驗結果如圖 5、圖 6所示。
圖5
方波刺激波形圖
刺激電流：(a)幅值10μA/cm2，寬度30 ms；(b)幅值40μA/cm2，寬度30 ms；(c)幅值10μA/cm2，寬度60 ms；(a1)、(b1)、(c1)分別為與之對應的硬件實現結果
Figure5. Waveform under square wave stimulationthe stimulus current: (a) amplitude 10μA/cm2, width 30 ms; (b) amplitude 40μA/cm2, width 30 ms; (c) amplitude 10μA/cm2, width 60 ms; (a1), (b1), (c1) are the corresponding hardware implementation waveform
圖6
正弦波刺激數值模擬波形圖
刺激電流：(a)峰值10μA/cm2，頻率30 Hz；(b)幅值10μA/cm2，頻率50 Hz；(c)幅值40μA/cm2，頻率30 Hz；(a1)、(b1)、(c1)分別為與之對應的硬件實現結果
Figure6. Numerical simulation waveform under sine wave stimulationThe stimulus current: (a) peak 10μA/cm2, frequency 30 Hz; (b) peak 10μA/cm2, frequency 50 Hz; (c) peak 40μA/cm2, frequency 30 Hz; (a1), (b1), (c1) are the corresponding hardware implementation waveform
圖 5是方波電流刺激下數值模擬波形圖及FPGA硬件實現波形圖。在方波刺激下，當刺激幅值不變而改變寬度時，動作電位的發放數會隨刺激寬度的加寬而變多；而當刺激寬度不變，改變刺激幅值時，動作電位發放數會隨幅值的增大而變多。由于篇幅所限，還有一些未列出的實驗結果，如當刺激幅值較小時，即當刺激電流低于閾值刺激時，或在神經元不應期進行刺激時，不產生動作電位，以及當刺激幅值較大時，會使動作電位產生振蕩，最后趨于停止。這些結果都與數值模擬結果一致。
圖 6是正弦波電流刺激數值模擬波形圖及硬件實現波形圖。可以看出，在正弦波刺激下，固定刺激正弦波的幅值而改變其頻率時，神經元會出現不同的放電模式。例如，當刺激正弦波峰值固定為10μA/cm2，頻率從30 Hz增大到50 Hz時，神經元在每個正弦波周期內均產生一個動作電位；而當頻率固定為30 Hz，峰值從10μA/cm2增大到40μA/cm2，神經元放電從一個周期一個動作電位變成了一個周期兩個動作電位，所以說，正弦刺激下神經元動作電位的發放與刺激的峰值和頻率均有關。
為了量化比較硬件實現結果與數值模擬結果的一致性，引入變量間的相關系數[24]，相關系數公式如式(11)，
| $\begin{gathered} {r_{xy}}=\frac{{\sum\limits_{k=1}^n {\left({\left({{x^k}-\left\langle {{x^k}} \right\rangle } \right)\left({{y^k}-\left\langle {{y^k}} \right\rangle } \right)} \right)} }}{{\sqrt {\sum\limits_{k=1}^n {{{\left({{x^k}-\left\langle {{x^k}} \right\rangle } \right)}^2} * } \sum\limits_{k=1}^n {{{\left({{y^k}-\left\langle {{y^k}} \right\rangle } \right)}^2}} } }} \hfill \\ \left({k=1, 2 \cdots n} \right) \hfill \\ \end{gathered} $ |
其中，xk、yk代表x、y兩變量的離散樣本，〈xk〉、〈yk〉是對x、y兩變量所有離散樣本數的平均。相關系數r介于-1到1之間，相關系數r的絕對值越大表示兩變量間的相關程度越高，即變量的一致性越高。當r的絕對值等于1時，表示兩變量完全相關。根據式(11)，對圖 5與圖 6的結果進行相關系數的計算。方波刺激只選取有刺激部分的結果樣本進行計算，正弦波刺激選取0~100 ms之間的樣本值進行計算。計算結果如表 3所示，可以發現無論是方波刺激還是正弦波刺激，FPGA實現的神經元動作電位響應結果與數值模擬結果的相關系數值基本都大于0.7，屬于高度線性相關，驗證了此設計方法的正確性。
表3
實驗結果的相關系數
Table3.
Correlation coefficient of the results
5 結束語
本文主要通過DSP Builder技術，對HH模型神經元進行了硬件實現，在FPGA上再現了神經元功能，將硬件實現的結果與數值模擬結果進行了對比并計算了相關系數。較高的相關系數值表明，本文的硬件實現方法是正確有效的。由于FPGA芯片資源的限制，本文只完成了單個神經元的FPGA實現，后續的研究工作就是在這個基礎上完成大規模神經元網絡的硬件實現，完成其特性分析研究，并為相關的工程應用提供一定的借鑒。
0 引言
19世紀50年代，Hodgkin與Huxley發表了一系列文章，首次提出了描述神經元動作電位發放與傳導的數學模型，即Hodgkin-Huxley(HH)模型[1-4]。由于這一杰出的貢獻，他們獲得了1963年諾貝爾生理學與醫學獎。此后，HH模型便成了神經元模型的典型代表。神經系統是一個復雜的非線性系統，現階段對神經元本身特性及神經系統的研究是神經信息學研究領域的熱點[5-8]，其中神經元是構成生物體神經系統的基本單位，故對神經元及神經元網絡的研究是神經系統研究的基礎。近期，通過硬件去實現神經元功能是一個研究創新點[9-10]，在硬件實現人工神經網絡方面，學者們已有很多的研究[11-14]，但是在以生物系統神經元為研究對象的計算神經科學領域，此方面的研究才剛剛開始。
在硬件實現神經元方面，早在1989年，Mead用模擬電路設計并制造了神經元網絡模型，也首次提出了神經形態芯片的創意[15]；1991年，Mahowald等[16]設計出一個基于HH模型的硅神經元，該神經元是最早的硅神經元之一；1998年，Toumazou等[17]用模擬電路實現了HH模型神經元；2011年，Folowosele等[18]用CMOS(complementary metal oxide semiconductor)工藝實現了Mihalas-Niebur模型神經元，再現了該神經元具有的豐富的非線性動力學特性。
如前所述，硬件實現神經元的方式均為搭建模擬電路的方式，相對于數字電路實現方式，模擬電路實現神經元的方式耗費大、耗時長，所以現階段用數字電路去實現神經元是主流。而具有并行處理方式的現場可編程門陣列(field programmable gate array, FPGA)的優勢在此時就顯現出來，但FPGA在硬件實現神經元方面的研究還較少。2004年，Graas等[19]在研究HH模型時使用了Xilinx生產的FPGA芯片，該方法可以極大地提高模擬復雜單神經元及其網絡的速度；2007年，Weinstein等[20-21]對FitzHugh-Nagumo神經元進行了硬件實現，并分析了其特性；2011年，張榮華等[22-23]對Morris-Lecar神經元模型進行了FPGA仿真，再現了Morris-Lecar神經元網絡的非線性動力學特性。
本文通過FPGA輔助開發工具箱DSP Builder與FPGA開發平臺QuartusⅡ聯合設計的方法，對HH模型神經元進行了硬件實現，在FPGA上再現了神經元功能，并對FPGA工作頻率進行了約束，對輸出進行了相應的設置。通過對硬件實現的神經元施加不同類型的電流刺激，將其動作電位響應與數值模擬結果進行分析對比，驗證了硬件實現的神經元具有和數值模擬結果一樣幅值、時間尺度的動作電位響應。
1 HH模型
在計算神經科學領域，HH模型是描述眾多神經元模型中最接近真實生物系統神經元的模型。單個神經元的HH模型由下面的“四個一階”微分方程組描述：
| $\left\{ \begin{gathered} C\frac{{{\text{d}}V}}{{{\text{d}}t}}={G_{{\text{Na}}}}{m^3}h\left({{E_{{\text{Na}}}}-V} \right) + {G_{\text{K}}}{n^4}\left({{E_{\text{K}}}-V} \right) + \hfill \\ {G_1} \pm \left({{E_1}-V} \right) + I \hfill \\ \frac{{{\text{d}}n}}{{{\text{d}}t}}={\alpha _n}\left({1-n} \right)-{\beta _n}n \hfill \\ \frac{{{\text{d}}h}}{{{\text{d}}t}}={\alpha _h}\left({1-h} \right)-{\beta _h}h \hfill \\ \frac{{{\text{d}}m}}{{{\text{d}}t}}={\alpha _m}\left({1-m} \right)-{\beta _m}m \hfill \\ \end{gathered} \right.$ |
式中，V是神經元膜電位，C是膜電容，I是外部刺激電流，GNa、GK分別為鈉離子通道、鉀離子通道最大電導，Gl為漏電導，m為鈉通道活化過程參數，h為鈉通道失活過程參數，n為鉀通道活化過程參數。α函數和β函數是與膜電位有關而與時間無關的速率常數，αm、βm、αh、βh、αn、βn的計算表達式如下[5]：
| $\begin{gathered} {\alpha _m}=\frac{{0.1\left({V + 40} \right)}}{{1-\exp \left({-\left({V + 40} \right)/10} \right)}}, \hfill \\ {\beta _m}=4\exp \left({-\left({V + 65} \right)/18} \right), \hfill \\ {\alpha _h}=0.07\exp \left({-\left({V + 65} \right)/20} \right), \hfill \\ {\beta _h}=\frac{1}{{1-\exp \left({-\left({V + 35} \right)/10} \right) + 1}}, \hfill \\ {\alpha _n}=\frac{{0.01\left({V + 55} \right)}}{{1-\exp \left({-\left({V + 55} \right)/10} \right)}}, \hfill \\ {\beta _n}=0.125\exp \left({-\left({V + 65} \right)/80} \right) \hfill \\ \end{gathered} $ |
HH模型各參數取值如表 1所示。
表1
HH模型的各參數取值
Table1.
Parameter values in HH model
V、m、n、h四個變量的初值選擇，V的初值為神經元細胞膜的靜息電位，本文選V的初值為-65 mV，m、n、h的初始值由下式算出：
| $\begin{gathered} {m_o}\left(V \right)=\frac{{{\alpha _m}}}{{{\alpha _m} + {\beta _m}}}, \hfill \\ {n_o}\left(V \right)=\frac{{{\alpha _n}}}{{{\alpha _n} + {\beta _n}}}, \hfill \\ {h_o}\left(V \right)=\frac{{{\alpha _h}}}{{{\alpha _h} + {\beta _h}}} \hfill \\ \end{gathered} $ |
2 FPGA
本文所用FPGA開發板為睿智FPGA開發板，FPGA芯片型號為Altera公司CycloneⅡ系列EP2C8Q208C8。FPGA一般使用硬件描述語言VHDL(Very-High-Speed Integrated Circuit Hardware Description Language)或Verilog HDL進行開發，開發平臺為QuartusⅡ。Altera公司還提供了FPGA輔助設計開發工具箱DSP Builder，本文選擇DSP Builder進行主系統的設計，對一些外部設計用硬件描述語言進行開發。所使用各軟件版本為QuartusⅡ11.0、Altera DSP Builder 11.0以及Matlab2009b。
3 DSP Builder系統建模
3.1 基本原理
HH模型是一個多變量耦合的復雜非線性系統，如式(1)所示，它由一個“四個一階”微分方程組描述。FPGA實現的數字神經元即是對HH方程進行硬件方式解析，以一個典型的微分方程來說明：
| $\left\{ \begin{gathered} y'\left(t \right)=f\left({t, y\left(t \right)} \right) \hfill \\ y\left({{t_0}} \right)={y_0} \hfill \\ \end{gathered} \right.$ |
式中，f函數為微分方程，y0是微分方程中變量y的初值。運用歐拉法對其進行離散求解，則有：y(n+1)=yn+ΔT*f(tn, yn)(n=0, 1, 2…)
在Simulink環境下，可根據上式的原理用DSP Builder庫文件對上式進行模型搭建，原理框圖如圖 1所示，主要是實現數字積分功能。
圖1
數字積分實現框圖
Figure1.
Block diagram of digital integral
3.2 HH模型離散化
首先選擇合適的方法對神經元模型進行離散化，本文選擇歐拉算法對HH模型進行離散。沒有選擇精度相對較高的龍格庫塔法去離散模型，是因為后期在硬件實現時要考慮FPGA的硬件資源消耗問題，如果使用龍格庫塔法離散的模型在FPGA上實現，勢必會造成大量的硬件資源消耗，甚至導致一些硬件資源較少的FPGA芯片無法完成運算。而相對來說歐拉法離散模型則需要較少的硬件消耗，且在步長選擇合適時，一樣能保持計算精度。用歐拉法對式(1)進行離散化處理，可以得到以下的差分方程。
| $\left\{ \begin{gathered} V\left({k + 1} \right)=V\left(k \right) + \Delta T * \left(\begin{gathered} \left({{G_{{\text{Na}}}}{m^3}h} \right)\left({{E_{{\text{Na}}}}-V\left(k \right)} \right) + \hfill \\ {G_{\text{K}}}{n^4}\left({{E_{\text{K}}}-V\left(k \right) + {G_1}\left({{E_1}-V\left(k \right)} \right) + I} \right) + D/C \hfill \\ \end{gathered} \right) \hfill \\ n\left({k + 1} \right)=n\left(k \right) + \Delta T * {\alpha _n}\left({1-n\left(k \right)} \right)-{\beta _n}n\left(k \right) \hfill \\ h\left({k + 1} \right)=h\left(k \right) + \Delta T * {\alpha _h}\left({1-h\left(k \right)} \right)-{\beta _h}h\left(k \right) \hfill \\ m\left({k + 1} \right)=m\left(k \right) + \Delta T * {\alpha _m}\left({1-m\left(k \right)} \right)-{\beta _m}m\left(k \right) \hfill \\ \end{gathered} \right.$ |
式中，V(k)、n(k)、m(k)、h(k)是4個變量的迭代序列，ΔT是離散化采樣時間。DSP Buildere根據以上離散化后的離散化模型完成具體的模型搭建。
3.3 DSP Builder建模
對HH模型神經元的硬件實現，本質上就是將離散化后的HH模型差分方程用數字信號處理的方式在FPGA上予以計算，結合3.1節所描述的基本原理以及圖 1的數字積分框圖，理解把握4個耦合變量的關系，下面以m變量的計算來說明模塊搭建過程，m變量的計算主要涉及三個方程：
| $\begin{gathered} \frac{{{\text{d}}m}}{{{\text{d}}t}}={\alpha _m}\left({1-m} \right)-{\beta _m}m, \hfill \\ {\alpha _m}=\frac{{0.1\left({V + 40} \right)}}{{1-\exp \left({-\left({V + 40} \right)/10} \right)}}, \hfill \\ {\beta _m}=4\exp \left({-\left({V + 65} \right)/18} \right) \hfill \\ \end{gathered} $ |
由離散化后的差分方程[式(5)]可知，在狀態k時刻，m變量的差分方程如下式：
| $m\left({k + 1} \right)=m\left(k \right) + \Delta T * {\alpha _m}\left({1-m\left(k \right)} \right)-{\beta _m}m\left(k \right)$ |
αm、βm的計算式變為：
| $\begin{gathered} {\alpha _m}=\frac{{0.1\left({V\left(k \right) + 40} \right)}}{{1-\exp \left({-\left({V\left(k \right) + 40} \right)/10} \right)}}, \hfill \\ {\beta _m}=4\exp \left({-\left({V\left(k \right) + 65} \right)/18} \right) \hfill \\ \end{gathered} $ |
通過式(7)可以由前一個狀態的值(m(k))算出下一個狀態的值(m(k+1))，逐次循環，可以算出所有離散點的值，根據式(7)，在DSP Builder中搭建計算模型，如圖 2所示。
圖2
m變量離散積分模塊
Figure2.
m variables discrete integral module
結合實現數字積分功能原理框圖(如圖 1所示)，圖 2為一個離散積分模塊。它可以完成輸入In1(dm(k)/dt)到輸出Out1(m(k+1))的計算。中間的虛線框內是一個初值選擇模塊，初值選擇模塊使用了一個比較輸入器(Multiplexer)、一個信號脈沖(Single Pulse)和一個固定值模塊，原理是在計算開始時，信號脈沖輸入值為0，則比較輸入器的0端導通，輸入固定值，此值即為變量的初值，在經過幾個延時之后，信號脈沖輸入值變為1，此時1端導通，完成m當前值的輸入。圖中Memory Delay起寄存器的作用。
m變量計算子系統如圖 3所示。圖中Subsystem的輸入In1為αm(1-m(k))-βmm(k)，此式結果由V(k)算出，圖中的Subsystem為前面搭建的離散積分模塊(圖 2)，前面的模塊就是對In1端αm(1-m(k))-βmm(k)的計算，查找表(Look up tables)LUT和LUT1分別完成膜電位V到αm、βm的映射。
圖3
m變量計算子系統圖
Figure3.
m variable calculation subsystem
以LUT為例(LUT負責計算αm的值)，說明查找表的參數設置，αm的計算式為：
| ${\alpha _m}=\frac{{0.1\left({V\left(k \right) + 40} \right)}}{{1-\exp \left({1-\left({V\left(k \right) + 40} \right)/10} \right)}}$ |
由公式可以看出α函數和β函數只與膜電位有關，α函數和β函數的計算式中含有指數函數以及除法，在DSP Builder中沒有合適的模塊來計算，所以在搭建模型時α函數和β函數的計算使用了查找表，作用就是用輸入V來計算α函數和β函數的值，先將輸入V轉換為地址，再由地址查找出在查找表中相應的值，查找表地址位寬為8位，共256個地址，輸入V值限定-80 mV(地址0)到50 mV(地址255)之間，地址(add)和膜電位V的關系式為：
| ${\text{add=}}\left({V-\left({-80} \right)} \right) * \frac{{{2^8}-1}}{{50-\left({-80} \right)}}=1.965 * \left({V + 80} \right)$ |
如圖 3所示，輸入add即是按式(10)將V處理后的地址值。根據式(10)，查找表LUT的MATLAB array項應該設置為：0.1*((-40+80):-130/2^8:(-40-50))./(exp(((-40+80):-130/2^8:(-40-50))/10)-1)。通過LUT可完成輸入V算出αm，βm的計算可由LUT1同理算出。由m變量子系統的搭建，可同理搭出n、h變量計算的子系統，有了m、n、h三個變量計算的子系統圖，可由HH模型離散化后的差分方程(式(4))綜合搭建出HH模型的總系統圖。
3.4 時序控制以及硬件資源消耗
在使用DSP Builder技術進行相關設計建模時，很重要的一點就是要把握好各變量的時序問題，否則很容易造成時序混亂，導致計算錯誤，為此文獻[20]提出了一種使用圖形化的設計方案，將各變量時序關系用流水線結構圖的形式表示出來，這種方法使各變量之間時序關系清楚明了，很有利于建模的成功。本文按照類似的方法，設計了HH模型建模時的流水線結構圖(圖 4)，并以此作為DSP Builder建模時的時序參考，實踐證明，這種方法是很有效的。
圖4
流水線結構圖
(a)主流水線結構圖；(b) (a)中虛線框內流水線結構圖
Figure4. Pipeline construction(a) the main pipeline structure; (b) pipeline architecture of dashed frame in (a)
芯片資源消耗量是FPGA相關設計工程中必須加以考慮的問題，因為設計方案所消耗的芯片資源消耗量直接關系到實驗可行度、實驗成本各方面的問題，而設計的計算精度和芯片消耗量又是互相制約的，這就需要找出一個平衡點，盡量減少資源消耗的同時保持計算的精度。如表 2所示，是本文設計中的FPGA資源消耗情況，在表中可以看出本次設計芯片內部的乘法器使用率達到了100%，基本邏輯單元與存儲器還有剩余。由于HH模型比較復雜，乘法運算比較多，在乘法器不足而又需要進行大量乘法運算時，可以用減少乘法器數據位寬的方式來實現。復雜的乘法運算可以通過查找表來完成，查找表在綜合時消耗的資源為Memory Bits與M4ks。要減少查找表綜合消耗的資源，可以通過減少查找表地址位寬以及減少查找表數據位寬的方式來實現。
表2
FPGA資源消耗量
Table2.
FPGA resource usage
4 結果與分析
首先對HH模型進行數值模擬仿真，對其施加不同的方波和正弦波兩種電流刺激，得到其數值模擬結果的動作電位波形。然后對相同參數設置下的硬件實現的神經元施加與數值模擬一致的電流刺激，在示波器上觀察其動作電位響應。將硬件實現的神經元動作電位響應與數值模擬結果進行對比，并取樣計算結果之間的相關系數。實驗結果如圖 5、圖 6所示。
圖5
方波刺激波形圖
刺激電流：(a)幅值10μA/cm2，寬度30 ms；(b)幅值40μA/cm2，寬度30 ms；(c)幅值10μA/cm2，寬度60 ms；(a1)、(b1)、(c1)分別為與之對應的硬件實現結果
Figure5. Waveform under square wave stimulationthe stimulus current: (a) amplitude 10μA/cm2, width 30 ms; (b) amplitude 40μA/cm2, width 30 ms; (c) amplitude 10μA/cm2, width 60 ms; (a1), (b1), (c1) are the corresponding hardware implementation waveform
圖6
正弦波刺激數值模擬波形圖
刺激電流：(a)峰值10μA/cm2，頻率30 Hz；(b)幅值10μA/cm2，頻率50 Hz；(c)幅值40μA/cm2，頻率30 Hz；(a1)、(b1)、(c1)分別為與之對應的硬件實現結果
Figure6. Numerical simulation waveform under sine wave stimulationThe stimulus current: (a) peak 10μA/cm2, frequency 30 Hz; (b) peak 10μA/cm2, frequency 50 Hz; (c) peak 40μA/cm2, frequency 30 Hz; (a1), (b1), (c1) are the corresponding hardware implementation waveform
圖 5是方波電流刺激下數值模擬波形圖及FPGA硬件實現波形圖。在方波刺激下，當刺激幅值不變而改變寬度時，動作電位的發放數會隨刺激寬度的加寬而變多；而當刺激寬度不變，改變刺激幅值時，動作電位發放數會隨幅值的增大而變多。由于篇幅所限，還有一些未列出的實驗結果，如當刺激幅值較小時，即當刺激電流低于閾值刺激時，或在神經元不應期進行刺激時，不產生動作電位，以及當刺激幅值較大時，會使動作電位產生振蕩，最后趨于停止。這些結果都與數值模擬結果一致。
圖 6是正弦波電流刺激數值模擬波形圖及硬件實現波形圖。可以看出，在正弦波刺激下，固定刺激正弦波的幅值而改變其頻率時，神經元會出現不同的放電模式。例如，當刺激正弦波峰值固定為10μA/cm2，頻率從30 Hz增大到50 Hz時，神經元在每個正弦波周期內均產生一個動作電位；而當頻率固定為30 Hz，峰值從10μA/cm2增大到40μA/cm2，神經元放電從一個周期一個動作電位變成了一個周期兩個動作電位，所以說，正弦刺激下神經元動作電位的發放與刺激的峰值和頻率均有關。
為了量化比較硬件實現結果與數值模擬結果的一致性，引入變量間的相關系數[24]，相關系數公式如式(11)，
| $\begin{gathered} {r_{xy}}=\frac{{\sum\limits_{k=1}^n {\left({\left({{x^k}-\left\langle {{x^k}} \right\rangle } \right)\left({{y^k}-\left\langle {{y^k}} \right\rangle } \right)} \right)} }}{{\sqrt {\sum\limits_{k=1}^n {{{\left({{x^k}-\left\langle {{x^k}} \right\rangle } \right)}^2} * } \sum\limits_{k=1}^n {{{\left({{y^k}-\left\langle {{y^k}} \right\rangle } \right)}^2}} } }} \hfill \\ \left({k=1, 2 \cdots n} \right) \hfill \\ \end{gathered} $ |
其中，xk、yk代表x、y兩變量的離散樣本，〈xk〉、〈yk〉是對x、y兩變量所有離散樣本數的平均。相關系數r介于-1到1之間，相關系數r的絕對值越大表示兩變量間的相關程度越高，即變量的一致性越高。當r的絕對值等于1時，表示兩變量完全相關。根據式(11)，對圖 5與圖 6的結果進行相關系數的計算。方波刺激只選取有刺激部分的結果樣本進行計算，正弦波刺激選取0~100 ms之間的樣本值進行計算。計算結果如表 3所示，可以發現無論是方波刺激還是正弦波刺激，FPGA實現的神經元動作電位響應結果與數值模擬結果的相關系數值基本都大于0.7，屬于高度線性相關，驗證了此設計方法的正確性。
表3
實驗結果的相關系數
Table3.
Correlation coefficient of the results
5 結束語
本文主要通過DSP Builder技術，對HH模型神經元進行了硬件實現，在FPGA上再現了神經元功能，將硬件實現的結果與數值模擬結果進行了對比并計算了相關系數。較高的相關系數值表明，本文的硬件實現方法是正確有效的。由于FPGA芯片資源的限制，本文只完成了單個神經元的FPGA實現，后續的研究工作就是在這個基礎上完成大規模神經元網絡的硬件實現，完成其特性分析研究，并為相關的工程應用提供一定的借鑒。
