神经核团异常放电的数值模拟与FPGA实现

王超，逯迈

兰州交通大学光电技术与智能控制教育部重点实验室，甘肃兰州730070

前言

功能相似的神经细胞聚集成界限明显的群落，称为神经核团或神经节。基底神经节，简称基底节，是大脑皮层下一系列神经核团的总称，是神经系统中调节机体运动的重要结构，与随意运动的产生、肌肉紧张和感觉信息的处理等直接相关，同时还被认为与脑的感知、认知、记忆、学习等机制有着密切的联系［1］。基底节主要包括纹状体、苍白球外侧（the External Segment of the Globus Pallidus, GPe）、苍白球内侧（the Internal Segment of the Globus Pallidus,GPi）、黑质（the Substantia Nigra, SN）、丘脑底核（Subthalamic Nucleus, STN）等核团。基底神经节的病变会导致帕金森病、亨廷顿舞蹈病、癫痫等运动紊乱性疾病，是脑神经科学研究的重点［2］。

建立基底节的数学模型并利用计算机软件进行数值模拟有助于模拟基底节中各神经核团的放电特性及其运动控制功能等［3］。1989年，Albin 等［4］提出宏观的基底节直接-间接通路模型，研究帕金森病的发病机制；2001年，Gurney等［5］提出基底节的选择-控制模型，模拟基底节的选择-控制功能；2004年，Rubin等［6］提出一个精细的基底节神经核团网络结构，并研究深部脑刺激治疗帕金森病的作用机理；2001年，Humphries 等［7］利用Integrate-Fire（IF）神经元模型搭建基底节核团网络，模拟基底节的动作选择功能和振荡放电；2013年，Thibeault等［8］利用Izhikevich神经元模型搭建基底节核团网络并模拟帕金森病不同放电状态。

近年来利用硬件模拟神经元及神经网络的放电特性也得到越来越多的关注。2004年，Grass 等［9］利用Xilinx 公司的FPGA 芯片实现了Hodgkin-Huxley（HH）神经元模型；2008年，Cassidy 等［10］提出利用Izhikevich神经元模型硬件实现大规模神经网络的方法；2013年，Paolo 等［11］将Izhikevich 模型搭建的神经网络应用于脑机接口；2015年，王金龙等［12］基于FPGA 芯片实现HH 模型，为硬件模拟神经元提供了新的参考。但是大部分研究集中于神经元及神经元网络的硬件实现，对于神经核团网络的硬件实现只有很少一部分。2006年，Prescott 等［13］设计制造基于基底节选择-控制模型的机器人，模仿基底节的选择与控制功能；2015年，Yang等［14］利用简化的突触模型在FPGA芯片实现一种基底节核团网络模型。

模拟神经核团的不同放电状态可以深入研究相关神经性疾病的发病机制，探究治疗方式的作用机理，如深部脑刺激、经颅磁刺激等。采用价格低廉、运行速度快的FPGA 芯片实现的神经核团网络具有重要的应用价值，如替代动物活体实验，降低医学实验成本；应用于医疗康复工程，替代受损的神经核团［15］；脑机接口技术，辅助神经性疾病患者的生活，为研究和治疗神经性疾病提供新的思路［16］；嵌入智能机器人，模拟人类简单的选择运动功能，实现类脑式的机器人操控，为智能机器人的设计提供新的方向。

本研究依据Rubin 等［6］提出的基底节神经核团网络结构，以Izhikevich 神经元模型为神经核团网络的节点，选用化学突触连接各神经元，运用Simulink软件和DSP Builder软件分别对神经核团网络进行建模。对比分析两种软件仿真结果的相对标准误差，验证建模的一致性。通过DSP Builder 软件与Quartus II 软件联合编译的方法，对神经核团网络进行FPGA 硬件实现，再现基底节神经核团正常放电和异常放电两种状态。

1 基底节神经核团网络模型

2004年，Rubin 等［6］提出一种基底节神经核团网络模型，研究各神经核团的放电波形和深部脑刺激作用机理，其提出的网络模型包含GPe、STN、GPi、丘脑皮层（Thalamic Cortex,TC）这4 种神经核团。如图1所示，红色箭头表示化学突触兴奋性输入，黑色箭头表示化学突触抑制性输入，绿色箭头表示输入到TC 神经元的皮层运动信号，3 个STN 神经元兴奋性输入到1 个GPe 神经元，2 个GPe 神经元抑制性输入到1 个STN 神经元，1 个GPe 神经元抑制性输入到1个GPi 和1 个STN 神经元兴奋性输入到同一个GPi神经元，最后由GPi 神经元汇总输入到TC 神经元［6］。在异常放电时，来自GPi节律性的抑制性突触输入减弱了TC 中继神经元对来自皮层的运动信号的响应能力，导致帕金森症的运动障碍。

