基于有限元法的硬脑膜外视皮层电刺激仿真研究

侯文生 章 毅 郑小林 吴小鹰 阴正勤 JIANG Yingtao BISWAJIT Das(重庆大学生物工程学院，重庆 400044)

2(第三军医大学西南医院眼科，重庆 400038)

3(Department of Electrical and Computer Engineering，University of Nevada，Las Vegas，NV 89154，USA)

引言

视觉是人类感知世界的主要途径，失明将带来极大的不便。视觉假体通过对视网膜、视神经、视皮层的电刺激重建视觉，近年来得到了广泛关注［1］。目前的视皮层视觉假体是在皮层表面［2］或皮层内［3］施加刺激，这破坏了皮层组织的完整性，同时存在脑脊液溢出的风险。一种新的思路是通过硬脑膜对视皮层进行刺激［4］，相对而言，它具有手术难度低、对组织伤害小等优点;但此时刺激电极与目标皮层间存在硬脑膜、脑脊液等组织，它们将改变皮层内的电场分布，从而对刺激方案的选择提出了新的要求。为从理论上认识经硬脑膜对视皮层的电刺激效应，首先分析硬脑膜外电刺激下脑组织内的电场分布，并采用电位对空间向量的二阶偏导数来表征此电场对神经元电活动的影响［5］。通过改变刺激参数，优化刺激方案，为实验提供理论指导。

其中σ 和ε 分别表示组织的电导率和介电常数，ω表示谐波的角频率，得到的V 乘以谐波的相位Φ 并取实部后即为组织内的电位，介电常数ε 的引入是考虑了生物组织的电容效应［7］。各组织的电导率和介电常数均视为各向同性均匀分布［8］，如表2所示［5，8］。

图1 经硬脑膜对猫视皮层电刺激仿真模型。(a)刺激电极阵列;(b)三维仿真模型;(c)经刺激电极轴线的二维截面。Fig. 1 The simulation model of epidural electrical stimulation applied on the cat’s visual cortex. (a)the stimulating electrodes array;(b)3D view of the model;(c)2D cross section through the axis of one electrode.

1 方法

1.1 脑组织内电场分布分析

采用基于有限元法的 COMSOL Multiphysics 3.5a(COMSOL Inc.，Burlington MA，USA)来分析脑组织内的电场分布。首先在考虑猫初级视皮层区域实际解剖结构的前提下［6］，参照文献［5］中的运动皮层区域硬脑膜外电刺激仿真模型，建立了一个由头皮、颅骨、硬脑膜、脑脊液、灰质(皮层)、白质6 层组织构成的模型，包括LS(Lateral Sulcus)和SS(Splenial Sulcus)这两条主要的脑沟。仿真模型中的电极阵列与动物实验中的电极阵列相似［4］，包括一个块状聚酰亚胺基底和四个圆柱状金电极，电极引脚通过金属导线与外接焊点(未画出)一一相连(图1(a))。电极位于硬脑膜上对应于初级视皮层区域的部位，中心点距大脑纵裂约2 mm［6］，上面覆压颅骨［5］。由于是对单侧大脑半球进行刺激，故模型一侧边界取为大脑镰与非刺激侧大脑半球的交界面。模型总体尺寸为15 mm ×10 mm ×15.2 mm(图1(b))，它是从2D 模型(图1(c))拉伸而得的，主要尺寸参数见表1［5～6］。

由于任何周期性刺激信号均可展开成一系列正弦或余弦信号(谐波)，通过分析谐波作用下的电场分布，可研究周期性刺激信号的作用效果，此时模型内各点电位可通过求解复数形式的Laplace 方程而得，

表1 模型尺寸参数Tab.1 Dimensions of the simulation model

采用四面体单元进行网格划分，刺激电极附近网格最密，并逐渐变稀疏;但为充分反映皮层内电场分布的变化，增加了电极正下方皮层内网格密度。为保证电位对空间向量二阶偏导数的连续性，单元阶次取为3，一共划分了206627 个单元。不同组织交界处电流密度连续;模型中除上表面(头皮表面)外的所有外边界均视为无穷远边界(V =0)，而头皮表面设置为绝缘边界(n·J =0)［5］。采用单电极或双电极进行单极性电流刺激，后者分为两种，即经距离较远的电极对(5 mm)或经距离较近的电极对(2 mm)进行刺激(称为远距离双电极刺激和近距离双电极刺激)，三种刺激模式下经每个电极进入模型的电流均为I。增大模型的网格密度或几何尺寸后，在单电极1 mA 直流刺激下，模型内最大电位的变化率均不超过2%，表明模型的网格密度和相关尺寸参数是有效的。

