基于有限元-神经元模型的大鼠硬膜外脊髓电刺激仿真研究

2015-09-16陈清阳王熠钊

中国生物医学工程学报 2015年3期

关键词：背根电势脊髓

陈清阳徐琦* 王熠钊孟勇

1(华中科技大学自动化学院，图像信息处理与智能控制教育部重点实验室，武汉 430074)2(华中科技大学同济医学院附属同济医院康复科，武汉 430030)

基于有限元-神经元模型的大鼠硬膜外脊髓电刺激仿真研究

陈清阳1徐琦1*王熠钊2孟勇1

1(华中科技大学自动化学院，图像信息处理与智能控制教育部重点实验室，武汉 430074)2(华中科技大学同济医学院附属同济医院康复科，武汉 430030)

临床研究证实，硬膜外脊髓电刺激(ESCS)能促进脊髓损伤后的运动功能恢复，但基本神经机制则需通过动物实验进行研究。在建立大鼠ESCS有限元-神经元组合模型的基础上，仿真L2节段单阴极刺激条件下脊髓组织内的电势分布，确定L1到S2节段背根、腹根纤维和不同深度背柱纤维激活阈值的变化规律，分析纤维位置和刺激脉宽对神经纤维激活阈值的影响。仿真结果发现，距电极最近的背根纤维激活阈值最低为0.41 V，浅层背柱纤维的激活阈值略高为0.47 V，最近腹根纤维的激活阈值最高为0.78 V，减小电极-纤维距离有利于脊髓纤维的选择性激活；不同深度背柱纤维的激活阈值随刺激脉宽增加而减小，但脉宽过大导致背柱纤维激活阈值的降幅变小，过长的脉宽使得激活阈值的降幅趋缓，合理选择刺激脉宽有利于激活深层背柱纤维，增加脊髓组织的激活区域，提高ESCS对于深层背柱纤维的募集能力。该仿真结果可为动物实验研究中合理选择刺激参数、提高刺激选择性能提供理论指导。

有限元模型；神经元模型；硬膜外脊髓电刺激；大鼠；脊髓纤维

引言

硬膜外脊髓电刺激(epidural spinal cord stimulation，ESCS)将电极植入脊髓硬膜外腔，通过刺激电极将电刺激脉冲传递至脊髓内直接刺激神经组织。近年利用计算机仿真成为研究神经康复方法作用机理的有效途径之一，国内研究人员针对神经电刺激用于视觉功能修复和触觉功能重建的机理进行了仿真研究，为神经假肢的设计提供了理论指导[1-2]。临床实验证实，ESCS具有镇痛和促进脊髓损伤后运动功能恢复的疗效。Holsheimer等采用有限元方法，建立了人体ESCS模型，用于指导镇痛治疗中刺激参数的选择[3]。Rattay等仿真分析了脊髓神经纤维的激活顺序，研究电极位置与目标肌肉激活的关系[4]。

大鼠作为主要实验动物，被广泛用于医学研究，探索神经电刺激疗法的作用机制[5]，但如何设计合理的刺激方案则需通过计算机模型的仿真结果进行指导[6]，已有文献的ESCS仿真都是针对人体脊髓，而大鼠脊髓的解剖结构与人体脊髓相比，在几何尺寸、神经纤维直径范围以及脑脊液厚度等方面均存在明显差异，所以人体脊髓模型的仿真结果不能用于指导大鼠电刺激实验方案的制订。本研究在建立大鼠脊髓腰骶段ESCS有限元模型的基础上，仿真单阴极刺激条件下脊髓组织内的电场分布，确定脊髓神经纤维的胞外电势，然后利用神经元模型计算脊髓神经纤维激活阈值，分析不同刺激参数的作用效果，得到不同神经纤维的募集顺序，用于指导刺激电极的设计和合理选择刺激参数。

1 方法

1.1大鼠的ESCS有限元模型

