帕金森病的脑深部电刺激治疗建模研究

Modeling Research on the Treatment of Deep Brain Electrical Stimulation in Parkinson’s Disease

摘要: 本文首先通过神经元Hodgkin-Huxley模型的基础上，融入化学突触模型，通过兴奋以及抑制型突触模型将多个单体神经元连接汇成神经核团，研究大脑皮层不受刺激时的神经核团内部电位。其次将多个核团(Cortex, Striatum/dMSN/iMSN, SNc, GPe, GPi/SNc, STN, Thalamus)按照基底神经节内部神经核团连接方式搭建，其间的兴奋或抑制方向一定。将去掉黑质SNc神经核团的基底神经节回路模型作为帕金森病态模型，仿真健康状态与帕金森病态的基底神经节回路模型，并阐述两种状态所特有的特征指标。接着基于PID自动控制理论，对帕金森病态的基底神经节靶点施加正弦高频电刺激，分析评价DBS在不同靶点位置的刺激效果的优劣，同时获取优化后的电刺激参数，为进行DBS的治疗提供参考。最后基于基底神经节回路模型，旨在寻找除已知的刺激靶点外，是否还存在有其他最优的靶点。通过本文的研究发现，STN靶点在治疗帕金森病态过程中以6 mV的误差结果，确定为最佳刺激靶点，而其他位置所选定的靶点iMSN在治疗过程中存在有17 mV的误差，为后期选取电刺激靶点提供了数据和模型依据，深度了解帕金森病态与正常状态的电位发放相关参数信号指标。

Abstract: In this paper, based on the Hodgkin-Huxley model of neurons, the chemical synaptic model is inte-grated into, and multiple monomer neurons are connected into neural nucleus clusters through ex-citatory and inhibitory synaptic models to study the internal potential of neural nucleus clusters when the cerebral cortex is not stimulated. Secondly, multiple nuclear clusters (Cortex, Stria-tum/dMSN/iMSN, SNc, GPe, GPi/SNc, STN, Thalamus) are constructed according to the connection method of the neural nuclear clusters inside the basal ganglia, and the direction of excitement or inhibition in the meantime is certain. The basal ganglion circuit model with the removal of the nigra SNc nucleus cluster is used as the model of the state of Parkinson’s disease, the basal ganglion cir-cuit model of the healthy state and the state of Parkinson’s disease is simulated, and the character-istic indicators unique to the two states are expounded. Then based on the PID automatic control theory, sinusoidal high-frequency electrical stimulation is applied to the basal ganglion targets of Parkinson’s disease, and the advantages and disadvantages of the stimulation effect of DBS at dif-ferent target locations are analyzed and evaluated. At the same time, the optimized electrical stim-ulation parameters are obtained to provide a reference for the treatment of DBS. Finally, based on the basal ganglion circuit model, the purpose is to find out whether there are other optimal targets in addition to the known stimulation targets. Through the research in this paper, it is found that the STN target is determined as the best stimulation target with an error result of 6 mV during the treatment of Parkinson’s disease, while the iMSN of the selected target in other locations has an er-ror of 17 mV during the treatment process, which provides a data and model basis for the selection of electrical stimulation targets in the later stage, and provides an in-depth understanding of the potential and related parameter signal indicators of the normal state of Parkinson’s disease.

1. 引言

帕金森病(Parkinson’s disease, PD)是一种严重的神经退行性疾病 [1] ，其临床表现包括运动发起困难、肌肉僵直和静止性震颤等运动症状，而且常常伴发睡眠紊乱及情绪情感障碍等非运动症状 [2] [3] [4] 。脑深部电刺激(deep brain stimulation, DBS)是一种治疗运动障碍性疾病和神经精神疾病的外科治疗手段，对于帕金森疾病有着良好的治疗效果，具有微创、可逆、可调节等优点，因此目前已成为中晚期帕金森疾病的最有效的治疗方案之一。DBS是通过立体定向方法进行精确定位，输入高频电刺激，从而改变相应核团的兴奋性，以达到改善帕金森病症状。