图1 基底节神经核团网络Fig.1 Network of basal ganglia nuclei

2 Simulink建模仿真

2.1 神经元模型

搭建基底节神经核团网络首先要选择一个合适的神经元模型。HH 模型过于复杂，消耗硬件资源多；IF 模型放电模式单一，不能表现神经核团的多种放电模式［17-18］。本研究选用的Izhikevich 神经元模型可以表现20 种放电模式，其中包括神经核团的簇放电和峰放电［19］。Izhikevich神经元数学模型如下：

其中，V是神经元膜电位，单位为mV；u是膜电位恢复变量，表示K+通道的激活和Na+通道的失活，膜电位恢复变量u为膜电位V提供负反馈；a、b、c、d为无量纲参数；Iapp为外部输入电流；Isyn为输入该神经元所有突触电流的总和。

2.2 化学突触模型

化学突触采用文献［20］中的模型，该模型可以实现神经核团网络中化学突触的抑制性输入和兴奋性输入，其表达式为：

其中，Isyn表示突触后电流；Gsyn表示耦合强度，单位为mS/cm2；S表示突触后膜上化学门控离子通道开放比例；ES是突触可逆电位，单位为mV；F(Vpre)表示突触前神经元膜电位的S型函数；α和β分别表示离子通道打开和关闭的速率；Vpre是突触前神经元的膜电位，单位为mV；θsyn是突触前可逆电位，单位为mV［21］。当ES=0 mV时，表示化学突触兴奋性输入；当ES=-75 mV，表示化学突触抑制性输入。表1为各神经元间突触系数Gsyn的值［14］。

表1 突触系数Gsyn（mS/cm2）Tab.1 Synaptic coefficient Gsyn(mS/cm2)

2.3 皮层运动信号

根据Rubin 等［6］和Yang 等［14］的研究，输入到TC核团的皮层运动信号ISM用式（3）模拟：

其中，H为阶跃函数；皮层电流最大值imax=30 pA；脉冲间隔时间参数分别为ρ=25 ms，δ=3 ms。波形如图2所示。

图2 皮层运动信号ISMFig.2 Cortical motor signal ISM

根据Thibeault 等［8］和Yang 等［14］的研究，对a、b、c、d参数进行不同设置，可以表现出GPe、STN、GPi、TC神经元不同放电特性。具体参数设置见表2和3。

表2 神经元正常放电参数Tab.2 Parameters for neurons in normal discharge

表3 神经元异常放电参数Tab.3 Parameters for neurons in abnormal discharge

2.4 Simulink仿真结果

图3为Simulink 软件仿真结果，分别为GPe、STN、GPi 和TC 神经元正常放电和异常放电的数值模拟波形。可以观察到GPe 和GPi 神经元异常放电时，放电波形从峰放电变为簇放电。当STN 神经元异常放电时，放电频率相对于正常放电更快。当TC神经元异常放电时，不能对运动信号做出正常响应，出现信号的缺失。

图3 GPe、STN、GPi和TC神经元正常放电和异常放电的Simulink模拟结果Fig.3 Simulink simulation results of GPe,STN,GPi and TC neurons in normal and abnormal discharges

3 DSP Builder建模仿真

DSP Builder 软件是Altera 公司提供的一种系统级开发工具，将建模仿真和硬件实现结合起来，在Matlab/Simulink 环境下以库文件的形式进行图形化的设计和仿真，并将Quartus II 软件作为底层设计工具置于后台。将建好的DSP Builder 设计文件（.mdl）通过Signal Compiler 转化成硬件描述语言HDL 文件（.vhd），进而可以在Quartus II 环境下进行编译和FPGA硬件配置。

3.1 神经元及突触建模

在硬件实现时，Izhikevich 原模型资源消耗相对过大，不利于神经核团网络的实现。2012年，Soleimani 等［22］利用四阶线性近似的方法简化了Izhikevich模型，模型如下：

根据Soleimani 等［22］提出的选择算法，设定m1=0.75，m2=0.375，m3=11，m4=62.5。

