静态直流电场作用下神经元放电动力学分析

伊国胜，王 江，魏熙乐，邓 斌，韩春晓

(1. 天津大学电气与自动化工程学院，天津 300072；

2. 天津职业技术师范大学自动化与电气工程学院，天津 300222)

无创式脑刺激技术如今已经被广泛应用到治疗各种神经和精神疾病中．其中常用的一种神经调制方法就是采用静态直流电场刺激中央神经系统，例如经颅直流电刺激(transcranial direct current stimulation，tDCS)．它采用弱直流电流，通过安置在脑部的两个电极形成静态直流电场，进而直接刺激脑皮层相应的区域[1]．研究表明，tDCS能够引起脑皮层神经可塑性的改变[2]，增强运动和认知功能[2]，同时能够调节工作记忆能力[3]．此外，临床 tDCS对多种神经精神疾病具有很好的疗效，如慢性疼痛综合症、中风、癫痫、重症抑郁、帕金森和渴望紊乱症等[4]．动物实验表明，阳极 tDCS通过引起静息态细胞膜电动势去极化而增加神经元放电速率；阴极 tDCS具有相反的效果，它通过引起静息态细胞膜电动势超极化而减小神经元放电速率[5]．但是 tDCS对脑功能的具体调制机制至今还是未知，同时对于静态直流电场是如何影响神经元活动的也缺乏理论研究和机制分析．

1977年，Mackey和 Glass等提出了动力学疾病的概念，指出当控制生理系统的参数超出正常范围时，便会出现不正常的动力学现象[6]．可见，神经系统动力学过程和其功能有着密切的关系．当控制系统动力过程的关键参数出现异常时，会导致系统动力学特性的异常，因而会影响神经系统的功能．研究发现，许多精神疾病如帕金森、阿兹海默症、抑郁症和强迫性精神障碍等，都与神经元的异常放电有关[6-9]．Devos等[10]提出脑调制技术的机制可能是引起病症的异常神经元放电节律能通过刺激变成正常的生理节律．Meissner等通过实验研究发现 DBS能够阻断底丘脑核(STN)中不正常的低频放电，证实了上述假说[6]．因此，从动力学角度研究静态直流电场作用下神经元的放电模式，可为探索神经调试技术如何改善神经系统动力学过程提供一定的理论基础．

探索静态直流电场对脑功能和神经活动的影响一直是研究学者和医学治疗所关心的问题．实验研究发现，外加电场或电流能够影响许多神经系统的兴奋性，如海龟的小脑、海马齿状细胞、锥体细胞和皮层神经元等[11]．1996年，Gluckman等研究发现静态直流电场加到大脑组织能够暂时抑制癫痫发作[11]．Durand和Bikson等通过对多种癫痫动物模型和病人的实验研究，发现低幅值的直流均匀电场和电流刺激产生的局部电场在控制脑切片的癫痫样神经元活动上同样有效[11]．此后，Bikson及其合作者以大鼠海马趾的切片为对象，通过实验研究了直流外电场对神经元兴奋性的影响，结果表明直流电场改变了顺向刺激引起的动作电位阈值；较大电场能够触发神经元兴奋而进入癫痫状态，并可以长时间地改变神经元的兴奋性[12]．但是，上述大部分研究结果在理论研究中很难复现，并且很难形成一个统一的理论框架．此外，其多数以实验研究为主，缺乏理论研究和机制分析．

本文从生物物理模型角度研究静态直流电场对脑活动和脑功能的影响．首先在Prescott模型[13]基础上对其进行改进，用以描述静态直流电场作用下单神经元的动力学行为；然后分析了静态直流电场作用下Hodgkin定义的3类神经元[14]的放电模式；最后从动力学角度研究分析了各种放电模式产生的起始动力学机制．

1 模型和方法

1.1 外电场对神经元细胞膜电压的影响

当人脑处于静态直流电场中时，静态直流电场会在神经元细胞膜上产生一个去极化感应电压ΔV ．Schwan研究发现，当神经元细胞是规则球形时，幅值为E的静态电场与其产生的细胞膜去极化感应电压ΔV之间的关系[6]为