采用COMSOL Multiphysics 4.3b(COMSOL Inc. Burlington MA, USA)，建立大鼠脊髓腰骶段的有限元模型。根据大鼠脊髓L2节段横截面的尼氏染色图片，得到脊髓灰质、白质、脑脊液、硬膜和椎骨的几何边界，确定大鼠脊髓结构的二维边界及几何尺寸[7]。由于椎骨电阻率较高，椎骨形状对电刺激的影响较小，为降低脊髓几何模型的复杂度，将椎骨内部和外部边界分别近似成圆形和方形，得到大鼠脊髓腰骶段的二维模型(如图1所示)，其中第二周围层用来作为接地层，第一周围层用于调整电极触点与接地层之间的阻抗[8]。为避免模型边界影响脊髓内的电势分布，将L1～ S2节段总长25 mm的大鼠脊髓模型沿轴向两端分别拉伸50 mm，得到总长为125 mm的大鼠脊髓三维模型[9]，由此该模型的尺寸为11 mm×11 mm×125 mm。

图1 大鼠脊髓腰骶段的二维几何模型(1-灰质；2-白质；3-脑脊液；4-硬膜；5-硬膜外腔；6-椎骨；7-第一周围层； 8-第二周围层)Fig.1 Two-dimension geometric model of the lumbosacral segment of rat spinal cord(1-Gray matter; 2-White matter; 3-Cerebrospinal fluid; 4-Dura mater; 5-Epidural space; 6-Vertebra; 7-First surrounding layer; 8-Second surrounding layer)

由于刺激频率较低(低于100 Hz)，因此组织电导率可视为常数，该脊髓模型中各组织电导率如表1所示[10]。刺激电极是采用聚酰亚胺对金触点进行绝缘封装的平板电极，圆形电极触点直径为1 mm，工作的单阴极触点置于L2节段中点，如图2所示。

表1 ESCS模型中的组织电导率Tab.1 Conductivity of compartments in the ESCS model

图2 大鼠脊髓电刺激的三维几何模型Fig.2 Three-dimension model for ESCS in a rat

采用自由剖分四面体方法进行网格剖分，由于电极及其周围区域电场变化较剧烈，该区域网格剖分较密集，而其他随电场变化减弱的区域使用较稀疏的网格剖分。为获得脊髓组织内的电势分布，也就是神经纤维的胞外电势大小，使用COMSOL Multiphysics中默认的迭代器求解Laplace公式[3]，即

·σV=0

(1)

式中，σ为组织电导率，V为组织内电势，是哈密顿算子。

考虑了电极-组织界面的压降，使用Robin边界条件，有

σV·n=g(Vmetal-V)

(2)

式中，n为边界对应的正交矢量，g为电极-组织界面的电导，Vmetal为电极触点的电势。

模型的顶部和底部采用Neumann边界条件，有

σV·n=0

(3)

即没有电流从该表面流出。

为验证有限拉伸长度和横截面尺寸的ESCS模型有效性，笔者仿真相同刺激电压作用下所建模型与横截面长宽或高度增加一倍后模型的电势分布，比较不同脊髓截面内电势分布的相关系数(correlation coefficient，CC)、放大因子(magnification factor，MAG)和均方根误差(root mean square error， RMSE)[11]，定量评价不同有限元模型的计算结果差异。

2 结果

在对所建立模型的有效性进行验证的基础上，就大鼠脊髓横截面的电势分布、不同纤维位置对背柱纤维激活阈值的影响、不同节段位置对背根和腹根纤维激活阈值的影响、脉宽对不同深度背柱纤维激活阈值的影响等方面进行仿真。

2.1大鼠ESCS有限元模型的有效性

当L2节段阴极触点的刺激电压为-5 V时，笔者分别采用脊髓有限元模型1(11 mm×11 mm×125 mm)、延伸长度增加一倍的模型2(11 mm×11 mm×225 mm)和横截面长宽各增加一倍的模型3(22 mm×22 mm×125 mm)，分别计算L2中点横截面内的电势分布，不同有限元模型的相关系数、放大因子和均方根误差如表3所示。由图3可知，有限元模型1的几何尺寸增加一倍后，电势分布的相关系数和放大因子均接近于1，均方根误差约为1 μV，远低于纤维激活阈值，表明该有限元模型边界能很好地近似脊髓电刺激的无限边界条件，可用于脊髓组织内电势分布的仿真计算。