对神经元和化学突触进行建模首先要选择合适的离散方法，考虑硬件资源的消耗，采用欧拉法对神经元模型和化学突触进行离散化处理。用欧拉法对Izhikevich模型进行处理，得到差分方程：

图4为神经元的建模过程，采用移位器和加法器替代乘法器，减少乘法器的使用，节约硬件资源［23］。移位器可以根据命令向左或向右移动数据位，移位操作为算术移位。用移位器和加法器的运算值趋近于模型中m1、m2、a、b和△T的值，以此替代乘法运算，如单位时间△T=0.02≈1/64+1/256。

图4 Izhikevich模型建模流程Fig.4 Establishment of Izhikevich model

运用DSP Builder软件对神经元模型和突触模型分别建模，然后用化学突触将各个神经元按照神经核团网络结构依次连接，构成基底节神经核团网络。在DSP Builder 建模时需要考虑时序，可以更准确地将突触的延时表现出来，取神经元模型计算1次的时间作为突触传导的延时。

3.2 DSP Builder软件仿真结果

考虑硬件资源的限制，本研究实现的神经核团网络由3 个GPe 神经元、3 个STN 神经元、3 个GPi 神经元、1 个TC 神经元和30 个化学突触构成，最后由GPi 神经元汇总输入到1 个TC 神经元。图5为DSP Builder软件的仿真结果，分别为GPe、STN、GPi和TC神经元正常放电和异常放电的仿真结果，与Simulink仿真结果基本一致。

图5 GPe、STN、GPi和TC神经元正常放电和异常放电的DSP Builder模拟结果Fig.5 DSP Builder simulation results of GPe,STN,GPi and TC neurons in normal and abnormal discharges

3.3 波形误差分析

为了量化评估Izhikevich线性简化模型搭建的神经核团网络的精确性，采用标准误差公式进行计算分析［22］。选取一个周期内两种软件的模拟结果，即0～150 ms 样本值进行对比分析。因为Simulink 软件没有时序控制模块，为了更好地对比两款软件的仿真结果，DSP Builder 软件仿真中不考虑突触延时，并对仿真时间进行调整。标准误差（Root Mean Square Error，RMSE）计算见式（6）：

其中，Vorigin为原模型电位，VPWL为线性简化模型电位。计算结果如表4所示，可以发现各神经元在不同状态的标准误差均小于0.1，其中TC 神经元的误差最大。造成误差的主要原因是线性简化模型中利用移位器代替了乘法器，a、b、c、d等值有略微改变，此外硬件实现中考虑资源的消耗，对输出波形的位数进行了限制。

表4 各神经核团的仿真结果标准误差Tab.4 Root mean squared error of simulation results of nuclei each neuron

4 神经核团的硬件实现

本研究所用FPGA 开发板为Altera 公司CycloneIV 系列EP4CE115F23I7，数模转换模块为DAC900E，采用架构在MATLAB/Simulink 上的DSP Builder 辅助设计工具箱对神经元、化学突触和神经核团网络进行建模，VerilogHDL 语言描述分频文件。所用各软件版本为Quartus II 13.0、DSP Builder 13.0以及MATLAB 2013。

将DSP Builder软件中建模完成的神经核团网络进行编译，然后导入Quartus II 软件加入分频文件并配置管脚，运行成功后导入FPGA 芯片，通过DA 芯片输出到示波器，观察到各神经元放电波形。图6为GPe、STN、GPi 和TC 神经元正常放电和异常放电时的硬件实现波形。图中的毛刺是外部干扰的原因，并不影响结果的正确性。图7为实现的神经核团网络实物图。

图6 GPe、STN、GPi和TC神经元正常放电和异常放电的硬件结果Fig.6 Hardware implementation results of GPe,STN,GPi and TC neurons in normal and abnormal discharges

图7 神经核团网络硬件模型Fig.7 Hardware model of the network of nuclei

5 讨论与分析

帕金森病是中老年人常见的一种神经退行性疾病。目前帕金森病的致病原因尚不清楚，但是研究表明与环境、遗传和衰老等多种因素共同作用有关。典型的临床症状有肌肉僵直，静止性震颤和运动迟缓等。其中，TC 神经元中继能力的丧失是帕金森病的主要症状之一［24］。Rubin 等［6］和Yang 等［14］在研究中提出一种丘脑神经元对皮层运动信号的中继可靠性指标（Reliability Index, RI），是评价帕金森状态的一个重要指标。