式中：R为神经元细胞半径，大约为10µm；θ为外部电场轴线和神经元细胞轴线之间的角度，当神经元完全处在感应电场中时，cosθ=1．但是实际神经元并不是规则的球形，它们具有大量的树突和轴突．因此，由感应电场引起的移动电荷会在细胞膜中累积，这样产生的细胞膜去极化感应电压就会增大．因此，考虑细胞几何形状和极化率时，应该引入极化长度λ[15]．极化长度λ的数值可以通过记录不同外部电场作用下神经元细胞膜电压来测得[15]．此时，稳定电场E和其产生的细胞膜去极化感应电压ΔV之间的关系为

式中τ为 Maxwell-Wagner时间常数，其物理含义是电荷以何种速度在细胞膜上累积[16]．根据式(2)可知，当感应电场为幅值是E的静态直流电场时，在细胞膜上产生的细胞膜去极化感应电压ΔV是直流的，即

为了研究方便，本文取极化长度λ=1mm，即静态直流电场等于其产生的细胞膜去极化感应电压．

1.2 Hodgkin 3类神经元兴奋性及Prescott模型

Hodgkin通过研究不同幅值直流电流对神经元兴奋性的影响，根据神经元放电频率和刺激电流之间的关系将神经元分成 3类[14]．Ⅰ类神经元在不同幅值电流刺激下会产生任意低频率的动作电位，其放电f-I曲线是连续的；Ⅱ类神经元不能产生任意低频率的动作电位，其放电 f-I曲线是不连续的；Ⅲ类神经元在不同幅值电流刺激下只能产生一个动作电位，不能产生重复放电，其放电 f-I曲线是没有定义的．因为计算神经元放电速率至少需要两个动作电位，但是Ⅲ类神经元至多只能产生一个动作电位．

Prescott等在 ML模型基础上，通过大量的实验研究构建了一个二维模型．虽然此模型简单，但是通过改变模型中一个参数βw，便能同时复现 Hodgkin定义的3类神经元兴奋性，并且可以解释其产生的动力学基础和生物物理基础．Prescott模型的具体形式[13]为

该模型含有两个变量：一个是快速激活变量V，代表细胞膜电压，控制钠离子电流INa的瞬态激活；另一个是慢速恢复变量w，是一个电压的函数，代表钾离子电流IK的慢激活变量，控制钾离子电流的激活．IS代表刺激电流，IL是漏电流．τw是慢变量w的时间常数，βw是模型的关键参数，模型中参数取值如表1所示．

在模型(4)中仅改变参数βw就能复现 Hodgkin定义的3类神经元兴奋性．图1中，当 βw=0 mV 时，神经元呈现Ⅰ类兴奋性，其 f-I曲线随着刺激电流增加是连续变化的；βw=−1 3mV时神经元呈现Ⅱ类兴奋性，其放电 f-I曲线是不连续的；βw=−2 1mV时神经元呈现Ⅲ类兴奋性，其放电频率始终为0,Hz. Prescott等利用此最小模型研究了 3类神经元产生的动力学基础．研究发现 3类兴奋性的产生是阈下动力学不匹配的反向电流之间非线性竞争的结果[13]．Ⅰ类兴奋性是通过极限环上的鞍点-节点分岔产生的；Ⅱ类兴奋性是通过亚临界霍夫分岔产生的；Ⅲ类兴奋性是通过穿越准分界线(QSC)产生单峰放电．图2是3类神经元在直流电流作用下的放电起始动态机制的相平面分析和分岔图[13]．

图1 Hodgkin 3类神经元放电f-I曲线Fig.1 Frequency-current(f-I)curve of Hodgkin’s three classes of neurons

1.3 静态直流电场作用下神经元模型

本文采用 Prescott模型建模研究脑皮层神经元在静态直流电场作用下的动态响应．当直流电场作用在神经元细胞上时，相当于其产生的细胞膜去极化感应电压ΔV作用在 Prescott模型上，将ΔV看作细胞膜电压V的一个外加扰动[6]．因此，可通过引入一附加项ΔV对静态直流电场 E的作用进行建模[6]．此时，神经元在静态直流电场的作用下，由于细胞膜的介质特性，膜两侧将在 V的基础上叠加一个跨膜的感应电压分量．故此时神经元总的膜电位为