表3不同有限元模型横截面内电势分布的相关系数、放大因子和均方根误差

Tab.3Correlationcoefficient,magnificationfactorandrootmeansquareerrorofthepotentialdistributioninthetransverseplanebetweenfiniteelementmodels

模型1和模型2模型1和模型3相关系数0.99860.9895放大因子1.08750.9987均方根误差/μV1.04381.0179

2.2大鼠脊髓横截面的电势分布

当刺激电压为-5 V时，大鼠脊髓L2节段中点横截面内的电势分布如图4所示。由图4可见，较高电势区域集中在脊髓背部浅层，即刺激阴极附近的电势幅值较高，变化也最为剧烈；而随着脊髓深度的增加，组织中的电势迅速减小，灰质内的电势分布接近于0 V，ESCS的刺激电场在脊髓浅层组织内扩散，较难作用于深层组织。

图4 单极刺激下L2节段中点横截面内电势分布Fig.4 The potential distribution of the transverse plane at the middle point of the spinal segment L2 in rat with the monopolar stimulation

2.3不同纤维位置对背柱纤维激活阈值的影响

在脉宽400 μs单极条件下，大鼠脊髓背柱不同位置纤维的激活阈值随电极-纤维距离的变化如图5所示。沿背腹方向，背柱纤维的激活阈值随电极-纤维距离的增大而迅速增加，最浅背柱纤维(电极-纤维距离为0.05 mm)的激活阈值为0.47 V，而最深背柱纤维(电极-纤维距离为1.25 mm)的激活阈值增加为8.61 V；沿中侧方向，浅层背柱纤维的激活阈值随电极-纤维距离的增大而增加缓慢，当电极-纤维距离从0.2 mm增加到0.8 mm时，激活阈值由0.53 V变为0.72 V。由于刺激电流在脊髓组织内扩散，浅层背柱纤维容易大面积激活，而深层背柱纤维不易被激活。

图5 不同位置背柱纤维的激活阈值随电极-纤维距离的变化Fig.5 Variation of the excitation threshold as a function of electrode-fiber distance for dorsal column fibers located in the dorsoventral and mediolateral direction

2.4不同节段位置对背根和腹根纤维激活阈值的影响

在脉宽400 μs单极条件下，L1到S2节段背根和腹根纤维激活阈值的变化如图6所示。随着脊髓根纤维进入脊髓位置与电极距离的增加，其激活阈值逐渐增大，距阴极最近的L2节段背根和腹根纤维激活阈值分别为0.41和0.78 V，而距阴极最远的S2节段背根和腹根纤维的激活阈值增大至6.51和7.83 V，相同脊髓节段的背根纤维激活阈值低于腹根纤维，可见背根纤维较腹根纤维更易被激活。

图6 从L1到S2节段脊髓根纤维激活阈值与脊髓节段的关系Fig.6 Variation of the excitation threshold as a function of segmental level for spinal root fibers from L1 to S2

2.5脉宽对不同深度背柱纤维激活阈值的影响

大鼠沿背腹方向，不同深度背柱纤维激活阈值与脉宽的关系如图7所示。随着脉宽增加，背柱纤维激活阈值逐渐变小，浅层背柱纤维激活阈值较低，而深层背柱纤维激活阈值较高。当脉宽为200 μs时，不同深度背柱纤维激活阈值由0.62 V增加到10.65 V；随着刺激脉宽增加，不同深度背柱纤维激活阈值的变化范围明显减小；当脉宽为1 600 μs时，不同深度背柱纤维激活阈值的变化范围为0.21 ～ 3.64 V。可见，采用长脉宽的刺激脉冲，有利于深层背柱纤维的激活。

图7 不同深度背柱纤维激活阈值与刺激脉宽之间的关系Fig.7 Variation of the excitation thresholds as a function of pulse width for dorsal column fibers located in dorsoventral direction

3 讨论与结论