表2 模型电学参数Tab.2 Electrical properties of the simulation model

1.2 外加电场对神经元电活动影响的分析

电位V 对任意空间向量u 的二阶偏导数也称为激活函数(Activating Function，AF)，能表征外加电场对与此向量平行分布的神经元电活动的影响，激活函数为正时，神经元将发生去极化，从而有可能产生兴奋;反之则发生超极化并抑制兴奋［5，9］。考虑到视皮层内的主要兴奋性神经元(锥体细胞)垂直皮层表面分布［10］，且主要关心刺激电极正下方皮层内神经元的电活动，只研究了沿z 轴方向的激活函数，此时u 平行于z 轴(图1)。因不同神经元兴奋时所需激活函数大小不一，故设定了一个阈值，将激活函数大于该阈值的点连接起来，得到的区域可称为组织活化区(Volume of Tissue Activated，VTA)［6］，认为刺激作用下该区域内的神经元均会产生兴奋。由于视皮层假体采用的是电极阵列进行多通道刺激，每个电极的组织活化区越小，相邻电极间相互影响就越小，从而提高了刺激的空间分辨率。

为方便观察仿真结果，经刺激电极轴线沿x、y轴选取了两个截面，并在距皮层表面0.15 mm 处选取了第三个截面(因距皮层表面约0.15 mm 的第Ι层内没有神经元分布［10］)，分别称为截面Ⅰ、Ⅱ、Ⅲ(图3(a))。通过研究这三个截面上激活函数及组织活化区的变化情况，可以直观了解外加电场对神经元电活动的影响，并定性分析刺激的作用效果。

图2 单相脉冲及其频谱。(a)脉冲波形;(b)复数频谱。Fig. 2 Waveform and spectrum of the monophasic pulse. (a)waveform;(b)complex spectrum.

1.3 刺激参数的选择

单相脉冲f(t)是传统的神经电刺激信号［11］，本研究亦主要探讨此种形式的刺激。其参数包括脉冲幅值A、脉冲宽度τ 及周期T (或频率f )，一个周期内的波形如图2(a)所示。将f(t)展开成指数形式的傅立叶级数，第n 个分量为:

其中n = -∞，…，-1，1，…，+∞

其复数频谱如图2(b)所示，各次谐波频率相差f，随着频率的增大，其幅值将以1/n 的速率衰减，而能量将以1/n2的速率衰减。

由于基波分量包含了原始刺激信号的主要能量及频率信息，为简化计算，将其作为原始刺激信号的代表，逐次改变各项参数，通过研究组织活化区的变化，定性了解不同刺激参数的作用效果。对单相脉冲而言，其基波分量的幅值为相位为0。在探讨不同刺激参数的作用效果之前，首先研究直流刺激下脑组织内电场分布及其对神经元电活动的影响，以求对硬脑膜外电刺激的机制有比较明确的认识。

2 结果

2.1 直流刺激下脑组织内电场分布及其对神经元电活动的影响

施加1 mA 直流刺激后，模型内电流密度分布如图3 所示。电流在流经脑脊液后将发生横向扩散，使单个电极的作用范围扩大。在单电极刺激(图3 中(b)～(d))下，由于组织内电导率各向同性均匀分布，同一组织内电流密度等值面近似球冠形分布;对双电极刺激而言，电流密度则以垂直于刺激电极对连线的平面对称分布(图3 中(e)和(f));离连线中心点越近，电流密度越小。

图3 1 mA 直流刺激下电流密度(mA/mm2)分布(箭头表示电流方向)。(a)截面Ⅰ、Ⅱ、Ⅲ示意图(不同刺激模式下的电流密度取自然对数，并绘制于相应截面上);(b)～(d)单电极刺激时截面Ⅰ～截面Ⅲ;(e)远距离双电极刺激时截面Ⅰ;(f)近距离双电极刺激时截面ⅡFig.3 Distribution of the current density (mA/mm2)generated by 1 mA DC current (Arrows indicating the direction of current). (a)schematic drawing of cross section Ⅰ，Ⅱ and Ⅲ. (The natural logarithm of the magnitude of current density under different stimulation mode were plotted on corresponding cross sections，including);(b)～(d)cross section Ⅰ～cross sectionⅢunder single electrode stimulation;(e)cross sectionⅠ under double electrodes stimulation with long distance;(f)cross section Ⅱunder double electrodes stimulation with short distance.