在静态直流电场作用下，外电场参数的引入相当于在细胞膜上叠加了一个直流跨膜电动势，并未引起模型基本结构的变化，而只是同时改变了钠离子电流、钾离子电流和漏电流各自的反电动势．因此，描述静态电场作用下神经元动态响应的Prescott模型为

极化长度λ=1mm，根据式(3)可知

即静态直流电场等于其产生的细胞膜去极化感应电压．由于本文只研究静态直流电场E对神经元放电模式的影响，故刺激电流 IS=0µA/cm2．模型(6)中部分参数取值见表1．

表1 模型(4)和模型(6)参数取值Tab.1 Parameters in model(4) and model(6)

本文中各种放电模式是通过 Matlab软件，采用龙格-库塔法以 0.01,ms的步长进行仿真得到的．文中各种相平面图和分岔图是利用 XPP软件，采用Euler方法以0.1,ms的步长进行仿真得到的．初始值选为当不存在静态直流电场刺激时，神经元最终稳定的静息值，即 V (0)=−6 4.9 mV 和w(0)=0．仿真时间持续4,000,ms，这样系统每个变量都能达到稳定状态.

图2 Hodgkin 3类兴奋性产生的相平面分析和分岔图Fig.2 Phase plane and bifurcation analysis of Hodgkin′s three classes of neurons

2 静态直流电场作用下神经元放电模式动力学分析

2.1 静态直流电场作用下神经元放电模式

3类神经元在不同静态直流电场作用下放电模式如图3～图5所示．图6是3类神经元放电频率f随静态直流电场幅值E变化曲线．

由图3～图5可见，当静态直流电场为正值时，3类神经元均处于静息状态，不产生动作电位，且稳态细胞膜电压值随静态直流电场幅值E的增大而增加．当静态直流电场为负值时，3类神经元会出现周期性放电现象．

图3 静态直流电场作用下Ⅰ类神经元放电模式Fig.3 Firing patterns of classⅠneuron exposed to static electric field

图5 静态直流电场作用下Ⅲ类神经元放电模式Fig.5 Firing patterns of classⅢneuron exposed to static electric field

图6 静态直流电场作用下3类神经元放电f-E曲线Fig.6 f-E curves of three classes of neurons exposed to static electric field

在静态直流电场作用下，Ⅰ类和Ⅱ类神经元放电模式变化过程类似．当E变化到某一临界值时，Ⅰ类和Ⅱ类神经元出现周期性放电行为．其放电频率f随着静态直流电场幅值 E的变化呈“抛物线”状，如图 6(a)和 6(b)所示．当神经元放电频率增加到最大值后，Ⅰ类和Ⅱ类神经元放电频率会随幅值 E的增加而逐渐减小最终变为 0．神经元膜电压会最终稳定在阈上膜电位处，如图 3(a)和图 4(a)中 E=−50 mV所示．说明此时刺激超过了神经元能够承受的正常生理范围，本文不对其进行分析．与Ⅰ类和Ⅱ类神经元不同的是，Ⅲ类神经元在静态直流电场E在[−2 2.29 mV,−20.62 mV]范围内会出现单峰放电行为，如图5(e)E=−2 2 mV 所示．3类神经元所能达到的最大的放电速率是不同的，Ⅰ类神经元的最大，Ⅱ类神经元次之，Ⅲ类神经元最小，如图6所示；其次，Ⅰ类神经元放电场 f-E曲线左侧是不连续的，右侧是连续变化的；Ⅱ类神经元放电 f-E曲线两侧均是不连续变化的；Ⅲ类神经元放电 f-E曲线左侧是连续变化的，右侧是不连续的．

2.2 静态直流电场作用下神经元放电起始动态机制相平面分析与分岔类型分析

为了揭示静态直流电场作用下 3类神经元产生各种放电模式所对应的放电起始动态机制，需要分析3类神经元不同放电模式产生的动力学基础．因为模型(6)是二维的，故神经元的放电起始动态机制可由相平面内V零线和w零线的相交方式来解释[17-18].零线是相空间中一切使变量的导数为 0的点相互连接而成的曲线[17]．