脊髓电刺激所募集的神经元结构主要有背柱、背根和腹根纤维，对于临床使用的刺激电压范围0 ～ 10 V，背根纤维首先被激活，这是由于背根纤维靠近脊髓入口处，轨迹出现明显弯曲，脊髓入口处脑脊液与脊髓白质的电导率显著不同，使沿背根纤维的电势分布在脊髓入口边界处的变化最为剧烈，神经纤维激励函数取最大值，容易引发动作电位的产生[15]；和腹根纤维相比，背柱纤维更靠近刺激电极，具有较低的激活阈值。

随着刺激脉宽增加，背柱纤维激活阈值逐渐减小，有利于激活深层背柱纤维，而不会明显增加浅层组织的电流密度；但过长的脉宽会使激活阈值的降幅趋缓，导致被刺激脊髓组织的注入电荷量增加，更容易产生组织损伤[16]。因此，在综合考虑生物安全性的情况下，合理选择刺激脉宽对于深层背柱纤维的激活具有重要价值。

本研究基于有限元-神经元模型仿真，得到脊髓背柱、背根和腹根纤维激活阈值随纤维位置和脉宽的变化规律，这与人体仿真以及临床实验的结果是一致的，而大鼠实验中刺激脊髓腰骶段的运动阈为(0.3±0.1)V[6]，与本研究仿真得到电极正下方背根纤维的激活阈值0.41V基本一致。但是，本研究的有限元-神经元模型仍存在以下局限：一是由于目前尚无大鼠脊髓组织的电导率数据，因此所建立的有限元模型参考已有文献给出的人体脊髓组织电导率，并忽略了组织的容性效应，使仿真得到的神经纤维胞外电势存在一定误差，从而产生神经纤维激活阈值的计算误差，但不会对神经纤维激活阈值的变化规律和不同种类神经纤维募集顺序的仿真结果产生实质性影响；二是神经元模型中大鼠背柱、背根和腹根纤维的相关参数，主要来自对已有数据的外推插值而非实际测量得到，使得模型仿真与动物实验的结果不能完全一致。

阵列电极(array electrode)是近年由电极阵列发展而来的一种新的刺激技术，在改进刺激选择性和控制能力方面具有优越性能。笔者通过有限元-神经元模型仿真，发现脊髓入口位置离刺激阴极附近的脊髓根纤维容易激活，背根纤维的激活阈值低于腹根；对于背柱纤维，浅层纤维比深层纤维容易激活，合理的长脉宽有利于深层背柱纤维的募集。该结论对于探索脊髓电刺激的作用机理具有一定的意义，可用于指导动物实验方案的制订以及阵列电极的刺激参数选择，为增强刺激效果和提高靶向控制能力提供理论依据。

[1] 夏楠, 侯文生, 章毅, 等. 经硬脑膜电刺激对视皮层神经活动影响的仿真研究 [J]. 中国生物医学工程学报, 2011,30(5):725-730.

[2] 徐飞, 张定国, 许恒. 电触觉的皮肤神经机理仿真及实验验证 [J]. 中国生物医学工程学报, 2014,33(1):22-27.

[3] Holsheimer J, Wesselink WA. Optimum electrode geometry for spinal cord stimulation:the narrow bipole and tripole [J]. Medical and Biology and Engineer and Computing, 1997,35(5):493-497.

[4] Rattay F, Minassian K, Dimitrijevic MR. Epidural electrical stimulation of posterior structures of the human lumbosacral cord: 2. quantitative analysis by computer modeling [J]. Spinal Cord, 2000,38(9):473-489.

[5] Gerasimenko YP, Lavrov IA, Courtine G,etal. Spinal cord reflexes induced by epidural spinal cord stimulation in normal awake rats [J]. Journal of Neuroscience Methods, 2006,157(2): 253-263.

[6] Zhou Hui, Xu Qi, He Jiping,etal. A fully implanted programmable stimulator based on wireless communication for epidural spinal cord stimulation in rats [J]. Journal of Neuroscience Methods, 2012,204(2): 341-348.

[7] Senqul G, Watson C, Tanaka L,etal. Atlas of spinal cords of the rat, mouse, marmoset, macaque, and human [M]. San Diego: Elsevier Academic Press, 2012:15-72.

[8] Lee D, Hershey B, Bradley K,etal. Predicted effects of pulse width programming in spinal cord stimulation: a mathematical modeling study [J]. Medical and Biology and Engineering and Computing, 2011,49(7):765-774.

