利用計算機模擬大腦的神經振蕩對于分析大腦的功能機制具有重要的意義。丘腦皮層神經群模型(TNMM)是一種通過建立丘腦和皮層耦合關系來反映神經活動機制的模型, 有助于我們理解大腦特定認知功能及各個節律腦電信號的神經振蕩活動。隨著神經群模型的復雜性和規模的不斷增加, 常規的計算機系統性能無法實現快速和大規模的模型仿真。本文提出了以FPGA為基礎的硬件計算仿真方法, 利用Altera公司的計算機輔助設計工具DSP Builder, 同時結合MATLAB/Simulink實現復雜神經群算法模型的搭建和硬件平臺移植。此方法充分利用FPGA的并行計算能力來實現大規模復雜神經群模型的快速仿真, 為神經群模型的計算機實現提供了新的方案和思路。
引用本文： 梁振虎, 周敬亮, 李小俚. 丘腦皮層神經群模型仿真及現場可編程門陣列實現研究. 生物醫學工程學雜志, 2016, 33(4): 616-625. doi: 10.7507/1001-5515.20160103 復制
引言
腦電是大腦神經活動的宏觀表現，在不同的腦功能活動和精神疾病狀態下會表現出特定的神經振蕩現象，例如睡眠、癲癇發作[1]、帕金森病[2]等都會表現出特有的神經振蕩活動[3-4]。因此，從神經振蕩角度研究腦電活動對于理解大腦的工作機制具有重要意義。尤其是采用數學和物理模型來模擬大腦的神經振蕩活動對于深入分析大腦的功能機制具有重要的價值[5-6]。
目前，針對大腦的神經活動，從微觀、介觀、宏觀等不同尺度均有建模。其中，微觀模型主要以單個神經元為單位進行建模[7]。從微觀層面對神經元的建模有助于了解大腦最基本的組成單元的工作原理，但是大腦由上千億個不同種類的神經細胞構成，很難確定每個神經元模型的參數，且各種神經元之間的連接和大腦神經系統的結構也異常復雜，大腦的功能也是眾多神經元群協同作用的結果，完成這項工程需要消耗巨大的計算機資源, 僅從微觀層次上建模對整個大腦或者大腦某一區域的研究來說是遠遠不夠的。因此，微觀模型不能真實地反映宏觀的腦電特征，只能用來仿真單個神經元[8-9]。介觀模型是建立微觀和宏觀的橋梁，以集總參數模型為主，對特定種類細胞組成的神經元群整體特性進行建模，例如單一的神經群模型(neural mass model, NMM)和平均場模型(mean-field cortical, MFC)，NMM和MFC起初是用來仿真嗅覺皮層的活動[10]以及自發的alpha節律腦電波[11]，隨后這些模型經過改進和擴展，能夠用較少的狀態變量來描述大量神經元群的平均電活動[12-13]。而宏觀模型則考慮了神經元群之間的耦合關系和腦區間的耦合關系[14-15]，比如多通道神經元群模型[16]。Wendling等[12]利用三個神經群模型的不同耦合方式，產生類癲癇數據，與真實的癲癇數據對比，研究癲癇的發作機制。通過對多通道神經振蕩的建模及其同真實腦電數據之間的對比分析，可以幫助理解大腦區域間的耦合關系以及癲癇、帕金森病、阿爾茨海默癥等的發病機制，驗證各種腦電數據處理算法的有效性。在三種尺度模型的研究方面，以介觀層面的研究為多。尤其是在經典神經群模型的基礎上演化而來的丘腦皮層神經群模型(thalamocortical neural mass model, TNMM)，在睡眠、記憶力鞏固以及節律產生等腦功能認知方面有著重要的作用[17-18]。通過簡單的數學模型可以在介觀層面仿真出大腦皮層與丘腦間的連接機制，對腦功能的研究具有重大價值。
隨著對神經群模型計算仿真研究的不斷深入，神經群規模和復雜度不斷增加，當前通用計算機的計算無法滿足對大規模的復雜神經群模型實現實時仿真的需求。基于并行計算和分布式的計算是解決計算復雜度的重要方法。尤其以現場可編程門陣列(Field Programmable Gate Array，FPGA)為基礎的硬件計算方法在并行計算和求解復雜非線性常微分方程的數值解方面具有一定的優勢，而且易于和MATLAB形成無縫連接。因此，本文利用FPGA模擬復雜的丘腦皮層神經群數學模型，采用Altera公司提供的計算機輔助設計工具DSP Builder聯合Simulink完成丘腦皮層神經群模型的設計，使用FPGA內部定制的NiosⅡ嵌入式處理器完成對硬件系統中外部設備的驅動程序設計，并實現對核心算法模塊參數的控制。
1 FPGA仿真平臺設計
FPGA大都是基于查找表技術，通過燒寫文件改變查找表內容的方法來實現對FPGA的重復配置。FPGA芯片主要由可編程輸入輸出單元、基本可編程邏輯單元、內嵌SRAM、豐富的布線資源、底層嵌入功能單元和內嵌專用單元等組成，可以完成極其復雜的時序與邏輯組合電路功能。基于FPGA的數字建模仿真平臺如圖 1所示。該平臺采用Altera公司的Cyclone Ⅲ系列FPGA芯片作為系統硬件核心，硬件系統中外部設備的驅動程序主要由FPGA內部定制的NiosⅡ嵌入式處理器完成，通過軟件編程實現人機交互以及對模型數據收發的控制，最大程度地利用FPGA內部邏輯單元，使處理器的外部接口更少。這樣可以大大提高系統的性能和穩定性，同時降低研制成本。
圖1
FPGA數字建模仿真平臺結構示意圖
Figure1.
Digital model and simulation platform of FPGA
DSP Builder是一個系統級設計工具，以MATLAB/Simulink的模塊庫形式出現，在Simulink中進行圖形化的設計和仿真，同時把算法仿真建模和硬件實現的設計工具連接起來，使算法的開發到硬件的實現可以無縫地過渡。神經群算法模型主要在DSP Builder中搭建完成，DSP Builder提供了一個從MATLAB/Simulink直接到FPGA硬件實現的設計接口，在Simulink中完成系統集成，然后通過Signal Compiler模塊將電路模型文件即Simulink模塊文件(.mdl)轉換成寄存器傳輸級(register transfer level, RTL)的VHDL代碼表述和工具命令語言(Tool Command Language, TCL)腳本，使其可以在FPGA中完成硬件電路的實現。同時通過QuartusⅡ中的SOPC Builder工具在FPGA內部定制NiosⅡ嵌入式處理器，將在DSP Builder中生成的模型IP核通過Avalon總線與NiosⅡ處理器相連接，從而將軟件的靈活性與硬件的高速性結合起來，完成與其它模塊的數據通信。
2 模型的建立
2.1 經典神經群模型
經典神經群模型[19-20]的基本思想是使興奮性和抑制性細胞群相互作用，產生相應的神經振蕩活動。單個神經群模型主要由錐體細胞(pyramidal cells, PC)、興奮性中間神經元(excitatory interneurons, EI)和抑制性中間神經元(inhibitory interneurons, Ⅱ)三種不同的神經群構成。經典神經群模型框圖如圖 2(a)所示。
圖2
丘腦皮層神經群模型構成原理框圖
(a)經典神經群模型; (b)大腦皮層神經群模型; (c)丘腦神經群模型
Figure2. Schematic diagram of TNMM(a) classical NMM; (b) the NMM of a cortical voxel; (c) the NMM of the thalamus
經典神經群的動力學演化主要包含兩部分。第一部分是將突觸前信息(動作電位的平均脈沖密度)轉變成突觸后信息(突觸后膜興奮或抑制電位)。
在興奮性和抑制性狀況下，線性傳遞函數的脈沖響應分別為：
| $ {h_e}\left( t \right) = Aat{e^{ - at}} $ |
| $ {h_i}\left( t \right) = Bat{e^{ - bt}} $ |
其中, 參數A和B分別為興奮性和抑制性平均突觸增益，a和b分別為興奮性和抑制性細胞平均突觸時間常數，通過調節這四個參數可以改變興奮和抑制突觸的靈敏度，確定模型的興奮和抑制特性。線性函數hi(t)和he(t)可由如下形式的微分方程描述：
| $ \left\{ \begin{gathered} \dot z\left( t \right) = {z_1}\left( t \right) \hfill \\ \dot z\left( t \right) = Ggx\left( t \right) - 2g{z_1}\left( t \right) - {g^2}z\left( t \right) \hfill \\ \end{gathered} \right. $ |
其中，G和g表示興奮性條件下的參數A和a或者抑制條件下的參數B和b，x(t)和z(t)分別為神經群模型的輸入和輸出數據。
第二部分是通過靜態非線性函數將平均膜電壓轉換成動作電位的平均密度，即平均點燃率，并作為線性變換函數的輸入。靜態非線性函數表示為：
| $ S\left( v \right) = 2{e_0}/\left( {1 + {e^{r\left( {{v_0} - v} \right)}}} \right) $ |
其中2e0表示最大點燃率，v0表示相對于點燃率e0的突觸后電位，r表示S型函數的彎曲程度，v表示突觸前平均膜電壓。
經典神經群模型可由以下微分方程表示：
| $ \left\{ \begin{gathered} {{\dot x}_1}\left( t \right) = {x_2}\left( t \right) \hfill \\ {{\dot x}_2}\left( t \right) = AaS\left( {{x_3}\left( t \right) - {x_5}\left( t \right)} \right) - 2a{x_2}\left( t \right) - {a^2}{x_1}\left( t \right) \hfill \\ {{\dot x}_3}\left( t \right) = {x_4}\left( t \right) \hfill \\ {{\dot x}_4}\left( t \right) = Aa\left[ {p + {C_2}S\left( {{C_1}{x_1}\left( t \right)} \right)} \right] - 2a{x_4}\left( t \right) - {a^2}{x_3}\left( t \right) \hfill \\ {{\dot x}_5}\left( t \right) = {x_6}\left( t \right) \hfill \\ {{\dot x}_6}\left( t \right) = Bb{C_4}S\left( {{C_3}{x_1}\left( t \right)} \right) - 2b{x_6}\left( t \right) - {b^2}{x_5}\left( t \right) \hfill \\ \end{gathered} \right. $ |
式中C1、C2、C3和C4表示錐體細胞群與興奮性和抑制性中間神經元群間的耦合強度參數。神經群模型的輸出為y=x3-x5, 用于模擬輸出的腦電數據。
2.2 丘腦皮層神經群模型
研究表明，大腦的認知功能與皮層區域間的相互作用密切相關。對區域間連通性的理解有助于我們了解大腦的工作機制，評估在特定認知功能方面不同皮層區域的作用。大腦中并不是只有一個神經群單獨工作，而是若干個神經群落之間相互影響相互耦合地工作[15, 21]。多個神經群模型之間的耦合有單向或雙向連接，根據連接距離分為同一區域內的短程連接(short range connections, SRCs)與不同區域間的長程連接(long range connections, LRCs)。在經典神經群模型的基礎上演化而來的丘腦皮層神經群模型通過SRCs與LRCs實現多區域神經群落間的連接，從而模擬腦區間的耦合與神經振蕩。
對于SRCs，同一區域的一個神經群模型的PC、EI和Ⅱ的膜電位轉換到另一個神經群落的PC動作電位對應的脈沖響應分別用hd2(t)、hd3(t)和hd4(t)表示：
| $ {h_{d2}}\left( t \right) = A{a_{d2}}t{e^{ - {a_{d2}}t}} $ |
| $ {h_{d3}}\left( t \right) = A{a_{d3}}t{e^{ - {a_{d3}}t}} $ |
| $ {h_{d4}}\left( t \right) = B{b_{d4}}t{e^{ - {b_{d4}}t}} $ |
對于LRCs，不同區域間只有PC對PC的連接，模型的延遲對應的脈沖響應用hd1(t)表示：
| $ {h_{d1}}\left( t \right) = A{a_{d1}}t{e^{ - {a_{d1}}t}} $ |
常數ad1、ad2、ad3和bd4表示在皮層區域中，神經群之間傳出連接的平均時延。每個神經群PC的輸出都可以用來模擬腦電數據。在連續系統中，神經群模型可以用Runge-Kutta微分方法進行求解。而在離散系統中，神經群模型則采用歐拉法將其離散化后進行求解。
丘腦皮層神經群模型由大腦皮層神經群模型和丘腦神經群模型相互連接而成。大腦皮層神經群模型主要由PC、EI和Ⅱ三種不同的神經元群構成，如圖 2(b)所示。大腦皮層神經群模型接受皮層其他區域興奮性輸入，中間神經元對PC起反饋調節作用，模型間的交互作用通過C1到C9九個連接參數表示。其中C1與C2、C3與C4分別表示PC與EI、PC與Ⅱ間的平均耦合強度參數，C5表示PC與PC間的平均連接參數，C6、C7、C8和C9表示調節同一皮層區域內部以及該區域與其它皮層區域連接性的耦合強度參數。丘腦神經群模型由兩個主要的神經元群組成，它們分別是：丘腦皮層延遲神經元群(thalamo-cortical relay neurons, TCR)和丘腦抑制網狀核神經元群(thalamic reticular nucleus, TRN), 如圖 2(c)所示。在丘腦模型中，TCR主要是將感覺信息向大腦皮層進行投射，同時它會接受來自外界的感覺性輸入以及皮層模塊的信息，而TRN接受來自其他網狀核神經元的抑制性輸入。C1t和C2t反映的是TCR和TRN之間的耦合強度參數，C3t、C4t和C5t反映的是TCR分別與大腦皮層區域PC、EI和Ⅱ之間的耦合強度參數。Kd表示神經元之間的連接強度，它的大小與神經元之間的距離有關，距離近的連接強度較強，距離遠的連接相對較弱。它們參數的大小是從真實的擴散加權磁共振成像(diffusion-weighted magnetic resonance imaging, DWMRI)數據估算出來的。P、Pth分別代表大腦皮層與丘腦內部的基本隨機活動，它們模擬了神經群內部特定的振蕩活動，并且服從高斯分布。
大腦皮層神經群模型可由以下微分方程表示：
| $ \left\{ \begin{gathered} {{\dot x}_1}\left( t \right) = {x_2}\left( t \right) \hfill \\ {{\dot x}_2}\left( t \right) = Aa\left[ {{C_5}S\left( {{x_1}\left( t \right) - {x_3}\left( t \right)} \right) + {K_{d3}}{C_8}S\left( {{x_{13}}\left( t \right)} \right) + {K_{d1}}{C_6}S\left( {{x_9}\left( t \right)} \right) + {K_{d2}}{C_7}S\left( {{x_{11}}\left( t \right)} \right) + } \right. \hfill \\ \left. {\;\;\;\;\;\;\;\;\;\;\;{K_{d1t}}{C_{3t}}S\left( {{y_7}\left( t \right)} \right) + {C_2}S\left( {{x_5}\left( t \right)} \right) + P} \right] - 2a{x_2}\left( t \right) - {a^2}{x_1}\left( t \right) \hfill \\ {{\dot x}_3}\left( t \right) = {x_4}\left( t \right) \hfill \\ {{\dot x}_4}\left( t \right) = Bb\left[ {{C_4}S\left( {{x_7}\left( t \right)} \right) + {K_{d4}}{C_9}S\left( {{x_{15}}\left( t \right)} \right)} \right] - 2b{x_4}\left( t \right) - {b^2}{x_3}\left( t \right) \hfill \\ {{\dot x}_5}\left( t \right) = {x_6}\left( t \right) \hfill \\ {{\dot x}_6}\left( t \right) = Aa\left[ {{C_1}S\left( {{x_1}\left( t \right) - {x_3}\left( t \right)} \right) + {K_{d2t}}{C_9}S\left( {{y_9}\left( t \right)} \right)} \right] - 2a{x_6}\left( t \right) - {a^2}{x_5}\left( t \right) \hfill \\ {{\dot x}_7}\left( t \right) = {x_8}\left( t \right) \hfill \\ {{\dot x}_8}\left( t \right) = Aa\left[ {{C_3}S\left( {{x_1}\left( t \right) - {x_3}\left( t \right)} \right) + {K_{d3t}}{C_{5t}}S\left( {{y_{11}}\left( t \right)} \right)} \right] - 2a{x_8}\left( t \right) - {a^2}{x_7}\left( t \right) \hfill \\ {{\dot x}_9}\left( t \right) = {x_{10}}\left( t \right) \hfill \\ {{\dot x}_{10}}\left( t \right) = A{a_{d1}}S\left( {{x_1}\left( t \right) - {x_3}\left( t \right)} \right) - 2{a_{d1}}{x_{10}}\left( t \right) - a_{d1}^2{x_9}\left( t \right) \hfill \\ {{\dot x}_{11}}\left( t \right) = {x_{12}}\left( t \right) \hfill \\ {{\dot x}_{12}}\left( t \right) = A{a_{d2}}S\left( {{x_1}\left( t \right) - {x_3}\left( t \right)} \right) - 2{a_{d2}}{x_{12}}\left( t \right) - a_{d2}^2{x_{11}}\left( t \right) \hfill \\ {{\dot x}_{13}}\left( t \right) = {x_{14}}\left( t \right) \hfill \\ {{\dot x}_{14}}\left( t \right) = A{a_{d3}}S\left( {{x_5}\left( t \right)} \right) - 2{a_{d3}}{x_{14}}\left( t \right) - a_{d3}^2{x_{13}}\left( t \right) \hfill \\ {{\dot x}_{15}}\left( t \right) = {x_{16}}\left( t \right) \hfill \\ {{\dot x}_{16}}\left( t \right) = B{b_{d4}}S\left( {{x_7}\left( t \right)} \right) - 2{b_{d4}}{x_{16}}\left( t \right) - b_{d4}^2{x_{15}}\left( t \right) \hfill \\ \end{gathered} \right. $ |
丘腦神經群模型用微分方程表示如下：
| $ \left\{ \begin{gathered} {{\dot y}_1}\left( t \right) = {y_2}\left( t \right) \hfill \\ {{\dot y}_2}\left( t \right) = {A_t}{a_t}\left[ {{P_{th}} + {K_{d1t}}{C_6}S\left( {{x_9}\left( t \right)} \right)} \right] - 2{a_t}{y_2}\left( t \right) - a_t^2{y_1}\left( t \right) \hfill \\ {{\dot y}_3}\left( t \right) = {y_4}\left( t \right) \hfill \\ {{\dot y}_4}\left( t \right) = {B_t}{b_t}{C_{2t}}S\left( {{y_5}\left( t \right)} \right) - 2{b_t}{y_4}\left( t \right) - b_t^2{y_3}\left( t \right) \hfill \\ {{\dot y}_5}\left( t \right) = {y_6}\left( t \right) \hfill \\ {{\dot y}_6}\left( t \right) = {A_t}{a_t}{C_{1t}}S\left( {{y_1}\left( t \right) - {y_3}\left( t \right)} \right) - 2{a_t}{y_6}\left( t \right) - a_t^2{y_5}\left( t \right) \hfill \\ {{\dot y}_7}\left( t \right) = {y_8}\left( t \right) \hfill \\ {{\dot y}_8}\left( t \right) = {A_t}{a_{d1t}}S\left( {{y_1}\left( t \right) - {y_3}\left( t \right)} \right) - 2{a_{d1t}}{y_8}\left( t \right) - a_{d1t}^2{y_7}\left( t \right) \hfill \\ {{\dot y}_9}\left( t \right) = {y_{10}}\left( t \right) \hfill \\ {{\dot y}_{10}}\left( t \right) = {A_t}{a_{d2t}}S\left( {{y_1}\left( t \right) - {y_3}\left( t \right)} \right) - 2{a_{d2t}}{y_{10}}\left( t \right) - a_{d2t}^2{y_9}\left( t \right) \hfill \\ {{\dot y}_{11}}\left( t \right) = {y_{12}}\left( t \right) \hfill \\ {{\dot y}_{12}}\left( t \right) = {A_t}{a_{d3t}}S\left( {{y_1}\left( t \right) - {y_3}\left( t \right)} \right) - 2{a_{d3t}}{y_{12}}\left( t \right) - a_{d3t}^2{y_{11}}\left( t \right) \hfill \\ \end{gathered} \right. $ |
根據解剖學信息可以對模型中的參數進行估計，表 1列出了參數的生理學意義以及產生delta、theta、alpha、beta和gamma節律波形的典型參數取值。通過調整丘腦皮層神經群模型中的興奮性和抑制性參數以及改變興奮性和抑制性細胞群之間的耦合關系，丘腦皮層神經群模型可以產生多種節律腦電數據。
表1
模型參數的生理學意義及其典型值
Table1.
Physiological interpretation and typical values of model parameters
3 神經群模型的FPGA實現
硬件計算的思想是將神經群的數學模型轉換成FPGA可執行的硬件描述語言，利用FPGA實現神經群的功能再現。神經群模型一般為高階微分方程組，為了實現神經群模型在DSP Builder環境下的數字仿真，我們首先要將神經群模型微分方程組進行離散化，建立差分方程表示的數學模型。常見的微分方程離散化方法有龍格庫塔法、歐拉法等。由于龍格庫塔法需要復雜的運算過程，應用FPGA實現必然消耗大量的硬件資源。因此，本文采用定步長的歐拉法對微分方程進行離散化[22]。
式(3)經歐拉方法離散化后得到：
| $ \left\{ \begin{gathered} \frac{{z\left( {n + 1} \right) - z\left( n \right)}}{{dt}} = {z_1}\left( t \right) \hfill \\ \frac{{{z_1}\left( {n + 1} \right) - {z_1}\left( n \right)}}{{dt}} = Ggx\left( n \right) - 2g{z_1}\left( n \right) - {g^2}z\left( n \right) \hfill \\ \end{gathered} \right. $ |
簡化后如下：
| $ \left\{ \begin{gathered} z\left( {n + 1} \right) = {z_1}\left( t \right)*dt + z\left( n \right) \hfill \\ {z_1}\left( {n + 1} \right) = Ggx\left( n \right)*dt\left( {1 - 2g*dt} \right){z_1}\left( n \right) - {g^2}z\left( n \right)*dt \hfill \\ \end{gathered} \right. $ |
在DSP Builder中構建微分方程，如圖 3所示。
圖3
微分方程組電路框圖及其封裝模塊
Figure3.
Circuit diagram of differential equations and the package modules
將大腦皮層神經群模型離散化后得到其離散化方程組如下：
| $ \left\{ \begin{gathered} \frac{{{x_1}\left( {n + 1} \right) - {x_1}\left( n \right)}}{{dt}} = {x_2}\left( n \right) \hfill \\ \frac{{{x_2}\left( {n + 1} \right) - {x_2}\left( n \right)}}{{dt}} = Aa\left[ {{C_5}S\left( {{x_1}\left( n \right) - {x_3}\left( n \right)} \right) + {K_{d3}}{C_8}S\left( {{x_{13}}\left( n \right)} \right) + {K_{d1}}{C_6}S\left( {{x_9}\left( n \right)} \right) + {K_{d2}}{C_7}S\left( {{x_{11}}\left( n \right)} \right) + } \right. \hfill \\ \left. {\;\;\;\;\;\;\;\;\;\;\;{K_{d1t}}{C_{3t}}S\left( {{y_7}\left( n \right)} \right) + {C_2}S\left( {{x_5}\left( n \right)} \right) + P} \right] - 2a{x_2}\left( n \right) - {a^2}{x_1}\left( n \right) \hfill \\ \frac{{{x_3}\left( {n + 1} \right) - {x_3}\left( n \right)}}{{dt}} = {x_4}\left( n \right) \hfill \\ \frac{{{x_4}\left( {n + 1} \right) - {x_4}\left( n \right)}}{{dt}} = Bb\left[ {{C_4}S\left( {{x_7}\left( n \right)} \right) + {K_{d4}}{C_9}S\left( {{x_{15}}\left( n \right)} \right)} \right] - 2b{x_4}\left( n \right) - {b^2}{x_3}\left( n \right) \hfill \\ \frac{{{x_5}\left( {n + 1} \right) - {x_5}\left( n \right)}}{{dt}} = {x_6}\left( n \right) \hfill \\ \frac{{{x_6}\left( {n + 1} \right) - {x_6}\left( n \right)}}{{dt}} = Aa\left[ {{C_1}S\left( {{x_1}\left( n \right) - {x_3}\left( n \right)} \right) + {K_{d2t}}{C_9}S\left( {{y_9}\left( n \right)} \right)} \right] - 2a{x_6}\left( n \right) - {a^2}{x_5}\left( n \right) \hfill \\ \frac{{{x_7}\left( {n + 1} \right) - {x_7}\left( n \right)}}{{dt}} = {x_8}\left( n \right) \hfill \\ \frac{{{x_8}\left( {n + 1} \right) - {x_8}\left( n \right)}}{{dt}} = Aa\left[ {{C_3}S\left( {{x_1}\left( n \right) - {x_3}\left( n \right)} \right) + {K_{d3t}}{C_{5t}}S\left( {{y_{11}}\left( n \right)} \right)} \right] - 2a{x_8}\left( n \right) - {a^2}{x_7}\left( n \right) \hfill \\ \frac{{{x_9}\left( {n + 1} \right) - {x_9}\left( n \right)}}{{dt}} = {x_{10}}\left( n \right) \hfill \\ \frac{{{x_{10}}\left( {n + 1} \right) - {x_{10}}\left( n \right)}}{{dt}} = A{a_{d1}}S\left( {{x_1}\left( n \right) - {x_3}\left( n \right)} \right) - 2{a_{d1}}{x_{10}}\left( n \right) - a_{d1}^2{x_9}\left( n \right) \hfill \\ \frac{{{x_{11}}\left( {n + 1} \right) - {x_{11}}\left( n \right)}}{{dt}} = {x_{12}}\left( n \right) \hfill \\ \frac{{{x_{12}}\left( {n + 1} \right) - {x_{12}}\left( n \right)}}{{dt}} = A{a_{d2}}S\left( {{x_1}\left( n \right) - {x_3}\left( n \right)} \right) - 2{a_{d2}}{x_{12}}\left( n \right) - a_{d2}^2{x_{11}}\left( n \right) \hfill \\ \frac{{{x_{13}}\left( {n + 1} \right) - {x_{13}}\left( n \right)}}{{dt}} = {x_{14}}\left( n \right) \hfill \\ \frac{{{x_{14}}\left( {n + 1} \right) - {x_{14}}\left( n \right)}}{{dt}} = A{a_{d3}}S\left( {{x_5}\left( n \right)} \right) - 2{a_{d3}}{x_{14}}\left( n \right) - a_{d3}^2{x_{13}}\left( n \right) \hfill \\ \frac{{{x_{15}}\left( {n + 1} \right) - {x_{15}}\left( n \right)}}{{dt}} = {x_{16}}\left( n \right) \hfill \\ \frac{{{x_{16}}\left( {n + 1} \right) - {x_{16}}\left( n \right)}}{{dt}} = B{b_{d4}}S\left( {{x_7}\left( n \right)} \right) - 2{b_{d4}}{x_{16}}\left( n \right) - b_{d4}^2{x_{15}}\left( n \right) \hfill \\ \end{gathered} \right. $ |
丘腦神經群模型離散化后得到其離散化方程組如下：
| $ \left\{ \begin{gathered} \frac{{{y_1}\left( {n + 1} \right) - {y_1}\left( n \right)}}{{dt}} = {y_2}\left( n \right) \hfill \\ \frac{{{y_2}\left( {n + 1} \right) - {y_2}\left( n \right)}}{{dt}} = {A_t}{a_t}\left[ {{P_{th}} + {K_{d1t}}{C_6}S\left( {{x_9}\left( n \right)} \right)} \right] - 2{a_t}{y_2}\left( n \right) - a_t^2{y_1}\left( n \right) \hfill \\ \frac{{{y_3}\left( {n + 1} \right) - {y_3}\left( n \right)}}{{dt}} = {y_4}\left( n \right) \hfill \\ \frac{{{y_4}\left( {n + 1} \right) - {y_4}\left( n \right)}}{{dt}} = {B_t}{b_t}{C_{2t}}S\left( {{y_5}\left( n \right)} \right) - 2{b_t}{y_4}\left( n \right) - b_t^2{y_3}\left( n \right) \hfill \\ \frac{{{y_5}\left( {n + 1} \right) - {y_5}\left( n \right)}}{{dt}} = {y_6}\left( n \right) \hfill \\ \frac{{{y_6}\left( {n + 1} \right) - {y_6}\left( n \right)}}{{dt}} = {A_t}{a_t}{C_{1t}}S\left( {{y_1}\left( n \right) - {y_3}\left( n \right)} \right) - 2{a_t}{y_6}\left( n \right) - a_t^2{y_5}\left( n \right) \hfill \\ \frac{{{y_7}\left( {n + 1} \right) - {y_7}\left( n \right)}}{{dt}} = {y_8}\left( n \right) \hfill \\ \frac{{{y_8}\left( {n + 1} \right) - {y_8}\left( n \right)}}{{dt}} = {A_t}{a_{d1t}}S\left( {{y_1}\left( n \right) - {y_3}\left( n \right)} \right) - 2{a_{d1t}}{y_8}\left( n \right) - a_{d1t}^2{y_7}\left( n \right) \hfill \\ \frac{{{y_9}\left( {n + 1} \right) - {y_9}\left( n \right)}}{{dt}} = {y_{10}}\left( n \right) \hfill \\ \frac{{{y_{10}}\left( {n + 1} \right) - {y_{10}}\left( n \right)}}{{dt}} = {A_t}{a_{d2t}}S\left( {{y_1}\left( n \right) - {y_3}\left( n \right)} \right) - 2{a_{d2t}}{y_{10}}\left( n \right) - a_{d2t}^2{y_9}\left( n \right) \hfill \\ \frac{{{y_{11}}\left( {n + 1} \right) - {y_{11}}\left( n \right)}}{{dt}} = {y_{12}}\left( n \right) \hfill \\ \frac{{{y_{12}}\left( {n + 1} \right) - {y_{12}}\left( n \right)}}{{dt}} = {A_t}{a_{d3t}}S\left( {{y_1}\left( n \right) - {y_3}\left( n \right)} \right) - 2{a_{d3t}}{y_{12}}\left( n \right) - a_{d3t}^2{y_{11}}\left( n \right) \hfill \\ \end{gathered} \right. $ |
在DSP Builder環境下，搭建的大腦皮層神經群模型及丘腦神經群模型如圖 4所示。
圖4
在DSP Builder中構建的大腦皮層神經群模型及丘腦神經群模型
(a)DSP Builder環境下的大腦皮層神經群模型; (b)DSP Builder環境下的丘腦神經群模型
Figure4. Neural mass model of the cortex and thalamus in DSP Builder(a) the NMM of a cortical voxel in DSP Builder; (b) the NMM of the thalamus in DSP Builder
仿真實驗中，采用歐拉法求解上述差分方程組，dt取0.001。根據表 1中的典型參數仿真出的δ、θ、α、β和γ波形如圖 5(b)所示。
圖5
平均突觸時間常數[a，b]對應的腦電各個節律信號
(a)平均突觸時間常數[
(a) areal distribution of the average synaptic time constants for
4 不同振蕩信號的參數選擇
表 1中列出了產生腦電各個節律振蕩信號的典型參數值，通過對比可以發現，在其他參數不變的情況下，增大C1的值，腦電的振蕩活動從beta向gamma頻段移動。同樣，保持C1不變，增大興奮性或抑制性的時間常數，振蕩活動從delta向theta頻段移動。從中可以發現，通過調整神經群的動態性能參數，可以產生從delta到gamma頻段的振蕩活動信號。然而，神經群模型中有很多參數，單獨研究每個參數的改變對模型輸出的影響是困難的。為此，我們通過改變興奮性時間常數a和抑制性時間常數b獲得產生相同腦電振蕩活動的參數區域。為了保證相應的振蕩活動，A/a、B/b的比值應為常數。對于每組a和b對應的振蕩活動信號，我們首先應該保證該信號為類正弦信號，進而求出該信號頻譜的峰值，將其分配到相應的腦電節律頻段區域。本文采用表 1中產生beta振蕩活動的典型參數值，參數a、b分別從5 s-1到200 s-1逐步增加，步進為5 s-1，其他參數保持不變，仿真結果如圖 5(a)所示。由圖 5可知，通過改變相應的參數，丘腦皮層神經群模型輸出的信號頻帶范圍會有相應的變化。圖 5(a)統計了腦電各個節律對應參數的變化范圍，如果要獲得多個節律的腦電波形，可以取相鄰腦電節律區域間的參數值，這樣輸出信號的頻帶會相對較寬，即包含多個節律的腦電信號。另一方面，經典神經群模型屬于單動力學神經群模型，單動力學神經群模型能夠得到不同頻率的窄帶信號，將單動力學神經群模型通過一定的方式連接起來，可以構建頻帶更寬的多動力學神經群模型，由此也可以得到含有多個節律的腦電信號。
5 結論
本文研究了丘腦皮層神經群模型，并采用FPGA實現了離散化快速計算和仿真，通過調整參數變化，給出了丘腦皮層神經群模型不同參數下對應的腦電節律信息。本文所建模型可以用于模擬不同腦功能狀態下的神經振蕩，生成正常與癲癇下的腦電數據以及各個節律下的腦電數據，并能對幅值進行調節。設計采用了MATLAB的Simulink與FPGA的DSP Builder相結合的方式，實現了軟件與硬件系統的有效融合，為神經群模型的計算機模擬提供了新的方案和思路。
引言
腦電是大腦神經活動的宏觀表現，在不同的腦功能活動和精神疾病狀態下會表現出特定的神經振蕩現象，例如睡眠、癲癇發作[1]、帕金森病[2]等都會表現出特有的神經振蕩活動[3-4]。因此，從神經振蕩角度研究腦電活動對于理解大腦的工作機制具有重要意義。尤其是采用數學和物理模型來模擬大腦的神經振蕩活動對于深入分析大腦的功能機制具有重要的價值[5-6]。
目前，針對大腦的神經活動，從微觀、介觀、宏觀等不同尺度均有建模。其中，微觀模型主要以單個神經元為單位進行建模[7]。從微觀層面對神經元的建模有助于了解大腦最基本的組成單元的工作原理，但是大腦由上千億個不同種類的神經細胞構成，很難確定每個神經元模型的參數，且各種神經元之間的連接和大腦神經系統的結構也異常復雜，大腦的功能也是眾多神經元群協同作用的結果，完成這項工程需要消耗巨大的計算機資源, 僅從微觀層次上建模對整個大腦或者大腦某一區域的研究來說是遠遠不夠的。因此，微觀模型不能真實地反映宏觀的腦電特征，只能用來仿真單個神經元[8-9]。介觀模型是建立微觀和宏觀的橋梁，以集總參數模型為主，對特定種類細胞組成的神經元群整體特性進行建模，例如單一的神經群模型(neural mass model, NMM)和平均場模型(mean-field cortical, MFC)，NMM和MFC起初是用來仿真嗅覺皮層的活動[10]以及自發的alpha節律腦電波[11]，隨后這些模型經過改進和擴展，能夠用較少的狀態變量來描述大量神經元群的平均電活動[12-13]。而宏觀模型則考慮了神經元群之間的耦合關系和腦區間的耦合關系[14-15]，比如多通道神經元群模型[16]。Wendling等[12]利用三個神經群模型的不同耦合方式，產生類癲癇數據，與真實的癲癇數據對比，研究癲癇的發作機制。通過對多通道神經振蕩的建模及其同真實腦電數據之間的對比分析，可以幫助理解大腦區域間的耦合關系以及癲癇、帕金森病、阿爾茨海默癥等的發病機制，驗證各種腦電數據處理算法的有效性。在三種尺度模型的研究方面，以介觀層面的研究為多。尤其是在經典神經群模型的基礎上演化而來的丘腦皮層神經群模型(thalamocortical neural mass model, TNMM)，在睡眠、記憶力鞏固以及節律產生等腦功能認知方面有著重要的作用[17-18]。通過簡單的數學模型可以在介觀層面仿真出大腦皮層與丘腦間的連接機制，對腦功能的研究具有重大價值。
隨著對神經群模型計算仿真研究的不斷深入，神經群規模和復雜度不斷增加，當前通用計算機的計算無法滿足對大規模的復雜神經群模型實現實時仿真的需求。基于并行計算和分布式的計算是解決計算復雜度的重要方法。尤其以現場可編程門陣列(Field Programmable Gate Array，FPGA)為基礎的硬件計算方法在并行計算和求解復雜非線性常微分方程的數值解方面具有一定的優勢，而且易于和MATLAB形成無縫連接。因此，本文利用FPGA模擬復雜的丘腦皮層神經群數學模型，采用Altera公司提供的計算機輔助設計工具DSP Builder聯合Simulink完成丘腦皮層神經群模型的設計，使用FPGA內部定制的NiosⅡ嵌入式處理器完成對硬件系統中外部設備的驅動程序設計，并實現對核心算法模塊參數的控制。
1 FPGA仿真平臺設計
FPGA大都是基于查找表技術，通過燒寫文件改變查找表內容的方法來實現對FPGA的重復配置。FPGA芯片主要由可編程輸入輸出單元、基本可編程邏輯單元、內嵌SRAM、豐富的布線資源、底層嵌入功能單元和內嵌專用單元等組成，可以完成極其復雜的時序與邏輯組合電路功能。基于FPGA的數字建模仿真平臺如圖 1所示。該平臺采用Altera公司的Cyclone Ⅲ系列FPGA芯片作為系統硬件核心，硬件系統中外部設備的驅動程序主要由FPGA內部定制的NiosⅡ嵌入式處理器完成，通過軟件編程實現人機交互以及對模型數據收發的控制，最大程度地利用FPGA內部邏輯單元，使處理器的外部接口更少。這樣可以大大提高系統的性能和穩定性，同時降低研制成本。
圖1
FPGA數字建模仿真平臺結構示意圖
Figure1.
Digital model and simulation platform of FPGA
DSP Builder是一個系統級設計工具，以MATLAB/Simulink的模塊庫形式出現，在Simulink中進行圖形化的設計和仿真，同時把算法仿真建模和硬件實現的設計工具連接起來，使算法的開發到硬件的實現可以無縫地過渡。神經群算法模型主要在DSP Builder中搭建完成，DSP Builder提供了一個從MATLAB/Simulink直接到FPGA硬件實現的設計接口，在Simulink中完成系統集成，然后通過Signal Compiler模塊將電路模型文件即Simulink模塊文件(.mdl)轉換成寄存器傳輸級(register transfer level, RTL)的VHDL代碼表述和工具命令語言(Tool Command Language, TCL)腳本，使其可以在FPGA中完成硬件電路的實現。同時通過QuartusⅡ中的SOPC Builder工具在FPGA內部定制NiosⅡ嵌入式處理器，將在DSP Builder中生成的模型IP核通過Avalon總線與NiosⅡ處理器相連接，從而將軟件的靈活性與硬件的高速性結合起來，完成與其它模塊的數據通信。
2 模型的建立
2.1 經典神經群模型
經典神經群模型[19-20]的基本思想是使興奮性和抑制性細胞群相互作用，產生相應的神經振蕩活動。單個神經群模型主要由錐體細胞(pyramidal cells, PC)、興奮性中間神經元(excitatory interneurons, EI)和抑制性中間神經元(inhibitory interneurons, Ⅱ)三種不同的神經群構成。經典神經群模型框圖如圖 2(a)所示。
圖2
丘腦皮層神經群模型構成原理框圖
(a)經典神經群模型; (b)大腦皮層神經群模型; (c)丘腦神經群模型
Figure2. Schematic diagram of TNMM(a) classical NMM; (b) the NMM of a cortical voxel; (c) the NMM of the thalamus
經典神經群的動力學演化主要包含兩部分。第一部分是將突觸前信息(動作電位的平均脈沖密度)轉變成突觸后信息(突觸后膜興奮或抑制電位)。
在興奮性和抑制性狀況下，線性傳遞函數的脈沖響應分別為：
| $ {h_e}\left( t \right) = Aat{e^{ - at}} $ |
| $ {h_i}\left( t \right) = Bat{e^{ - bt}} $ |
其中, 參數A和B分別為興奮性和抑制性平均突觸增益，a和b分別為興奮性和抑制性細胞平均突觸時間常數，通過調節這四個參數可以改變興奮和抑制突觸的靈敏度，確定模型的興奮和抑制特性。線性函數hi(t)和he(t)可由如下形式的微分方程描述：
| $ \left\{ \begin{gathered} \dot z\left( t \right) = {z_1}\left( t \right) \hfill \\ \dot z\left( t \right) = Ggx\left( t \right) - 2g{z_1}\left( t \right) - {g^2}z\left( t \right) \hfill \\ \end{gathered} \right. $ |
其中，G和g表示興奮性條件下的參數A和a或者抑制條件下的參數B和b，x(t)和z(t)分別為神經群模型的輸入和輸出數據。
第二部分是通過靜態非線性函數將平均膜電壓轉換成動作電位的平均密度，即平均點燃率，并作為線性變換函數的輸入。靜態非線性函數表示為：
| $ S\left( v \right) = 2{e_0}/\left( {1 + {e^{r\left( {{v_0} - v} \right)}}} \right) $ |
其中2e0表示最大點燃率，v0表示相對于點燃率e0的突觸后電位，r表示S型函數的彎曲程度，v表示突觸前平均膜電壓。
經典神經群模型可由以下微分方程表示：
| $ \left\{ \begin{gathered} {{\dot x}_1}\left( t \right) = {x_2}\left( t \right) \hfill \\ {{\dot x}_2}\left( t \right) = AaS\left( {{x_3}\left( t \right) - {x_5}\left( t \right)} \right) - 2a{x_2}\left( t \right) - {a^2}{x_1}\left( t \right) \hfill \\ {{\dot x}_3}\left( t \right) = {x_4}\left( t \right) \hfill \\ {{\dot x}_4}\left( t \right) = Aa\left[ {p + {C_2}S\left( {{C_1}{x_1}\left( t \right)} \right)} \right] - 2a{x_4}\left( t \right) - {a^2}{x_3}\left( t \right) \hfill \\ {{\dot x}_5}\left( t \right) = {x_6}\left( t \right) \hfill \\ {{\dot x}_6}\left( t \right) = Bb{C_4}S\left( {{C_3}{x_1}\left( t \right)} \right) - 2b{x_6}\left( t \right) - {b^2}{x_5}\left( t \right) \hfill \\ \end{gathered} \right. $ |
式中C1、C2、C3和C4表示錐體細胞群與興奮性和抑制性中間神經元群間的耦合強度參數。神經群模型的輸出為y=x3-x5, 用于模擬輸出的腦電數據。
2.2 丘腦皮層神經群模型
研究表明，大腦的認知功能與皮層區域間的相互作用密切相關。對區域間連通性的理解有助于我們了解大腦的工作機制，評估在特定認知功能方面不同皮層區域的作用。大腦中并不是只有一個神經群單獨工作，而是若干個神經群落之間相互影響相互耦合地工作[15, 21]。多個神經群模型之間的耦合有單向或雙向連接，根據連接距離分為同一區域內的短程連接(short range connections, SRCs)與不同區域間的長程連接(long range connections, LRCs)。在經典神經群模型的基礎上演化而來的丘腦皮層神經群模型通過SRCs與LRCs實現多區域神經群落間的連接，從而模擬腦區間的耦合與神經振蕩。
對于SRCs，同一區域的一個神經群模型的PC、EI和Ⅱ的膜電位轉換到另一個神經群落的PC動作電位對應的脈沖響應分別用hd2(t)、hd3(t)和hd4(t)表示：
| $ {h_{d2}}\left( t \right) = A{a_{d2}}t{e^{ - {a_{d2}}t}} $ |
| $ {h_{d3}}\left( t \right) = A{a_{d3}}t{e^{ - {a_{d3}}t}} $ |
| $ {h_{d4}}\left( t \right) = B{b_{d4}}t{e^{ - {b_{d4}}t}} $ |
對于LRCs，不同區域間只有PC對PC的連接，模型的延遲對應的脈沖響應用hd1(t)表示：
| $ {h_{d1}}\left( t \right) = A{a_{d1}}t{e^{ - {a_{d1}}t}} $ |
常數ad1、ad2、ad3和bd4表示在皮層區域中，神經群之間傳出連接的平均時延。每個神經群PC的輸出都可以用來模擬腦電數據。在連續系統中，神經群模型可以用Runge-Kutta微分方法進行求解。而在離散系統中，神經群模型則采用歐拉法將其離散化后進行求解。
丘腦皮層神經群模型由大腦皮層神經群模型和丘腦神經群模型相互連接而成。大腦皮層神經群模型主要由PC、EI和Ⅱ三種不同的神經元群構成，如圖 2(b)所示。大腦皮層神經群模型接受皮層其他區域興奮性輸入，中間神經元對PC起反饋調節作用，模型間的交互作用通過C1到C9九個連接參數表示。其中C1與C2、C3與C4分別表示PC與EI、PC與Ⅱ間的平均耦合強度參數，C5表示PC與PC間的平均連接參數，C6、C7、C8和C9表示調節同一皮層區域內部以及該區域與其它皮層區域連接性的耦合強度參數。丘腦神經群模型由兩個主要的神經元群組成，它們分別是：丘腦皮層延遲神經元群(thalamo-cortical relay neurons, TCR)和丘腦抑制網狀核神經元群(thalamic reticular nucleus, TRN), 如圖 2(c)所示。在丘腦模型中，TCR主要是將感覺信息向大腦皮層進行投射，同時它會接受來自外界的感覺性輸入以及皮層模塊的信息，而TRN接受來自其他網狀核神經元的抑制性輸入。C1t和C2t反映的是TCR和TRN之間的耦合強度參數，C3t、C4t和C5t反映的是TCR分別與大腦皮層區域PC、EI和Ⅱ之間的耦合強度參數。Kd表示神經元之間的連接強度，它的大小與神經元之間的距離有關，距離近的連接強度較強，距離遠的連接相對較弱。它們參數的大小是從真實的擴散加權磁共振成像(diffusion-weighted magnetic resonance imaging, DWMRI)數據估算出來的。P、Pth分別代表大腦皮層與丘腦內部的基本隨機活動，它們模擬了神經群內部特定的振蕩活動，并且服從高斯分布。
大腦皮層神經群模型可由以下微分方程表示：
| $ \left\{ \begin{gathered} {{\dot x}_1}\left( t \right) = {x_2}\left( t \right) \hfill \\ {{\dot x}_2}\left( t \right) = Aa\left[ {{C_5}S\left( {{x_1}\left( t \right) - {x_3}\left( t \right)} \right) + {K_{d3}}{C_8}S\left( {{x_{13}}\left( t \right)} \right) + {K_{d1}}{C_6}S\left( {{x_9}\left( t \right)} \right) + {K_{d2}}{C_7}S\left( {{x_{11}}\left( t \right)} \right) + } \right. \hfill \\ \left. {\;\;\;\;\;\;\;\;\;\;\;{K_{d1t}}{C_{3t}}S\left( {{y_7}\left( t \right)} \right) + {C_2}S\left( {{x_5}\left( t \right)} \right) + P} \right] - 2a{x_2}\left( t \right) - {a^2}{x_1}\left( t \right) \hfill \\ {{\dot x}_3}\left( t \right) = {x_4}\left( t \right) \hfill \\ {{\dot x}_4}\left( t \right) = Bb\left[ {{C_4}S\left( {{x_7}\left( t \right)} \right) + {K_{d4}}{C_9}S\left( {{x_{15}}\left( t \right)} \right)} \right] - 2b{x_4}\left( t \right) - {b^2}{x_3}\left( t \right) \hfill \\ {{\dot x}_5}\left( t \right) = {x_6}\left( t \right) \hfill \\ {{\dot x}_6}\left( t \right) = Aa\left[ {{C_1}S\left( {{x_1}\left( t \right) - {x_3}\left( t \right)} \right) + {K_{d2t}}{C_9}S\left( {{y_9}\left( t \right)} \right)} \right] - 2a{x_6}\left( t \right) - {a^2}{x_5}\left( t \right) \hfill \\ {{\dot x}_7}\left( t \right) = {x_8}\left( t \right) \hfill \\ {{\dot x}_8}\left( t \right) = Aa\left[ {{C_3}S\left( {{x_1}\left( t \right) - {x_3}\left( t \right)} \right) + {K_{d3t}}{C_{5t}}S\left( {{y_{11}}\left( t \right)} \right)} \right] - 2a{x_8}\left( t \right) - {a^2}{x_7}\left( t \right) \hfill \\ {{\dot x}_9}\left( t \right) = {x_{10}}\left( t \right) \hfill \\ {{\dot x}_{10}}\left( t \right) = A{a_{d1}}S\left( {{x_1}\left( t \right) - {x_3}\left( t \right)} \right) - 2{a_{d1}}{x_{10}}\left( t \right) - a_{d1}^2{x_9}\left( t \right) \hfill \\ {{\dot x}_{11}}\left( t \right) = {x_{12}}\left( t \right) \hfill \\ {{\dot x}_{12}}\left( t \right) = A{a_{d2}}S\left( {{x_1}\left( t \right) - {x_3}\left( t \right)} \right) - 2{a_{d2}}{x_{12}}\left( t \right) - a_{d2}^2{x_{11}}\left( t \right) \hfill \\ {{\dot x}_{13}}\left( t \right) = {x_{14}}\left( t \right) \hfill \\ {{\dot x}_{14}}\left( t \right) = A{a_{d3}}S\left( {{x_5}\left( t \right)} \right) - 2{a_{d3}}{x_{14}}\left( t \right) - a_{d3}^2{x_{13}}\left( t \right) \hfill \\ {{\dot x}_{15}}\left( t \right) = {x_{16}}\left( t \right) \hfill \\ {{\dot x}_{16}}\left( t \right) = B{b_{d4}}S\left( {{x_7}\left( t \right)} \right) - 2{b_{d4}}{x_{16}}\left( t \right) - b_{d4}^2{x_{15}}\left( t \right) \hfill \\ \end{gathered} \right. $ |
丘腦神經群模型用微分方程表示如下：
| $ \left\{ \begin{gathered} {{\dot y}_1}\left( t \right) = {y_2}\left( t \right) \hfill \\ {{\dot y}_2}\left( t \right) = {A_t}{a_t}\left[ {{P_{th}} + {K_{d1t}}{C_6}S\left( {{x_9}\left( t \right)} \right)} \right] - 2{a_t}{y_2}\left( t \right) - a_t^2{y_1}\left( t \right) \hfill \\ {{\dot y}_3}\left( t \right) = {y_4}\left( t \right) \hfill \\ {{\dot y}_4}\left( t \right) = {B_t}{b_t}{C_{2t}}S\left( {{y_5}\left( t \right)} \right) - 2{b_t}{y_4}\left( t \right) - b_t^2{y_3}\left( t \right) \hfill \\ {{\dot y}_5}\left( t \right) = {y_6}\left( t \right) \hfill \\ {{\dot y}_6}\left( t \right) = {A_t}{a_t}{C_{1t}}S\left( {{y_1}\left( t \right) - {y_3}\left( t \right)} \right) - 2{a_t}{y_6}\left( t \right) - a_t^2{y_5}\left( t \right) \hfill \\ {{\dot y}_7}\left( t \right) = {y_8}\left( t \right) \hfill \\ {{\dot y}_8}\left( t \right) = {A_t}{a_{d1t}}S\left( {{y_1}\left( t \right) - {y_3}\left( t \right)} \right) - 2{a_{d1t}}{y_8}\left( t \right) - a_{d1t}^2{y_7}\left( t \right) \hfill \\ {{\dot y}_9}\left( t \right) = {y_{10}}\left( t \right) \hfill \\ {{\dot y}_{10}}\left( t \right) = {A_t}{a_{d2t}}S\left( {{y_1}\left( t \right) - {y_3}\left( t \right)} \right) - 2{a_{d2t}}{y_{10}}\left( t \right) - a_{d2t}^2{y_9}\left( t \right) \hfill \\ {{\dot y}_{11}}\left( t \right) = {y_{12}}\left( t \right) \hfill \\ {{\dot y}_{12}}\left( t \right) = {A_t}{a_{d3t}}S\left( {{y_1}\left( t \right) - {y_3}\left( t \right)} \right) - 2{a_{d3t}}{y_{12}}\left( t \right) - a_{d3t}^2{y_{11}}\left( t \right) \hfill \\ \end{gathered} \right. $ |
根據解剖學信息可以對模型中的參數進行估計，表 1列出了參數的生理學意義以及產生delta、theta、alpha、beta和gamma節律波形的典型參數取值。通過調整丘腦皮層神經群模型中的興奮性和抑制性參數以及改變興奮性和抑制性細胞群之間的耦合關系，丘腦皮層神經群模型可以產生多種節律腦電數據。
表1
模型參數的生理學意義及其典型值
Table1.
Physiological interpretation and typical values of model parameters
3 神經群模型的FPGA實現
硬件計算的思想是將神經群的數學模型轉換成FPGA可執行的硬件描述語言，利用FPGA實現神經群的功能再現。神經群模型一般為高階微分方程組，為了實現神經群模型在DSP Builder環境下的數字仿真，我們首先要將神經群模型微分方程組進行離散化，建立差分方程表示的數學模型。常見的微分方程離散化方法有龍格庫塔法、歐拉法等。由于龍格庫塔法需要復雜的運算過程，應用FPGA實現必然消耗大量的硬件資源。因此，本文采用定步長的歐拉法對微分方程進行離散化[22]。
式(3)經歐拉方法離散化后得到：
| $ \left\{ \begin{gathered} \frac{{z\left( {n + 1} \right) - z\left( n \right)}}{{dt}} = {z_1}\left( t \right) \hfill \\ \frac{{{z_1}\left( {n + 1} \right) - {z_1}\left( n \right)}}{{dt}} = Ggx\left( n \right) - 2g{z_1}\left( n \right) - {g^2}z\left( n \right) \hfill \\ \end{gathered} \right. $ |
簡化后如下：
| $ \left\{ \begin{gathered} z\left( {n + 1} \right) = {z_1}\left( t \right)*dt + z\left( n \right) \hfill \\ {z_1}\left( {n + 1} \right) = Ggx\left( n \right)*dt\left( {1 - 2g*dt} \right){z_1}\left( n \right) - {g^2}z\left( n \right)*dt \hfill \\ \end{gathered} \right. $ |
在DSP Builder中構建微分方程，如圖 3所示。
圖3
微分方程組電路框圖及其封裝模塊
Figure3.
Circuit diagram of differential equations and the package modules
將大腦皮層神經群模型離散化后得到其離散化方程組如下：
| $ \left\{ \begin{gathered} \frac{{{x_1}\left( {n + 1} \right) - {x_1}\left( n \right)}}{{dt}} = {x_2}\left( n \right) \hfill \\ \frac{{{x_2}\left( {n + 1} \right) - {x_2}\left( n \right)}}{{dt}} = Aa\left[ {{C_5}S\left( {{x_1}\left( n \right) - {x_3}\left( n \right)} \right) + {K_{d3}}{C_8}S\left( {{x_{13}}\left( n \right)} \right) + {K_{d1}}{C_6}S\left( {{x_9}\left( n \right)} \right) + {K_{d2}}{C_7}S\left( {{x_{11}}\left( n \right)} \right) + } \right. \hfill \\ \left. {\;\;\;\;\;\;\;\;\;\;\;{K_{d1t}}{C_{3t}}S\left( {{y_7}\left( n \right)} \right) + {C_2}S\left( {{x_5}\left( n \right)} \right) + P} \right] - 2a{x_2}\left( n \right) - {a^2}{x_1}\left( n \right) \hfill \\ \frac{{{x_3}\left( {n + 1} \right) - {x_3}\left( n \right)}}{{dt}} = {x_4}\left( n \right) \hfill \\ \frac{{{x_4}\left( {n + 1} \right) - {x_4}\left( n \right)}}{{dt}} = Bb\left[ {{C_4}S\left( {{x_7}\left( n \right)} \right) + {K_{d4}}{C_9}S\left( {{x_{15}}\left( n \right)} \right)} \right] - 2b{x_4}\left( n \right) - {b^2}{x_3}\left( n \right) \hfill \\ \frac{{{x_5}\left( {n + 1} \right) - {x_5}\left( n \right)}}{{dt}} = {x_6}\left( n \right) \hfill \\ \frac{{{x_6}\left( {n + 1} \right) - {x_6}\left( n \right)}}{{dt}} = Aa\left[ {{C_1}S\left( {{x_1}\left( n \right) - {x_3}\left( n \right)} \right) + {K_{d2t}}{C_9}S\left( {{y_9}\left( n \right)} \right)} \right] - 2a{x_6}\left( n \right) - {a^2}{x_5}\left( n \right) \hfill \\ \frac{{{x_7}\left( {n + 1} \right) - {x_7}\left( n \right)}}{{dt}} = {x_8}\left( n \right) \hfill \\ \frac{{{x_8}\left( {n + 1} \right) - {x_8}\left( n \right)}}{{dt}} = Aa\left[ {{C_3}S\left( {{x_1}\left( n \right) - {x_3}\left( n \right)} \right) + {K_{d3t}}{C_{5t}}S\left( {{y_{11}}\left( n \right)} \right)} \right] - 2a{x_8}\left( n \right) - {a^2}{x_7}\left( n \right) \hfill \\ \frac{{{x_9}\left( {n + 1} \right) - {x_9}\left( n \right)}}{{dt}} = {x_{10}}\left( n \right) \hfill \\ \frac{{{x_{10}}\left( {n + 1} \right) - {x_{10}}\left( n \right)}}{{dt}} = A{a_{d1}}S\left( {{x_1}\left( n \right) - {x_3}\left( n \right)} \right) - 2{a_{d1}}{x_{10}}\left( n \right) - a_{d1}^2{x_9}\left( n \right) \hfill \\ \frac{{{x_{11}}\left( {n + 1} \right) - {x_{11}}\left( n \right)}}{{dt}} = {x_{12}}\left( n \right) \hfill \\ \frac{{{x_{12}}\left( {n + 1} \right) - {x_{12}}\left( n \right)}}{{dt}} = A{a_{d2}}S\left( {{x_1}\left( n \right) - {x_3}\left( n \right)} \right) - 2{a_{d2}}{x_{12}}\left( n \right) - a_{d2}^2{x_{11}}\left( n \right) \hfill \\ \frac{{{x_{13}}\left( {n + 1} \right) - {x_{13}}\left( n \right)}}{{dt}} = {x_{14}}\left( n \right) \hfill \\ \frac{{{x_{14}}\left( {n + 1} \right) - {x_{14}}\left( n \right)}}{{dt}} = A{a_{d3}}S\left( {{x_5}\left( n \right)} \right) - 2{a_{d3}}{x_{14}}\left( n \right) - a_{d3}^2{x_{13}}\left( n \right) \hfill \\ \frac{{{x_{15}}\left( {n + 1} \right) - {x_{15}}\left( n \right)}}{{dt}} = {x_{16}}\left( n \right) \hfill \\ \frac{{{x_{16}}\left( {n + 1} \right) - {x_{16}}\left( n \right)}}{{dt}} = B{b_{d4}}S\left( {{x_7}\left( n \right)} \right) - 2{b_{d4}}{x_{16}}\left( n \right) - b_{d4}^2{x_{15}}\left( n \right) \hfill \\ \end{gathered} \right. $ |
丘腦神經群模型離散化后得到其離散化方程組如下：
| $ \left\{ \begin{gathered} \frac{{{y_1}\left( {n + 1} \right) - {y_1}\left( n \right)}}{{dt}} = {y_2}\left( n \right) \hfill \\ \frac{{{y_2}\left( {n + 1} \right) - {y_2}\left( n \right)}}{{dt}} = {A_t}{a_t}\left[ {{P_{th}} + {K_{d1t}}{C_6}S\left( {{x_9}\left( n \right)} \right)} \right] - 2{a_t}{y_2}\left( n \right) - a_t^2{y_1}\left( n \right) \hfill \\ \frac{{{y_3}\left( {n + 1} \right) - {y_3}\left( n \right)}}{{dt}} = {y_4}\left( n \right) \hfill \\ \frac{{{y_4}\left( {n + 1} \right) - {y_4}\left( n \right)}}{{dt}} = {B_t}{b_t}{C_{2t}}S\left( {{y_5}\left( n \right)} \right) - 2{b_t}{y_4}\left( n \right) - b_t^2{y_3}\left( n \right) \hfill \\ \frac{{{y_5}\left( {n + 1} \right) - {y_5}\left( n \right)}}{{dt}} = {y_6}\left( n \right) \hfill \\ \frac{{{y_6}\left( {n + 1} \right) - {y_6}\left( n \right)}}{{dt}} = {A_t}{a_t}{C_{1t}}S\left( {{y_1}\left( n \right) - {y_3}\left( n \right)} \right) - 2{a_t}{y_6}\left( n \right) - a_t^2{y_5}\left( n \right) \hfill \\ \frac{{{y_7}\left( {n + 1} \right) - {y_7}\left( n \right)}}{{dt}} = {y_8}\left( n \right) \hfill \\ \frac{{{y_8}\left( {n + 1} \right) - {y_8}\left( n \right)}}{{dt}} = {A_t}{a_{d1t}}S\left( {{y_1}\left( n \right) - {y_3}\left( n \right)} \right) - 2{a_{d1t}}{y_8}\left( n \right) - a_{d1t}^2{y_7}\left( n \right) \hfill \\ \frac{{{y_9}\left( {n + 1} \right) - {y_9}\left( n \right)}}{{dt}} = {y_{10}}\left( n \right) \hfill \\ \frac{{{y_{10}}\left( {n + 1} \right) - {y_{10}}\left( n \right)}}{{dt}} = {A_t}{a_{d2t}}S\left( {{y_1}\left( n \right) - {y_3}\left( n \right)} \right) - 2{a_{d2t}}{y_{10}}\left( n \right) - a_{d2t}^2{y_9}\left( n \right) \hfill \\ \frac{{{y_{11}}\left( {n + 1} \right) - {y_{11}}\left( n \right)}}{{dt}} = {y_{12}}\left( n \right) \hfill \\ \frac{{{y_{12}}\left( {n + 1} \right) - {y_{12}}\left( n \right)}}{{dt}} = {A_t}{a_{d3t}}S\left( {{y_1}\left( n \right) - {y_3}\left( n \right)} \right) - 2{a_{d3t}}{y_{12}}\left( n \right) - a_{d3t}^2{y_{11}}\left( n \right) \hfill \\ \end{gathered} \right. $ |
在DSP Builder環境下，搭建的大腦皮層神經群模型及丘腦神經群模型如圖 4所示。
圖4
在DSP Builder中構建的大腦皮層神經群模型及丘腦神經群模型
(a)DSP Builder環境下的大腦皮層神經群模型; (b)DSP Builder環境下的丘腦神經群模型
Figure4. Neural mass model of the cortex and thalamus in DSP Builder(a) the NMM of a cortical voxel in DSP Builder; (b) the NMM of the thalamus in DSP Builder
仿真實驗中，采用歐拉法求解上述差分方程組，dt取0.001。根據表 1中的典型參數仿真出的δ、θ、α、β和γ波形如圖 5(b)所示。
圖5
平均突觸時間常數[a，b]對應的腦電各個節律信號
(a)平均突觸時間常數[
(a) areal distribution of the average synaptic time constants for
4 不同振蕩信號的參數選擇
表 1中列出了產生腦電各個節律振蕩信號的典型參數值，通過對比可以發現，在其他參數不變的情況下，增大C1的值，腦電的振蕩活動從beta向gamma頻段移動。同樣，保持C1不變，增大興奮性或抑制性的時間常數，振蕩活動從delta向theta頻段移動。從中可以發現，通過調整神經群的動態性能參數，可以產生從delta到gamma頻段的振蕩活動信號。然而，神經群模型中有很多參數，單獨研究每個參數的改變對模型輸出的影響是困難的。為此，我們通過改變興奮性時間常數a和抑制性時間常數b獲得產生相同腦電振蕩活動的參數區域。為了保證相應的振蕩活動，A/a、B/b的比值應為常數。對于每組a和b對應的振蕩活動信號，我們首先應該保證該信號為類正弦信號，進而求出該信號頻譜的峰值，將其分配到相應的腦電節律頻段區域。本文采用表 1中產生beta振蕩活動的典型參數值，參數a、b分別從5 s-1到200 s-1逐步增加，步進為5 s-1，其他參數保持不變，仿真結果如圖 5(a)所示。由圖 5可知，通過改變相應的參數，丘腦皮層神經群模型輸出的信號頻帶范圍會有相應的變化。圖 5(a)統計了腦電各個節律對應參數的變化范圍，如果要獲得多個節律的腦電波形，可以取相鄰腦電節律區域間的參數值，這樣輸出信號的頻帶會相對較寬，即包含多個節律的腦電信號。另一方面，經典神經群模型屬于單動力學神經群模型，單動力學神經群模型能夠得到不同頻率的窄帶信號，將單動力學神經群模型通過一定的方式連接起來，可以構建頻帶更寬的多動力學神經群模型，由此也可以得到含有多個節律的腦電信號。
5 結論
本文研究了丘腦皮層神經群模型，并采用FPGA實現了離散化快速計算和仿真，通過調整參數變化，給出了丘腦皮層神經群模型不同參數下對應的腦電節律信息。本文所建模型可以用于模擬不同腦功能狀態下的神經振蕩，生成正常與癲癇下的腦電數據以及各個節律下的腦電數據，并能對幅值進行調節。設計采用了MATLAB的Simulink與FPGA的DSP Builder相結合的方式，實現了軟件與硬件系統的有效融合，為神經群模型的計算機模擬提供了新的方案和思路。