图 7解释了Ⅰ类神经元在静态直流电场E作用下，放电起始动态机制的相平面分析与分岔类型分析．当不存在直流电场时，V零线和w零线有3个交点，如图 7(a)中 E=0 mV 所示．此时系统有 3个平衡点：1个稳定的平衡点和 2个不稳定的平衡点．系统所有状态最终都收敛到这个稳定的平衡态，神经元不会出现动作电位，对应于静息状态．随着施加负向直流电场 E，V零线不断向右上方移动．最左侧的稳定平衡点和其相邻的不稳定平衡点逐渐吸引，当 E增加到一个临界值时，二者结合在一起被破坏然后消失，同时出现一个稳定的极限环，此时对应周期放电产生．图 7(a)中箭头表示相空间中膜电压轨迹运动方向．此过程中，平衡点被破坏和周期性放电的产生是通过极限环上鞍点-节点分岔完成的，如图 7(b)所示．随着负向静态电场幅值的继续增大，当静态直流电场 E 在[-40.7,mV，-16.65,mV]范围内，V零线和w零线只交于一个不稳定平衡点，并且稳定的极限环均存在，故此时神经元继续保持周期放电状态，如图7(a)中E=−38 mV ．当E＜−40.70 mV时，V零线和w零线交点由一个不稳定平衡点变为稳定的平衡点，与此同时稳定极限环消失，如图 7(a)中E=−50 mV所示．此时系统所有状态最终都收敛到这个稳定的平衡态，神经元不会出现周期性放电现象．此平衡点对应的膜电压属于阈上范围，故神经元最终稳定在阈上的膜电位处，如图 3(a)所示．此时，直流电场幅值过大，超过了神经元所能够承受的正常生理范围．此过程涉及到的极限环消失和稳定平衡点的出现是通过亚临界霍夫分岔完成的，如图7(b)所示．

图7 静态直流电场作用下Ⅰ类神经元放电起始动态机制Fig.7 Dynamical mechanism of action potential initiation exhibited by classⅠneuron under the effects of static DC electric field

图8解释了Ⅱ类神经元在静态直流电场E作用下，放电起始动态机制的相平面分析与分岔类型分析．当不存在直流电场时，V零线和w零线交于一个稳定的平衡点，如图8(a)E=0 mV 所示．此时系统所有状态最终都收敛到这个稳定的平衡态，神经元不会出现动作电位，对应于静息状态．随着负向直流电场E幅值的不断增加，V零线不断向右上方移动．当 E增加到一个临界值时，稳定平衡点会通过亚临界霍夫分岔变为一个不稳定的平衡点和一个稳定的极限环，对应周期放电产生．图8(a)中箭头表示相空间中膜电压轨迹运动方向．随着负向静态电场幅值的继续增大，当静态直流电场 E在[-45.79,mV，-17.63,mV] 范围内，V零线和w零线只交于一个不稳定平衡点，并且稳定的极限环始终存在，故此时神经元继续保持周期放电状态，如图 8(a)中 E=−4 2 mV所示．当E＜−45.79 mV时，V零线和w零线交点会通过亚临界霍夫分岔由一个不稳定平衡点变为稳定的平衡点，与此同时稳定的极限环消失，如图 8(a)中E=−5 0 mV所示．此时神经元不会出现周期性放电现象，并且膜电压最终稳定在阈上的膜电位处，如图4(a)所示．

图8 静态直流电场作用下Ⅱ类神经元放电起始动态机制Fig.8 Dynamical mechanism of action potential initiation exhibited by classⅡneuron under the effects of static DC electric field