1 mA 直流刺激下模型内相关电学物理量如表3 所示。可见在多电极刺激模式下，当电极间距足够大时(超过5 mm)，脑组织内电学量分布与单电极刺激相近，即只作用于其周围局部组织;减小刺激电极间距，相邻电极的相互影响增加。此外，电流从刺激电极传导至皮层后，电流密度下降了约180 倍，电位下降了30 倍左右，表明要达到与皮层刺激相似的效果，需要施加更大的刺激强度。值得注意的是，模型内最大电流密度仅仅集中在刺激电极与硬脑膜接触面的边缘部分，随后将迅速衰减，因而刺激强度稍大时，组织并不会受到很大损害。

表3 直流刺激下脑组织内电场分布Tab.3 Distribution of electric field under DC stimulation

为便于进一步分析，将三个截面的交点称为基点(Pb)，截面Ⅰ与Ⅱ的交线称为基线(Lb)，并用LⅠ、LⅡ表示Ⅰ与Ⅲ、Ⅱ与Ⅲ的交线，统称为特征截面和特征截线(图4(a))。1 mA 直流刺激下，皮层及白质内部分区域的激活函数如图4 所示，其最大值位于基线与皮层表面的交点处，并以此点为中心衰减。对单电极刺激(图4 中(b)～(d))而言，激活函数等值面近似球冠状分布;双电极刺激下，仅刺激电极正下方局部区域内保持球冠形分布，在其它区域将发生一定改变;特别是位于两电极连线中心平面附近的区域，由于大部分电流叠加后抵消，导致该处电场强度(电位的一阶导数)很小，进而导致激活函数值几乎接近0 分布(见图4(f))。需要注意的是，基点处激活函数值在单电极刺激下为787.439 mV/mm2，远距离双电极刺激时几乎不变(785.413 mV/mm2)，仅近距离双电极刺激时略有下降(779.547 mV/mm2)。这意味着即使在多电极刺激模式下，各电极下方的神经元电活动仍然由此电极上施加的刺激决定;当电极间距较大(超过5 mm)时，单个电极的作用范围几乎不受其它电极的影响。

理论上可指定任意正数为激活函数阈值来确定神经元兴奋与否［9］，但取值越小，电极的作用范围越大，从而在多电极刺激模式下，相邻电极的组织活化区将发生重叠(图5(a))。为避免此种情况的出现，选取110 mV/mm2作为1 mA 直流刺激下神经元兴奋时所需的激活函数阈值，可使近距离双电极刺激下，两个电极的组织活化区刚好接触而不重叠(图5(b))，其等值面将与截面III 形成一个封闭区域(即组织活化区)。单电极刺激的组织活化区近似于以基点为底面圆心、基线为高线、垂直指向白质表面的球缺(图5(c)，(d))。而双电极刺激下，当电极距离较近时，组织活化区将发生一定形变，靠近皮层表面的区域产生缺失;距离较远时两个电极各自的组织活化区仍然近似球缺形分布(图5(e)，(f))。在组织活化区内对常数1 进行积分，可得单电极刺激下的组织活化区体积为1.053 mm3，远距离双电极刺激下的组织活化区体积为1.040 mm3，近距离双电极刺激下的组织活化区体积为1.207 mm3。