由于脑电信号微弱而复杂，目前传统的开环DBS仍广泛应用于当前临床治疗中，而闭环DBS还处于初步研究阶段。并随着应用算法等技术的进步以及对神经生理学的深入了解，闭环DBS中的有一些障碍正逐渐被克服。然而，直接以PD患者作为对象，直接测试控制算法或参数的调节是不现实的，计算模型解决了动物和临床实验结果受多种因素影响难以再现的问题，实现了在大量数据下表征帕金森状态的生物特征识别。帕金森的计算模型具有再现帕金森状态的能力有助于探索DBS的治疗机制，为闭环DBS 的研究提供平台。虽然计算模型中的DBS 参数可能无法准确地转换到患者身上，但计算模型可以快速、有效地测试和识别潜在的算法，然后经过动物实验后才能转入临床实践。基于计算模型的闭环DBS 研究可以使模型作为一个虚拟病人进行刺激参数调节，刺激模式，刺激位点，闭环控制算法等测试从而更系统、更有效地确定最佳参数设置，达到最佳自适应刺激的效果。

E. Rubin和D. Terman等人基于电导的Hodgkin-Huxley (HH)类神经元方程 [5] ，从单个神经元的生理特征和突触连接建立了包括STN、GPe、GPi和TH细胞在内的BG网络计算模型。E. Rubin和D. Terman等人构建的模型具有生理意义和完整的数学分析，是PD计算模型研究中的一项突破性进展。

本文的目的研究脑深部电刺激治疗的最优靶点选择和DBS参数(刺激强度和刺激频率等)优化以保证DBS的最佳治疗效果。本文以单个神经元模型为基础，逐渐完善神经系统之间的联系，将多个神经元组合成团，研究健康状态和PD状态下Thalamus核团的信号，并进行对比，然后模拟DBS治疗PD寻找到电刺激参数以及电刺激靶点的最优解。

2. 模型假设与思路

2.1. 模型假设

在建模的过程中，在不影响模型意义与计算精度的前体下，为了使模型简单明确，建立了如下假设：

1) 神经元单体一致。

2) 神经元之间的连接通过化学突触模型中的兴奋性突触模型(AMPA)以及抑制型突触模型(GABA)来模拟。

3) 神经核团之间根据图1中各个核团之间的兴奋、抑制关系来采用相应的化学突触模型来连接。

4) 大脑皮层收到的外接刺激等同于高频正弦电流刺激。

5) 在没有任何刺激的情况下，神经元膜电位为−65 mV。

6) 本文所有模型不考虑温度影响。

2.2. 参数表(见表1)

Table 1. Parameter table

表1. 参数表

2.3. 模型搭建

在静息状态下，神经元膜内外的离子浓度不同形成神经细胞的膜电位。当神经系统受到外界刺激时，膜电位产生的动作电位可以形成电位发放.这些动作电位的峰发放和簇发放形成神经系统的信息传递编码，典型的神经元膜电位可由Hodgkin-Huxley模型 [6] 描述，模型的等效电路图如下图1：

Figure 1. HH model equivalent circuit

图1. HH模型等效电路

HH模型是一个准确描述神经元的动作电位产生和传播的定量模型。它包含离子通道激活和失活、动作电位等神经元特性，在离子层面对神经元的电活动进行阐释。HH模型可以用如下的四阶微分方程表述 [7] ：

$\begin{array}{l}C\frac{\text{d}V}{\text{d}t}=-g_{N_a}{m}^{3}h\left(V-V_{N_a}\right)-g_K{n}^4\left(V-V_K\right)-g_L\left(V-V_L\right)+I_{external}+I_{synapse}\\ \frac{\text{d}m}{\text{d}t}=-\left(a_m+b_m\right)m+a_m,\\ \frac{\text{d}h}{\text{d}t}=-\left(a_h+b_h\right)h+a_h,\\ \frac{\text{d}n}{\text{d}t}=-\left(a_n+b_n\right)n+a_n\end{array}$ (1)

如图1所示，将细胞膜的电荷储存能力模拟为电容，其电容用C表示，C = 1 μF，将不同的离子通道模拟为电阻，将电动势模拟为膜内外离子浓度差异产生的电化学电位。在Hodgkin-Huxley模型中，

1) V(t)表示神经元的膜电位。

2) m(t)，h(t)，n(t)描述细胞膜内外离子通道的电导特性。