图 9解释了Ⅲ类神经元在静态直流电场E作用下的起始动态机制相平面分析与分岔类型分析．当不存在直流电场时，V零线和w零线交于一个稳定的平衡点，对应于静息状态，神经元不会出现动作电位，如图 9(a)E=0 mV所示．当静态直流电场 E在[−20.62 mV,−22.29 mV]范围内时，神经元通过穿越准分界线(QSC)产生单峰放电，如图 9(b)所示．此时在静态电场的作用下，处于准分界线下方的阈下膜电压轨线会绕过准分界线的最右侧头部“*”处，然后回到稳定平衡点，从而产生一个单峰放电．图 9(a)和图9(b)中箭头表示相空间中膜电压轨迹运动方向．当E继续增加到一个临界值时，稳定平衡点会通过亚临界霍夫分岔变为一个不稳定的平衡点和一个稳定的极限环，对应周期放电产生．随着负向静态电场幅值的继续增大，当静态直流电场E在[−46.16 mV,−22.30 mV]范围内时，V零线和w零线只交于一个不稳定平衡点，并且稳定的极限环始终存在，此时神经元继续保持周期放电状态，如图 9(a)中E=−45 mV所示．当E＜−46.16 mV时，V零线和w零线交点会通过鞍点-节点分岔由一个不稳定平衡点变为3个平衡点(分别是1个稳定的平衡点和2个不稳定的平衡点)，与此同时稳定的极限环消失，如图9(a)中E=−50 mV 所示．此时系统所有状态最终都收敛到这个稳定的平衡态，神经元不会出现周期性放电现象．此平衡点对应的膜电压属于阈上范围，故神经元最终稳定在阈上的膜电位处，如图5(a)中所示．

图9 静态直流电场作用下，Ⅲ类神经元放电起始动态机制Fig.9 Dynamical mechanism of action potential initiation exhibited by classⅢneuron under the effects of static DC electric field

当产生极限环上鞍点-节点分岔时，随着静态直流电场幅值的变化，细胞膜电压轨线(如图 7（a)中E=−17 mV 对应的虚线）能够缓慢地通过V零线和w零线之间的狭小空间，因此神经元能够缓慢地放电，故此时神经元放电 f-E曲线是连续的，如图 6(a)右侧连续部分和图 6(c)左侧连续部分．当产生亚临界霍夫分岔时，由于不存在极限环上鞍点-节点分岔时V零线和w零线之间狭小空间的约束，因此神经元不能够进行慢速放电，故神经元放电f-E曲线是不连续的，如图 6(a)左侧不连续部分和图 6(c)右侧不连续部分及图6(b)所示．

3 结 语

本文通过对 Prescott二维最小神经元模型进行改进，研究了静态直流外电场对Hodgkin定义的3类神经元放电模式的影响．研究发现，在负向直流电场作用下，Ⅰ类神经元会通过极限环上鞍点-节点分岔使稳定平衡点变成稳定极限环，从而产生周期性放电，然后通过亚临界霍夫分岔使稳定极限环转变成稳定平衡点，结束周期性放电；Ⅱ类神经元会通过亚临界霍夫分岔使稳定平衡点变成稳定极限环，从而产生周期性放电，然后再通过亚临界霍夫分岔使稳定极限环转变成稳定平衡点，结束周期性放电；Ⅲ类神经元会首先通过穿越准分界线产生单峰放电，此时放电 f-E曲线是没有定义的，然后通过亚临界霍夫分岔使稳定平衡点变成稳定极限环，从而产生周期性放电，最后通过极限环上鞍点-节点分岔使稳定极限环转变成稳定平衡点，结束周期性放电．随着负向直流电场幅值的增加，3类神经元放电频率先增加后减小，均呈“抛物线”形状．

人脑暴露在静态直流外电场下，根据本文采用的生物物理模型研究结果可知，其可以改变神经元的放电特性，进而对神经元活动产生影响．本文的研究结果证实了静态直流电场可以用在神经系统疾病治疗领域，用以改变病态的神经元放电节律．本文的研究结果可以为探究 tDCS的治疗机理提供一定的理论基础，从而更好地改善这种基于静态直流电场刺激的治疗方案．此外，本文的研究结果不仅可以解释静态直流外电场作用下神经元放电特性的实验现象，还可以用来完善外电场作用下神经元模型，从而为进一步研究各种无创式电磁刺激疗法(如 TMS、rTMS和tDCS等)和电磁辐射对人体的影响提供一定的理论基础．