[9] 杨安峰, 王平. 大鼠的解剖和组织 [M]. 北京: 科学出版社, 1985: 159-161.

[10] McIntyre CC, Richardson AG, Grill WM. Modeling the excitability of mammalian nerve fibers:influence of afterpotentials on the recovery cycle [J]. Journal of Neurophysiology, 2002,87(2): 995-1006.

[11] Raspopovic S, Capogrosso M, Micera S. A computational model for the stimulation of rat sciatic nerve using a transverse intrafascicular multichannel electrode [J]. IEEE Transaction on Neural Systems and Rehabilitation Engineering, 2011,19(4): 333-344.

[12] Suh YS, Chung K, Goggeshall RE. A study of axonal diameters and areas in lumbosacral roots and nerves in the rat [J]. The Journal of Comparative Neurology, 1984,222(1): 473-481.

[13] Ladenbauer J, Minassian K, Hofstoetter US,etal. Stimulation of the human lumbar spinal cord with implanted and surface electrodes: a computer simulation study [J]. IEEE Transaction on Neural Systems and Rehabilitation Engineering, 2010,18(6): 637-645.

[14] Struijk JJ, Holsheimer J, Spincemaille GH,etal. Theoretical performance and clinical evaluation of transverse tripolar spinal cord stimulation [J]. IEEE Transactions on Rehabilitation Engineering, 1998,6(3): 277-285.

[15] Struijk JJ, Holsheimer J, Boom HBK. Excitation of dorsal root fibers in spinal cord stimulation: a theoretical study [J]. IEEE Transactions on Bimomedical Engineering, 1993,40(7): 632-639.

[16] Butterwick A, Vankov A, Huie P,etal. Tissue damage by pulsed electrical stimulation [J]. IEEE Transactions on Bimomedical Engineering, 2007,54(12): 2261-2267.

Simulation of Epidural Stimulation in Rat Spinal Cord Based on Finite Element and Neuron Model

Chen Qingyang1Xu Qi1*Wang Yizhao2Meng Yong1

1(Key Laboratory of Image Processing and Intelligent Control of Education Ministry, School of Automation, Huazhong University of Science and Technology, Wuhan 430074, China)2(Department of Rehabilitation Medicine, Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan 430030, China)

Clinical research demonstrated that epidural spinal cord stimulation (ESCS) was able to facilitate motor function recovery after spinal cord injury, but the underlying neural mechanism must be explored by animal experiments. In this paper, based on the construction of finite element and neuron model of ESCS, the distribution of electric potentials was simulated in the rat spinal cord with a stimulating cathode located at the midpoint of L2spinal segment. The excitation thresholds of typical spinal fibers were calculated to analyze the influence of the fiber location and pulse width on the recruitment of nerve fibers. The simulation results showed that the dorsal root fiber entering the spinal cord at the cathode level had the lowest excitation threshold of 0.41 V, whereas the excitation threshold of the ventral root fiber at the same level was 0.78 V and much higher than that of the most superficial dorsal column fiber with the excitation threshold of 0.47 V. The decrease of the fiber-electrode distance was able to improve the selectivity of fiber activation. It was also shown that the excitation threshold of dorsal column fibers at different depths within the midsagittal plane of spinal cord decreased with the increase of pulse width, whereas too long pulse width corresponded to the decreased slope of the varied activation threshold for dorsal column fibers. Therefore, the choice of an appropriate pulse width is essential to the activation of deep dorsal column fibers and the increasing activated area of the stimulated tissue. Our simulation results would offer guidance for the selection of stimulation parameters and the improvement of the stimulation selectivity during animal experiments.

finite element model; neuron model; epidural spinal cord stimulation; rat; spinal cord fiber

10.3969/j.issn.0258-8021. 2015. 03.007

2014-12-15，录用日期:2015-05-04

国家自然科学基金(60874035)；国家重点基础研究(973)项目(2013CB329500)

R318.38

0258-8021(2015) 03-0308-07

*通信作者(Corresponding author)， E-mail: xuqi@hust.edu.cn