其中，Nerrors表示丘脑元没有准确响应皮层运动信号脉冲的个数，NSM表示皮层运动信号脉冲的总个数。当RI＞0.9 时，表示正常状态；RI＜0.6 时，表示帕金森状态；其余为中间状态。

图8、图9为实现的神经核团网络在异常放电和正常放电时TC 神经元对皮层运动信号的响应波形。上栏波形为皮层运动信号，下栏波形为TC 对运动信号的响应。

图8 丘脑皮层神经元正常放电Fig.8 TC neuron in normal discharge

选取其中20 个运动信号脉冲个数进行计算（如图9红框）。正常放电时，RI=1.0＞0.9，TC 神经元没有出现信号缺失；异常放电时，RI=0.3＜0.6，TC 神经元出现信号的缺失，证明实现的基底节神经核团表现为帕金森状态。

图9 丘脑皮层神经元异常放电Fig.9 TC neuron in abnormal discharge

RI 可以验证硬件实现的基底节神经核团网络可以展现帕金森病的一种放电状态。但是由于RI过于单一，仅限于理论上的验证，还要结合各神经核团的放电波形等其他因素进一步验证。

帕金森病与苍白球核团和STN 核团的异常放电有着明显联系。在帕金森病相关临床研究和动物模型生理实验中发现，GPe 和GPi神经元出现明显的爆发式簇放电［25］，并且GPi低频爆发式簇放电与帕金森症的震颤有关［26］。从模拟的神经核团放电波形可以观察到，在异常状态时，GPe 和GPi 神经元再现了帕金森病理模型的簇放电。同时可以观察到STN 神经元的放电波形出现异常，放电频率加快。这也符合病理模型中对STN 核团的研究，STN 神经元放电活动的改变引发明显的行为异常。基底节SN 致密部多巴胺能分泌的减少影响GPe 核团对STN 核团的抑制，使STN 神经元连续非正常放电，并且放电频率增加，导致输出核团GPi 的兴奋，过度抑制皮层的相关运动区域，进而导致人体产生运动迟缓，动作保持困难和震颤等运动功能障碍［27］。

现在大多数研究集中于实现单个神经元和链式神经网络等简单的神经网络或者稍微复杂的小世界神经网络［28］，但是对于功能性的神经核团网络研究较少，并且更多的研究限于数值模拟仿真。如Thibeault 等［8］通过数值模拟研究基底节神经核团放电特性并研究脑深部电刺激对神经核团网络的影响，但是没有涉及硬件实现。Humphries 等［7］、Prescott等［13］选用IF 模型或者人工神经元等相对简单的神经元模型硬件实现神经核团网络，但是由于神经元模型简单，不能表现神经核团多种放电模式。本研究则选用动力特性较为复杂的Izhikevich模型搭建神经核团网络，不仅通过软件模拟神经核团的两种放电状态，而且结合FPGA进行了硬件实现。

6 结论

多种神经性疾病的产生与神经核团的异常放电有着重要联系，数值模拟和FPGA 硬件模拟神经核团异常放电为研究神经核团放电状态与神经性疾病间的联系提供了更多的研究方式。本研究在Simulink、DSP Builder 和Quartus II 联合仿真编译的基础上，通过单片大资源FPGA 模拟了基底节神经核团网络异常放电和正常放电两种放电模式。采用相对标准误差公式量化分析了线性简化Izhikevich 模型的误差，验证了两种软件建模的一致性。在神经核团网络正常状态时，RI=1.0，表明TC神经元如实地中转了传入的每一个运动信号；当核团网络处于异常状态时，RI=0.3，小于文献［14］所实现的神经核团帕金森状态的RI 值（0.55），表明TC 神经元对运动信号的中转出现明显的缺失。此外，结合相关帕金森病理实验记录进行分析，证明实现的神经核团网络可以模拟帕金森病的一种放电状态。

由于FPGA 芯片资源的限制，本研究只实现了一个小规模的神经核团网络，并且实现的神经核团放电模式较少，仅有峰放电和簇放电。后续的研究工作可以在这个基础上完成大规模、多种放电功能的神经核团网络的硬件实现。用移位器和加法器代替乘法器，可以大大节约FPGA 资源的消耗，为以后大规模神经元网络的实现提供新的参考。此外，硬件实现的基底神经核团可以应用于替代动物活体实验、脑机接口和医疗康复工程的研究中。