［1］ Nitsche M A，Cohen L G，Wassermann E M，et al.Transcranial direct current stimulation：State of the art 2008[J].Brain Stimulation，2008，1(3)：206-223.

［2］ Datta A，Bikson M，Fregni F. Transcranial direct current stimulation in patients with skull defects and skull plates：High-resolution computational FEM study of factors altering cortical current flow[J].Neuroimage，2010，52(4)：1268-1278.

［3］ Ohn S H，Park C I，Yoo W K，et al. Time-dependent effect of transcranial direct current stimulation on the enhancement of working memory[J].Neuroreport，2008，19(1)：43-47.

［4］ Wahner T，Valero-Cabre A，Pascual-Leone A. Noninvasive human brain stimulation[J].Annu Rev Biomed Eng，2007，9：527-565.

［5］ Oostendorp T F，Hengeveld Y A，Wolters C H，et al.Modeling transcranial DC stimulation[C]//The30th Annual International IEEE EMBS Conference.Vancouver，Canada，2008：4226-4229.

［6］ Modolo J，Thomas A W，Stodilka R Z，et al. Modulation of neuronal activity with extremely low-frequency magnetic fields：Insights from biophysical modeling[C] //The2010IEEE5th International Conference on Bio-Inspired Computing Theories and Applications(BICTA).Beijing，China，2010：1356-1364.

［7］ Cziger B，Csikos D，Hidasi Z，et al. Quantitative EEG in early Alzheimer disease patients—Power spectrum and complexity features[J].International Journal of Psycho-physiology，2008，68(1)：75-90.

［8］ Bolwig T G，Hansen E S，Hansen A，et al. Toward a better understanding of the pathophysiology of OCD SSRI responders：QEEG source localization[J].Acta Psychiatr Scand，2007，115(3)：237-242.

［9］ Smeal R M，Ermentrout G B，White J A. Phaseresponse curves and synchronized neural networks[J].Philos Trans R Soc Lond B Biol Sci， 2010，365(1551)：2407-2422.

［10］ Devos D，Lefebvre L. Effect of deep brain stimulation and L-Dopa on electrocortical rhythms related to movement in Parkinson's disease[J].Prog Brain Res，2006，159：331-349.

［11］ 车艳秋. 外电场下神经元的分岔研究[D]. 天津：天津大学电气与自动化工程学院，2008.Che Yanqiu. Bifurcation in Neuron Models Exposed to External Electric Fields[D]. Tianjin：School of Electrical Engineering and Automation，Tianjin University，2008(in Chinese).

［12］ 李会艳. 神经系统的非线性动力学分析与控制[D]. 天津：天津大学管理与经济学部，2006.Li Huiyan. Analysis and Control of Nonlinear Dynamics of Neural Systems[D]. Tianjin：School of Management，Tianjin University，2006(in Chinese).

［13］ Prescott S A，de Koninck Y，Sejnowski T J. Biophysical basis for three distinct dynamical mechanisms of action potential initiation[J].PLoS Computational Biology，2008，4(10)：e1000198.

［14］ Izhikevich E M.Dynamical Systems in Neuroscience：The Geometry of Excitability and Bursting[M].London：The MIT Press，2005.

［15］ Radman T，Ramos R L，Brumberg J C，et al. Role of cortical cell type and morphology in sub- and suprathreshold uniform electric field stimulation[J].Brain Stimulation，2009，2(4)：215-228.

［16］ Bedard C，Kroger H，Destexhe A. Model of low-pass filtering of local field potentials in brain tissue[J].Phys Rev E，2006，73：051911.

［17］ 金启涛，王 江，伊国胜，等. 经颅磁刺激感应外电场作用下最小神经元模型放电起始动态机理分析[J].物理学报，2012，61(11)：118701-1-10.Jin Qitao，Wang Jiang，Yi Guosheng，et al. Action potential initial dynamical mechanism analysis in a minimum neuron model exposure to TMS induced electric field[J].Acta Phys Sin，2012，61(11)：118701-1-10(in Chinese).

［18］ Gai Y，Doiron B，Rinzel J. Slope based stochastic resonance：How noise enables phasic neurons to encode slow signals[J].PLoS Computational Biology，2010，6(6)：e1000825.