图4 1 mA 直流刺激下激活函数(mV/mm2)分布(纵坐标表示特征截线上各点的激活函数值，横坐标表示相应点到基点的距离)。(a)特征截面和特征截线(激活函数绘制于相应截面和截线上);(b)单电极刺激时截面I 和截线LⅠ内;(c)单电极刺激时截面II 和截线Lb 内;(d)单电极刺激时截面III 和截线LⅡ内;(e)近距离双电极刺激时截面II 和截线LⅡ内;(f )远距离双电极刺激时截面III 和截线LⅠ内Fig. 4 Distribution of activating function (mV/mm2 )generated by 1 mA DC current (The vertical axis in the vignette represents the activating function value，while the lateral axis represents the distance to the base point). (a)characteristic cross sections and characteristic lines (The activating function was plot on correlative cross sections and lines;(b) cross section Ⅰ and line LⅠ under single electrode stimulation;(c)cross section Ⅱand line Lb under single electrode stimulation;(d)cross section Ⅲand line LⅡunder single electrode stimulation;(e)cross section Ⅱand line LⅡunder double electrodes stimulation with short distance;(f)cross section Ⅲ and line LⅠ under double electrodes stimulation with long distance.

图5 不同激活函数阈值下各刺激模式的组织活化区分布。(a)阈值过小导致组织活化区重叠;(b)选择合适阈值使近距离双电极刺激组织活化区接触而不重叠(此阈值下不同刺激模式的组织活化区投影于相应截面上);(c)单电极刺激时在截面Ⅲ的投影;(d)单电极刺激时在截面Ⅱ的投影;(e)远距离双电极刺激时在截面Ⅲ的投影;(f)远距离双电极刺激时在截面Ⅰ的投影Fig.5 Distribution of VTA with different threshold of activating function for each stimulation mode. (a)overlapping of VTA caused by small threshold;(b)proper threshold was chosen to make the VTA of double electrodes with short distance just contacting. And with this threshold (VTA within different stimulation mode was projected on correlative cross sections);(c)the projection of VTA on cross section Ⅲ under single electrode stimulation;(d)the projection of VTA on cross section Ⅱunder single electrode stimulation;(e)the projection of VTA on cross section Ⅲunder double electrodes stimulation with long distance; (f) the projection of VTA on cross section Ⅰ under double electrodes stimulation with long distance.

2.2 刺激参数对刺激效果的影响

在研究参数变化对刺激效果的影响时，对单电极刺激而言，仍然是固定某个阈值，研究组织活化区体积随参数变化的改变情况。对双电极刺激而言，选取一个使各自的组织活化区刚好接触而不重叠的激活函数阈值，研究该阈值随刺激参数变化的改变情况;需要注意的是神经元兴奋时所需激活函数值是固定的，两电极组织活化区刚好接触时的阈值越大，其相互影响越严重。

根据前期动物实验所用参数，选择A =1 mA，τ= 0.2 ms，f = 200 Hz (记为A0，τ0，f0)的单相脉冲为基准刺激模式，其基波信号幅值为0.08 mA。此时远距离双电极刺激下激活函数阈值为3.6 mV/mm2;近距离双电极刺激的激活函数阈值为8.8 mV/mm2，这也是单电极刺激所需的激活函数阈值，单电极刺激的组织活化区体积仍然为1.053 mm3。分别改变A、τ、f，组织活化区体积VVTA及激活函数阈值AFthrl(远距离双电极刺激)或AFthrs(近距离双电极刺激)的相对变化如图6 所示。可以看出，无论是增大刺激幅值、刺激脉宽还是刺激频率，单电极刺激的组织活化区体积及双电极刺激的激活函数阈值均增大，这意味着将有更多神经元被激活，其原因在于这三个参数的增加都会导致刺激信号单个周期能量的增大，从而引起激活函数值增加和单个电极的组织活化区扩大。此外，幅值改变对组织活化区的影响大于脉宽和频率对组织活化区的影响，这是由于在基波分量幅值表达式中刺激幅值A对其影响较大;频率改变对组织活化区造成的影响略小于脉宽改变的影响，这可能与生物组织的电容效应有关［6］。

图6 刺激参数改变时组织活化区相对变化率。(a)单电极刺激下组织活化区体积;(b)远距离双电极刺激下激活函数阈值;(c)近距离双电极刺激下激活函数阈值Fig.6 The relative variation ratio of VTA with changes of stimulus parameters. (a)volume of VTA under single electrode stimulation;(b)threshold of activating function for double electrodes stimulation with long distance;(c)threshold of activating function for double electrodes stimulation with short distance

3 讨论与结论

与传统的生物组织电场分析［5］相比，研究中考虑了生物组织的电容效应;建立的仿真模型包括了主要的解剖结构。虽然需要通过MRI 切片重建模及考虑各向异性的电学特性，才能得到与实际电场分布很相近的仿真结果［11］，但由于仅需研究初级视皮层附近局部区域内的电场分布，所建立的模型尺寸也远小于猫脑直径(约5 cm)［7］，故模型内不同组织层厚度均匀，其电学参数也各向同性均匀分布，所得仿真结果是可以接受的。需要注意的是，本研究没有考虑电极与周围组织接触处的双电层电容和法拉第阻抗［12］，认为两者电位相等。这是因为在小电流刺激下，可以不考虑法拉第阻抗［12］;又由于采用的是恒流刺激，可进一步忽略双电层电容［6］，仅脑组织的电阻及电容特性会对其内的电场分布及刺激效果产生影响。

通过刺激函数来研究外加刺激对神经元电活动的影响，虽然神经元的兴奋性不仅与激活函数的大小相关，还取决于神经元自身的尺寸、外形等特性，但在短时脉冲作用下，激活函数仍是评估外加电场对神经元电活动影响的有效指标［5］。此外，虽然引起视觉的是视皮层上相关神经集群的电活动［10］，但研究单个神经元的电活动是研究集群电活动的基础。因而将激活函数作为定性了解外加刺激作用效果的指标是合理的。

主要研究了单相脉冲各项参数对刺激效果的影响，这是由于此种刺激模式对神经活动具有很好的选择性［12］，但长时间使用会对组织造成伤害。一种新的刺激方式是双相(biphasic)脉冲，此时单个周期内电流沿不同方向流过组织，从而消除了电荷的累积效应，降低了对组织和电极的损害［12］;下降型三角脉冲、指数脉冲及Guass 脉冲等非方波刺激形式由于对神经活动具有较好的选择性也得到了一定应用，但其作用效果还需更多深入研究［13］。

本研究通过计算机仿真，发现经硬脑膜可以实现对视皮层的有效刺激，只是需要比直接皮层刺激约大30 倍的刺激强度。在多电极刺激下，只要电极间距超过5 mm，单个电极的作用效果就几乎不受其它电极的影响，降低幅值、脉宽和频率等参数，可以改善单相方波脉冲的刺激效果，其中幅值的减小对空间分辨率的提高效果最为显著。所得结论为实验提供了一定理论指导。

［1 ］ 任秋实. 视觉假体的研究进展与面临的挑战［J］. 生命科学，2009，21(2):234 -240.

［2 ］ Chowdhury V，Morler JW，Coroneo MT. Surface stimulation of the brain with a prototype array for a visual cotrx prothesis［J］.Clinical Neuroscience，2004，11(7):750 -755.

［3 ］ Warren DJ，Fernandez E，Normann RA. High-resolution twodimensional spatial mapping of cat striate cortex using a 100-microelectrode array ［J］. Neuroscience，2001，105 (1):19-31.

［4 ］ 郭冰冰，侯文生，阴正勤，等. 一种经硬脑膜的多通道视觉诱发电位实验测试系统的设计［J］. 仪器仪表学报，2009，30(6，增刊):204 -207.

［5 ］ Wongsarnpigoon A， Grill WM. Computational modeling of epidural cortical stimulation［J］. J Neural Eng，2008，5:446 -454.

［6 ］ Butson CR，McIntyre CC. Tissue and electrode capacitance reduce neural activation volumes during deep brain stimulation［J］. Clinical Neurophysiology，2005，116:2490 -2500.

［7 ］ Kalia M，Whitteridge D. Visual areas in the splenial sulcus of the cat［J］. J Physiol，1973，232:275 -283.

［8 ］ Gabriely S，Lau RW，Gabriel C. The dielectric properties of biological tissues:Ⅱ. Measurements in the frequency range 10 Hz to 20 GHz［J］. Phys Med Biol，1996，41:2251 -2269.

［9 ］ Rattay F. The basic mechanism for the electrical stimulation of the nervous system ［J］. Neuroscience，1999，89(2):335-346.

［10］ 寿天德. 视觉信息处理的脑机制［M］. 上海:上海科技教育出版社，1999.2.

［11］ Gilad O，Horesh L，Holder DS. Design of electrodes and current limits for low frequency electrical impedance tomography of the brain［J］. Med Bio Eng Comput，2007，45:621 -633.

［12］ Merrill DR，Bikson M，Jefferys J. Electrical stimulation of excitable tissue:design of efficacious and safe protocols ［J］.Journal of Neuroscience Methods，2005，141:171 -198.

［13］ Sahin M，Tie Yanmei. Non-rectangular waveforms for neural stimulation with practical electrodes［J］. J Neural Eng，2007，4:227 -233.