3) g_Na = 120 mS/cm²，g_k = 36 mS/cm²，g_L = 0.3 Ms/cm²分别对应钠离子、钾离子和泄漏电流关于细胞膜的电导系数的最大值。

4) V_Na = 50 mV，V_k = −77 mV，V_L = −54.5 mV分别对应钠离子、钾离子和泄漏电流的反向电压。

5) 其中a函数和b函数是和膜电位相关的门变量率函数 [8] ，其表达式如下所示：

$a_m=\frac{0.1\left(V+40\right)}{1-{\text{e}}^{-0.1(V+40)}},a_h=0.07{\text{e}}^{-0.05(V+65)},a_n=\frac{0.01\left(V+55\right)}{1-{\text{e}}^{-0.1(V+55)}}$ (2)

$b_m=4{\text{e}}^{-(V+65)/18},b_h=\frac{1}{1+{\text{e}}^{-0.1(V+35)}},b_n=0.125{\text{e}}^{-(V+65)/80}$ (3)

6) I_external对应外界对神经元的刺激影响。

7) I_synapse对应神经元之间的化学突触电流，一般包含兴奋突触和抑制突触两类。

8) 神经元的电位发放就是神经元的动作电位–时间关系图。

2.4. 数值模拟外界刺激

1) 直流电刺激：

根据Hodgkin-Huxley模型参数，使用matlab/simulink软件将HH神经元模型搭建完成。根据参考文献 [8] ，设置直流电刺激，得出动作电位–时间图如下图2(a)，放大后得到的单个频率图如下图2(b)。我们将该结果与参考文献 [9] 中的模型结果进行对比，可以发现动作电位–时间图基本一致，且单个频率图的趋势也类似，图2(c)、图2(d)为文献中模型的结果图。通过该对比结果，表明我们根据 [8] 中的资料所搭建出的HH神经元模型准确，可以用作下一步的基底神经节回路模型的搭建。

Figure 2. Action potential under DC stimulation-time chart. (a) The action potential-time diagram of the model in this article; (b) Enlarged view of a single frequency in this article; (c) Reference action potential-time chart; (d) Reference single frequency enlargement diagram

图2. 直流刺激下动作电位–时间图。(a) 本文模型动作电位–时间图；(b) 本文单个频率放大图；(c) 参考文献动作电位–时间图；(d) 参考文献单个频率放大图

依据电位发放分类方式进行判断，特征指标为：在该直流刺激下的电位发放为峰发放，通过单个频率示意图中的数据可以得出该峰发放的振幅为75.22971，频率为50 HZ。

2) 交流电刺激：

设置电刺激为交流刺激( $I=9.6\sin \left(0.05\pi \cdot t\right)$ )对神经元进行刺激，得到的动作电位–时间图如下图3(a)所示，放大后的单个频率如下图3(b)所示。

图3(b)的特征指标为：振幅为85.537，静息间隔为28 ms，激活时间为12 ms，峰峰间距为12 ms，簇发放周期为50 ms。

Figure 3. (a) Action potential-time chart under AC stimulation; (b) Enlarged view of a single frequency

图3. (a) 交流电刺激下动作电位–时间图；(b) 单个频率放大图

3. 基底神经节神经回路理论模型

3.1. 兴奋型突触与抑制型突触模型

在本次建模中，神经元之间的突触连接分为兴奋突触以及抑制突触两种，其中，兴奋型突触模型的代表为AMPA，抑制型突触模型的代表为GABA。其中，无论是兴奋型突触或抑制型突触，基本原理都归于突触前与突触后的电压差作为初始条件，最终通过兴奋或抑制模型，得出电流对下一个神经元进行刺激作用 [10] 。

通过 [8] 中的兴奋型突触模型以及抑制型突触模型的相关参数资料，搭建该部分突触模型，结合第一问中所搭建的单个神经元模型组成基底神经节回路模型。突触电流 $I=g_{AMPA}\cdot r\cdot \left(V_{post}-E_{AMPA}\right)$ ，其中 $g_{AMPA}$ 为最大突触电导， $E_{AMPA}$ 为逆转电位。门控变量r代表突触离子通道开放的比例。我们假设r服从一级动力学。具体的突触模型以及具体参数如下面的函数表达式所示 [11] 。

