基于神经元Hodgkin-Huxley模型的脑深部电刺激治疗帕金森病的建模研究
帕金森病作为一种全球常见的精神退行性疾病，日趋成为中老年人正常生活的一大威胁。目前缓解帕金森病症状的治疗方法主要有三种：药物治疗、手术治疗和脑深度刺激。脑深度刺激作为一种副作用小、安全性高的新方法受到越来越多患者以及研究者的重视。脑深度刺激治疗帕金森病的机理来自基底神经节。本文使用神经元Hodgkin-Huxley模型为基础，分别构建直接通路，间接通路，健康回路以及帕金森病回路，并对其电位发放状态进行研究。之后对脑深度刺激治疗帕金森病的影响因素进行分析，考虑了不同靶点以及刺激参数对帕金森病治疗效果的影响。最后对直接/间接通路中除 STN和 GPi外是否存在其他最优电刺激靶点进行讨论。
问题1：从Hodgkin-Huxley模型出发，通过数值求解发现当神经细胞受到直流刺激时，如果刺激较小，产生的作用将很快衰减，细胞重回极化状态；足够大的直流电会使膜电位产生周期放电，放电波峰产生的时间、速度、高度和间隔与外界刺激的强度有直接关系；当施加交流刺激时，反应电位和刺激电流幅值、周期及信号形式均存在相关关系。
问题2：基于单个神经元的H-H模型，考虑化学突触，建立了神经团模型；参考已有研究，选取了常见的神经团块信号传递关系，并对交、直流电流刺激作用下的直接通路和间接通路上的神经团块进行放电分析得到电位响应特征曲线。
问题 3：健康回路和 PD 回路之间的显著差别在于 BG 黑质核团多巴胺神经元是否可以正常工作，本问中分别建立了两种回路模型，探究了 SNc 在健康回路中的调控作用，对比发现SNc对于神经回路的正常工作起到重要作用。发现在 PD 回路中，SNc 作用减弱，不能发挥正常作用。
问题 4：本问对于帕金森症的 DBS 治疗靶点进行了对比研究，选取 STN 和 GPi/SNr两个重要位置给予刺激，比较两个靶点的治疗效果，探究了治疗效果与刺激电流参数之间的相关关系。
问题 5：为了探究是否存在对帕金森症治疗效果更佳的神经刺激靶点，分析了健康回路中SNc的作用，并选取GPe作为靶点分析了GPe对 STN的抑制作用，评估了刺激 GPe靶点作为治疗手段的可行性。
关键词：脑深度刺激、Hodgkin-Huxley模型、帕金森病回路、基底神经节
1 问题重述 
1.1 问题背景
帕金森病（Parkinson disease，PD）是仅次于阿尔茨海默病的第二常见的神经退行性疾病，PD 多发生于老年人,且随着年龄增长发病率增高。年龄≥65岁的老年人 PD 发病率约1%,年龄≥85岁的老年人约5%[1]。超过 90% 的 PD 患者是特发性的，没有明确的病因。
PD 主要病理特点为中脑黑质多巴胺能神经元严重缺失和纹状体多巴胺神经递质减少[2]。
临床表现的特征是静止性震颤，肌强直，运动迟缓，姿势步态障碍等运动症状。同时患
者可伴有抑郁、便秘和睡眠障碍等非运动症状。帕金森病的诊断主要依靠病史、临床症
状及体征，一般的辅助检查多无异常改变。目前缓解帕金森病症状的治疗方法主要有：药物治疗、手术治疗和脑深部刺激（DBS）三种[3]。药物治疗用于早期帕金森疾病，手术治疗适用性较差且切除后不可逆。DBS 通过精确定位，选取脑内特定的靶点植入刺激电极，通过输入高频电刺激，改变相应核团的兴奋性，达到改善治疗帕金森病症状的效果。DBS 治疗帕金森病的靶点包括丘脑底核（STN）和苍白球内侧核（GPi/SNc）的脑深部电刺激等[4][5]。
脑深部电刺激治疗帕金森病的机理来自基底神经节（Basal ganglia，BG），BG结构和丘脑-皮层（Thalamus-Cortex）结构如图 1.1。与基底神经节相关的神经核团（Cortex，Striatum/dMSN/iMSN，SNc，GPe，GPi/SNc，STN，Thalamus）内部包括大量的神经元，核团之间相互连接，发生神经信息传递，经典的基底神经节神经团块结构中，神经信息的传导包括两条相互平行的信号传递通路[6]：直接通路（Cortex→Str→GPi/SNr）和间接通路（Cortex→Str→GPe→STN→GPi/SNr）。
图1.1 基底神经节内部神经核团连接（箭头方向代表兴奋，方块方向表示抑制，绿色路径
是直接通路，红色路径是间接通路，蓝色路径是多巴胺）
直接通路和间接通路分别通过丘脑兴奋或抑制运动皮层（图 1.2）。虽然帕金森病的病理机制目前仍不十分清楚，但临床主流的观点是[7]：基底神经节黑质核团（SNc）多巴胺能神经元的退化状体
STN 直接参与基底神经节中直接通路的运动发起与间接通路的运动抑制之间的平衡。STN-DBS对改善PD运动障碍神经机制，优化DBS参数对PD治疗策略等具有重要意义。同样，GPi-DBS对PD运动控制也产生直接影响。
图1.2 基底神经节内部神经核团的健康状态和PD状态
显然，帕金森病的脑深部刺激治疗的运动控制机理十分复杂，研究脑深部电刺激治疗的最优靶点选择和 DBS 参数（刺激强度和刺激频率等）优化，可以从基底神经节黑质核团的损害入手，模型分析基底神经节中神经元电位发放的特征指标；或者研究基底神经节中直接通路与间接通路之间的信息传递平衡，也许是解决 PD难题的突破口[8]。
1.2 问题提出
基于上述研究背景，本文需研究和解决以下问题：
问题1：利用给出的神经元Hodgkin-Huxley模型（附件1），数值模拟外界刺激（包括直流刺激和交流刺激）情况下，单个神经元的电位发放情况，并给出神经元电位发放的特征指标。
问题2：根据问题 1的神经元Hodgkin-Huxley模型，结合附件1中神经元之间的突触连接理论，建立基底神经节神经回路的理论模型，计算基底神经节内部神经元的电位发放（每个神经团块可以简化为5-10个神经元）。
问题 3：根据建立的基底神经节回路模型，理论分析正常状态（图 1.2 中 Healthy 回路）和帕金森病态（图 1.2 的 PD 回路中去掉黑质 SNc）基底神经节回路电位发放的特征指标。
问题 4：利用建立的基底神经节回路模型，对帕金森病态的基底神经节靶点添加高频电刺激，可以模拟脑深部电刺激治疗帕金森病的状态。请模型确定最佳刺激靶点，是刺激靶点STN，还是刺激靶点 GPi；请模型优化刺激的参数，如电刺激强度，电刺激频率和电刺激模式等。 
问题 5：在直接通路的神经通路中，或者间接通路的神经通路中，模型回答脑深部电刺激治疗是否存在其它最优电刺激靶点。
2 模型假设
(1) 每个神经元的特征相同，即数学模型中每个神经元的常数都保持不变。
(2) 神经元所处的外坏境相同，即在分析神经元之间的相互作用时，不考虑环境造成 的电位差异，只考虑神经元产生的信号引起的电位变化。
(3) 神经元之间只考虑兴奋与抑制两种关系，不考虑其他关系。
(4) 对帕金森病态的基底神经节靶点只分别施加直流刺激与交流刺激。
3 符号说明及中英文对照表
论文中用到的一些符号及其意义汇总如下表：
表3.1 符号说明表
| 符号 | 意义 |
神经元膜电容
神经元膜电位
钠离子通道特性
钾离子通道特性
泄露电流特性
钠离子关于细胞膜的电导系数最大值
钾离子关于细胞膜的电导系数最大值
泄露电流关于细胞膜的电导系数最大值
钠离子反向电压
钾离子反向电压
泄露电流反向电压
开关函数
外界神经元的刺激影响
神经元之间的化学突触电流
兴奋突触电流
抑制突触电流 
最大电导
受体开放状态的比例
突触前电压
突触后电压
逆转电位
突触模型常系数
表3.2 中英文对照表
| 中文术语 | 英文术语 |
| 帕金森病 脑深度刺激 基底神经节 丘脑底核 苍白球外段 苍白球内段 黑质 黑质致密部 纹状体 直接通路 间接通路 超直接通路 丘脑 皮层 | Parkinson disease，PD deep brainstimulation，DBS basal ganglia，BG subthalamicnucleus，STN external globus pallidus，GPe Striatum，Str. direct pathways indirect pathways hyperdirect pathways Thalamus cortex |
4 问题分析
本文从单个神经元的 H-H 模型出发，通过突触连接理论构建单神经团模型，考虑到神经团之间的兴奋、抑制关系，会产生相互平行的信号传递通路，即直接通路和间接通路作用机制。对于健康回路和帕金森病态的神经回路，会产生信号传导通路的差异，分别构建正常回路和 PD 回路数值模型并施加直流、交流刺激作用于不同的神经团块靶点，总结分析不同刺激条件下的特征指标（包括但不限于动作电位-时间关系图的振幅、频率、平均频率、静息间隔、激活时间、峰峰间距、簇发放周期等）。比较刺激不同靶点时 PD回路的特征指标，取修正效果相对最佳的刺激靶点和电流作为治疗方案。
图4.1 问题求解流程图
5 问题1建模与求解 
5.1 基于电导的神经元 Hodgkin-Huxley 模型
Hodgkin-Huxley 中的神经元模型信号作用关系如下图所示：
图5.1 H-H模型信号传递关系图
神经系统中的神经细胞通过突触连接，在静息时细胞膜内外将会存在由于离子浓度不同而引起的电位差。神经系统受到外界刺激时，膜电位产生节律放电。神经元电学模型有很多种，比如 Integrate-and-Fire 模型，FitzHugh-Nagumo 模型，Rall 电缆模型，Hodgkin-Huxley模型等，本文选用 Hodgkin-Huxley 模型[9]对神经元电位变化进行模拟。
Hodgkin-Huxley 模型控制方程如下：
(5.1)
图5.2 Hodgkin-Huxley 神经元轴突的等效电路和电位发放
在H-H模型中，
| (1) (2) (3) |
| 表示神经元膜电容， |
| ； | |||
|
| 表示神经元的膜电位； | ||||||
| ， | ， | 分别描述细胞膜内外的钠离子通道、钾离子通道和泄露电流特 | |||||
性;
| (4) 常系数 |
| 对应钠离子关于细胞膜的电导系数的最大值， | ||
|
| 对 应 钾 离 子 关 于 细 胞 膜 的 电 导 系 数 的 最 大 值 ， | |||
|
| 对应泄漏电流关于细胞膜的电导系数的最大值; | |||
| (5) | ， | ， | ，分别对应钠离子、钾离子和 | |
泄露电流的反向电压;
(6) 模型中的开关函数分别是：
| (7) (8) |
| (5.2) | |
|
| 对应外界对神经元的刺激影响； | ||
|
| 对应神经元之间的化学突触电流，一般包含兴奋突触和抑制突触两类； | ||
(9) 神经元的电位发放就是神经元的动作电位-时间关系图。
5.2 数值方法介绍
确定选用的神经元模型之后需要选择合适的数学方法对 H-H 模型中的常微分方程组进行求解。龙格-库塔法(Runge-Kutta)是一类求解非线性常微分方程的重要迭代算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。
一阶常微分方程可以写作：
|
| (5.3) |
借助差分可以表达为:
|
| (5.4) (5.5) (5.6) | ||||
| 推出（近似等于，极限为 |
| ） | |||
|
| |||||
| 另外根据微分中值定理，存在 |
| ,使得 | |||
|
| |||||
其中
|
| (5.7) |
称为平均斜率，龙格库塔方法就是求得 的一种算法。 
在此基础上，经过复杂的数学推导，可以得出截断误差为
的四阶龙格库塔公式：
(5.8)
在数值计算软件 Matlab 中内建有成熟的基于龙格-库塔法的数值求解函数 ode23、ode45、ode15 等，其在求解原理上并无太大差异，仅仅存在求解精度及速度上的差别。考虑到问题的复杂性，为了加快计算进度，采用 ode15s 方法进行求解，通过与相关文献[9]上的神经元放电曲线对比并无差异，可见精度满足使用要求，本文后续计算若非特别说明，均为ode15s数值求解结果。
5.3 结果分析
通过 Matlab 数值求解 Hodgkin-Huxley 神经元模型，在直流刺激条件下，单个神经元的电位发放情况如下图所示（以放电峰的个数为基准调整刺激电流大小，以期得到临界放电状态改变时对应的电流特征）：
5.3.1 神经元电位变化曲线-直流刺激
图5.3 直流刺激
图5.4 直流刺激
代码
42. [T,y]=ode15s(@(T,y) NeuronDC(T,y,C,gNa,gK,gL,VNa,VK,VL,I_extern
al(i),I_synapse),Region,Origin,opts);
43. % Time.Ti=T;
44. eval(['Time(1).T',num2str(i),'=T']);
45. % Potential.Pi=y;
46. eval(['Potential(1).P',num2str(i),'=y']);
| 47. 48. % Plot 49. 50. % single plot 51. figure (i) 52. 53. set(gcf,'Units','normalized','Position',[0 0 1 1]); 54. % plot(Time.Ti,Potential(:,1).Pi,'k','LineWidth',3.5); 55. eval(['plot(Time(1).T',num2str(i),',Potential(1).P',num2str(i), '(:,1),''k'',''LineWidth'',3.5)']); 56. set(gca,'FontSize',20); 57. set(gca, 'LooseInset', [0,0,0,0]); 58. xlim auto 59. ylim auto 60. title=strcat('I_{dc}=',num2str(I_external(i)),' \muA/cm^{2}'); 61. xlabel({'Time t/ms',title},'FontSize',25) 62. ylabel('Potential V/mV','FontSize',25) 63. 64. Picture=strcat('DC',num2str(i)); 65. print(gcf,Picture,'-dpng','-r900'); 66. 67. end 68. 69. 70. 71. % AC 72. 73. [rows,~]=size(w); 74. Time=struct('T1',{},'T2',{},'T3',{},'T4',{},'T5',{},'T6',{}); 75. Potential=struct('P1',{},'P2',{},'P3',{},'P4',{},'P5',{},'P6',{}); 76. for i=1:6 77. 78. I_external3=10*cos(w(i)*It); 79. [T,y]=ode15s(@(T,y) NeuronAC(T,y,C,gNa,gK,gL,VNa,VK,VL,It,I_ext ernal3,I_synapse),Region,Origin,opts); 80. % Time.Ti=T; 81. eval(['Time(1).T',num2str(i),'=T']); 82. % Potential.Pi=y; 83. eval(['Potential(1).P',num2str(i),'=y']); 84. 85. % Plot 86. 87. % single plot 88. figure (i) 89. set(gcf,'Units','normalized','Position',[0 0 1 1]); 90. % plot(Time.Ti,Potential(:,1).Pi,'k','LineWidth',3.5); 91. eval(['plot(Time(1).T',num2str(i),',Potential(1).P',num2str(i), '(:,1),''k'',''LineWidth'',3.5)']); 92. set(gca,'FontSize',20); 93. set(gca, 'LooseInset', [0,0,0,0]); |