1) 兴奋型突触模型——谷氨酸能突触(AMPA)

$\{\begin{array}{c}{I}_{synapse}^{AMPA}=g_{AMPA}\ast r\ast \left(V_{post}-E_{AMPA}\right)\\ \frac{\text{d}r}{\text{d}t}=a\ast S\left(V_{pre}\right)\ast \left(1-r\right)-b\ast r\\ S\left(V_{pre}\right)=\frac{1}{1+\exp \left(-\frac{V_{pre}-V_p}{k_p}\right)}.\end{array}$ (4)

式中： $g_{AMPA}=0.35~1nS$ 是最大电导，r受体开放状态的比例， $V_{pre}$ 和 $V_{post}$ 是突触前和突触后电压， ${\mathrm{E}}_{AMPA}=0mV$ 是逆转电位， $\alpha =1.1\times {10}^6{\text{M}}^{-1}{\text{s}}^{-1}$ ， $\beta =190{\text{s}}^{-1}$ ， $V_P=2\text{mV}$ ， $K_P=5\text{mV}$ 。

2) 抑制型突触模型——GABA能突触(GABA)

$\{\begin{array}{c}{I}_{synapse}^{GABA}=g_{GABA}*r*\left(V_{post}-E_{GABA}\right),\\ \frac{\text{d}r}{\text{d}t}=a*S\left(V_{pre}\right)*\left(1-r\right)-b*r,\\ S\left(V_{pre}\right)=\frac{1}{1+\exp \left(-\frac{V_{pre}-V_p}{k_p}\right)}.\end{array}$ (5)

式中： $g_{AMPA}=0.25~1.2nS$ 是最大电导，r受体开放状态的比例， $V_{pre}$ 和 $V_{post}$ 是突触前和突触后电压， ${\mathrm{E}}_{AMPA}=-80mV$ 是逆转电位， $\alpha =5\times {10}^6{\text{M}}^{-1}{\text{s}}^{-1}$ ， $\beta =180{\text{s}}^{-1}$ ， $V_P=2\text{mV}$ ， $K_P=5\text{mV}$ 。

3.2. 神经元连接方式

神经元之间主要依靠化学突触和电突触进行突触连接，本文主要考虑化学突触。在神经元化学突触连接方式中，由于神经元之间的联系较为复杂，我们主要考虑了以下几种连接方式：串联型、分层结构、全连接型，如图4所示。

串联型：上一个神经元的轴突与下一个神经元的树突相连接，为一一对应关系。该类型较为简单，神经元单体之间的相互联系较少，便于计算，但与实际神经元电位发放情况的一致性较差。如图4(a)所示。

分层结构：将多个神经元按照层级进行区分，如图4(b)为1-4-1的分层结构，第一层的神经元分别向第二层的4个神经元传递电流信号，第二层的4个神经元紧接着共同连接第三层的1个神经元。如下图4(b)所示。

全连接型：所有神经元都与其他任一神经元由相互连接关系，相互传递刺激信号，整个神经核团模型更加完整，更符合实际神经元连接方式，但整体计算量增加偏大，影响整个模型的运行速度。如下图4(c)所示。

Figure 4. Ganglion connection diagram: (a) Series type; (b) layered type; (c) Fully connected type

图4. 神经节连接方式图：(a) 串联型；(b) 分层型；(c) 全连接型

通过以上三种连接方式对比，再结合模型复杂程度，我们最终决定选用分层结构作为神经元之间的连接方式，由于每一个神经核团内部有5~10个神经元，故我们使用6个神经元作为一个神经核团，按照1-4-1的方式进行连接，其中，每一个神经元都以兴奋性突触模型及抑制性突触模型连接，通过电位差计算得出的微小电流相加之和，作为下一个神经元的刺激信号。以此来形成基底神经节神经回路的理论模型，作为后续研究健康状态与帕金森病态之间的联系与差别。

3.3. 回路中的神经元电位发放

将基底神经节神经回路模型搭建完成后，将GPe、GPi/SNr、STN以及Thalamus中添加示波器以此观测该模块中的神经元电位发放，为增加观测的完整度，我们同时取神经核团内部的第1层中的一个神经元，第2层中的两个神经元以及第3层中的一个神经元作为观测目标，分别为x-1；x-2.1；x-2.2；x-3 (x为GPe、GPi/SNr、STN、Thalamus)，根据相关文献 [12] 中，在没有GPi或皮层输入的情况下，丘脑–皮层细胞被认为在−60 mV附近处于静止状态。

最终模型跑出来的结果中可以发现，如下图5所示，神经元电位在前期会有一个较大的波动，之后便都处于−65 mV附近处于静止状态，与文献资料 [11] 中的数据基本一致，可以认为该模型搭建结果正确。

从下图5中的结果中可以发现，GPe、GPi/SNr、STN以及Thalamus中所选取的四个神经节单元所表现出的电位发放基本一致，前期的电位波动会有细微区别，100 ms之后基本处于静止状态，保持在−65 mV附近。并且，可以发现，同位置的神经节单元的电位发放也较为相似，通过该模型可以发现，由于每一个核团内部都为6个神经元单元且连接方式一致，故其核团内部相同位置的神经元的电位特征都极为相似。

Figure 5. Neuronal potential distribution diagram

图5. 神经元电位发放图

4. 健康状态以及帕金森病态之间的神经节回路电位发放对比

4.1. 模型设置

当大脑皮层受到刺激后，整个基底神经节神经回路处于兴奋状态。同样的，在GPe、GPi/SNr、STN以及Thalamus中添加示波器以此观测该模块中的神经元电位发放，这次仅针对神经核团的电位发放而不需要观测内部的神经元电位发放。

4.2. PD状态下基底神经节

PD状态下，基底核网络计算模型中各个核团神经元的放电模式都是一种无规则的状态，并且TH神经元能够准确地中继大脑皮层感觉运动区的输入信号。图6表示在1000 ms的仿真时间内基底核网络计算模型在帕金森病态下的放电情况。为了便于分析，每种类型的神经元细胞都进行多次频率放大研究，如下图6右侧图所示。帕金森病态下，TH神经元细胞对于大脑皮层感觉运动区的每个输入脉冲都能有与之对应的关系，且TH神经元的幅值大概在70 mV左右。STN神经元细胞在帕金森病态情况下，放电的频率较低且稀疏无规则，相同的，GPe和GPi的放电频率也较低，同时GPe神经元的放电也是稀疏、无规则的放电模式。在没有黑质SNc释放多巴胺的情况下，STN和GPe神经元细胞组成的兴奋性–抑制网络进行自发性振荡。并且由于STN神经元细胞对GPi神经元细胞有兴奋作用、以及GPe神经元细胞对GPi神经元细胞有抑制作用，两者对GPi神经元细胞的放电情况有比较大的影响，也会影响GPi神经元细胞对外界输入的接收情况。

Figure 6. Distribution of potentials of each nuclear group in the PD state

图6. PD状态下的各核团电位发放情况

GPe为簇发放，其特征指标：振幅为88.888，静息间隔为40.03 ms，激活时间为30.54 ms，峰峰间距为10.43 ms，簇发放周期为85.49 ms。

GPi/SNr为簇发放，其特征指标：振幅为98.32，静息间隔为52.44 ms，激活时间为44.53 ms，峰峰间距为14.02 ms，簇发放周期为92.55 ms。

STN为簇发放，其特征指标：振幅为76.93，静息间隔为53.39 ms，激活时间为46.53 ms，峰峰间距为10.82 ms，簇发放周期为100.84 ms。

Thalamus为簇发放，其特征指标：振幅为76.32，静息间隔为51.34 ms，激活时间为50.28 ms，峰峰间距为11.32，簇发放周期为99.24 ms。

4.3. 健康状态下的基底神经节

而在健康状态下，当大脑皮层收到刺激后，各个核团神经元都能接受到相应的输入信号，图7表示了1000 ms的仿真时间内基底核网络计算模型在健康状态下的放电情况。

Figure 7. Distribution of potentials of each nuclear group in a healthy state

图7. 健康状态下的各核团电位发放情况

在各核团电位放电情况下，可以发现，相对于帕金森病态下的放电情况，健康状态下的下丘脑在稳定接收大脑皮层受刺激后发出电信号的同时，其表现出来的电位发放频率更高，更稳定，幅值基本稳定在−75 mV附近。与此同时，由于SNc稳定释放出多巴胺对整个回路进行平衡，可以发现，健康状态下的基底神经节各核团的电位发放相对于帕金森病态的电位发放更加稳定，且频率更高。整体回路中的直接通路、间接通路以及超直接通路在稳态电刺激的影响下相互平衡，共同维持基底神经节各核团的电位发放，为正常活动提供了强大的保障。

GPe为峰发放，其特征指标：振幅为74.35，频率为67.43 HZ。

GPi/SNr为峰发放，其特征指标：振幅为74.44，频率为31.33 HZ。

STN为峰发放，其特征指标：振幅为74.46，频率为51.71 HZ。

Thalamus为峰发放，其特征指标：振幅为74.54，频率为42.70 HZ。

5. 刺激靶点的选择与电刺激参数优化

现如今的帕金森治疗手段主要也都以STN以及SPi/SNr作为靶点 [13] 来进行电刺激治疗，由于STN与黑质间存在非常紧密的双向神经纤维联系，STN既可能调控黑质中的非多巴胺能神经元，也可能影响其中的多巴胺能神经元。而STN发出的神经纤维除投射至GPi外，还投射至GPe、黑质网状部、纹状体和脚桥核等核团。并且，STN还与基底神经节中的GPe和黑质网状部间存在直接的双向环路联系。这些STN的传入传出联系成为STN神经元活动，特别是STN运动调控和非运动调控功能的神经解剖学基础 [14] 。

5.1. PID算法

PID算法是控制行业最经典、最简单、而又最能体现反馈控制思想的算法。PID控制器问世至今已有近70年历史，它以其结构简单、稳定性好、工作可靠、调整方便而成为工业控制主要技术之一。当被控对象结构和参数不能完全掌握，或不到精确数学模型时，控制理论其它技术难以采用时，系统控制器结构和参数必须依靠经验和现场调试来确定，这时应用PID控制技术最为方便。即当我们不完全了解一个系统和被控对象﹐或不能有效测量手段来获系统参数时，最适合用PID控制技术。PID控制算法，实际中也有PI和PD控制。PID控制器就是系统误差，利用比例、积分、微分计算出控制量进行控制，如图8所示。在本题中，我们通过PID实现基底神经节回路的闭环控制，在丘脑可靠性的时变基础上建立刺激规则，通过PID控制器根据丘脑可靠性的估计水平按比例调整刺激，从而实现下丘脑从帕金森状态下恢复可靠性，实现对帕金森的治疗。

Figure 8. PID schematic

图8. PID原理图

如今DBS-STN [15] 也成为了使用最为频繁的治疗帕金森病症的方式之一，然而，根据模型回路搭建情况发现，若对GPi进行电刺激，也可以通过整个回路使得下丘脑部分收到相对更兴奋的信号，故我们在本次模拟治疗帕金森病症的过程中，选择分别刺激靶点STN以及靶点GPi，通过PID调节算法，不断更新出最佳的参数，最终将调节后的电位图与健康状态的下丘脑电位进行对比以及做差，得出关于更好情况的刺激靶点以及电流优化情况。

5.2. 刺激靶点GPi/SNr

由于基底神经节回路模型中缺少由SNc所释放的多巴胺对整个回路的电位发放进行调节，导致下丘脑Th所产生的电位发放与健康状态下的Th相差较大，从而导致帕金森病症的产生，根据文献资料中的提示，在治疗帕金森病症时，通常通过电刺激对相关靶点，例如GPi/SNr、STN神经核团，以此分别凭借其对下丘脑Th的抑制以及对GPi/SNr的兴奋，最终导致更多的电刺激进入到下丘脑Th中。

我们在电信号进入靶点GPi/SNr前，通过施加高频正弦电信号，将更多的电刺激导入GPi/SNr神经核团内，以此达到兴奋下丘脑Th的最终目的，在施加高频正弦电信号的过程中，加入PID控制，使得帕金森病状的下丘脑Th电位发放在相关电信号的影响下，逐渐与健康状态下的下丘脑Th电位发放一致。通过这一调整过程的不断进行，可以得到相关调节电流的正弦函数，即为最终优化得到的最佳电流。

如下图9所示，通过刺激靶点GPi/SNr，通过PID算法自动控制，最终得到的下丘脑Th电位发放如图9(a)所示，通过对比发现，该放电电位与健康状态下的下丘脑放电电位基本一致，频率及振幅等参数也近乎相同，可以看出刺激靶点GPi/SNr之后得到的结果较为理想，根据两幅电位图进行对比后发现，如下图9(d)所示，其误差基本处于10 mV左右，该模型的结果表明可以在靶点GPi/SNr施加相应的电刺激，优化后的正弦电流如图9(c)所示，其函数公式为I = 4.8sin(50*π) + 10.4。

Figure 9. The result of stimulating the target GPi/SNr

图9. 刺激靶点GPi/SNr的结果

5.3. 刺激靶点STN

与刺激靶点GPi/SNr一致，我们选择另一个靶点对其施加相应的电刺激，通过PID自动控制，可以得到如下图10所示的结果，同理，将其中左上角的电位图与健康状态下的神经节回路电位发放图10(a)进行对比后，可以发现，其频率振幅也与其十分相近，通过做差分析，如图10(d)所示，可以得出这两幅电位图的误差为6 mV左右，相较于刺激靶点GPi/SNr中的误差值，可以得出相关结论：刺激靶点GPi/SNr和刺激STN都可以使得下丘脑Th中的电位发放与健康状态下的电位发放基本一致，但相较而言，在相同的控制原理的基础上，刺激STN可以获得更理想的结果，故最优的靶点为GPi/SNr。其中，控制完成后得到的正弦电流如图10(c)所示，其函数为I = 4.2sin(50π) + 10.2。

Figure 10. The result of stimulating the target STN

图10. 刺激靶点STN的结果

5.4. 其余最优靶点选择

根据PD状态图发现，刺激STN核团与GPi和核团后起到的效果都是减少对Thalamus的抑制作用，从而使得Motor Cortex兴奋，缓解帕金森症状，因此，在对其他靶点进行刺激后也应该起到减少对Thalamus的抑制作用，所以根据PD状态图中兴奋、抑制的关系可以发现，除了GPi核团与STN核团外，可以选择的刺激靶点就只有Striatum核团中的iMSN部分。

因此，我们将iMSN作为刺激靶点，我们通过PID自动控制器，将帕金森病态的神经元回路开放电位以健康状态开放电位为目标进行调节，最终得到的电位图如下图11(a)所示，健康状态下的Thalamus波动如图11(b)所示，刺激电流如图11(c)所示，可以发现，该电位图的频率与健康状态神经元回路的电位图频率相差较大，且误差图如图11(d)所示，误差基本处于17 mV左右，故表明在该PID控制算法下，刺激iMSN靶点的效果不理想，远不如刺激GPi/SNr或STN所计算出的下丘脑Th核团电位。

Figure 11. The result of stimulating the target iMSN

图11. 刺激靶点iMSN的结果

最终可以得出相关结论，在直接通路或间接通路的神经通路中，脑深部电刺激治疗不存在其他最优刺激靶点。

6. 结论

本文通过搭建神经元Hodgkin-Huxley模型，经过直流与交流电刺激对神经元模型进行有效验证，通过搭建化学突触模型中的兴奋型突触模型——谷氨酸能突触(AMPA)以及抑制型突触模型——GABA能突触(GABA)，将单体神经元连接成为复杂的神经核团，用于搭建基底神经节神经回路的理论模型，依据基底神经节内部神经核团连接方式进行一致性连接，最终搭建出整个治疗帕金森病的仿真模型。

该模型通过众多的参数共同调节，将神经团块中的神经元数量定为6个后，计算了基底神经节内部神经元的电位发放情况，得出与文献中的内部电位一致的结果，肯定了该模型的真实性，通过去掉黑质SNc，观察到正常状态与帕金森病态中下丘脑Th的电位发放情况，在整个仿真过程中，同时提取了GPi、GPe、STN以及Thalamus四个神经核团的电位情况，且通过数值确定了正常状态与帕金森病态两种情况下的电位发放特征指标。最后通过在直接通路以及间接通路中寻找其余可靠刺激靶点，确定了iMSN可作为允许的刺激靶点，但从最后优化后的电位发放情况中可以发现，靶点iMSN的刺激效果远不如靶点STN与GPi，其误差大于17 mV。