夏 楠 侯文生* 章 毅 鄭小林 吳小鷹 陰正勤 JIANG Yingtao

1（重慶大學(xué)生物工程學(xué)院，重慶 400044）2（第三軍醫(yī)大學(xué)西南醫(yī)院眼科，重慶 400038）3（Department of Electrical and Computer Engineering，University of Nevada，Las Vegas，NV 89154，USA）

引言

經(jīng)硬腦膜對視皮層電刺激為視覺功能修復(fù)提供了一種新的技術(shù)思路［1］，即通過置于顱骨下、硬腦膜上的微電極陣列對視皮層進(jìn)行刺激，相對植入視皮層內(nèi)部的視覺假體［2］而言，它具有手術(shù)難度低、對視皮層組織損傷小等優(yōu)點。但經(jīng)硬腦膜的刺激電流需要流經(jīng)硬腦膜、腦脊液等組織才能到達(dá)目標(biāo)視皮層組織，認(rèn)識和掌握刺激電流對目標(biāo)皮層的作用效果非常必要。目前許多學(xué)者采用有限元法、邊界元法或有限體元法等數(shù)值方法來求解組織內(nèi)的電場分布［3］，嘗試用激活函數(shù)判定電場對神經(jīng)元電活動影響的可行性［4］；本課題組在之前工作中也通過建立貓視皮層的有限元模型，利用激活函數(shù)來表征電極的作用范圍，即主觀選取某個閾值來定義組織活化區(qū)，定性的討論經(jīng)硬腦膜對目標(biāo)皮層施加電刺激時電極的作用效果，初步證實了經(jīng)硬腦膜刺激的可行性［5］。

NEURON是一款由耶魯大學(xué)Carnevale和Hines開發(fā)的神經(jīng)計算軟件，為單個神經(jīng)元和神經(jīng)網(wǎng)絡(luò)的理論仿真提供一個方便快捷的研究平臺，CameronC.McIntyre和 WarrenM.Grill等人借助此平臺研究了深部腦刺激中高頻電刺激對電極周圍神經(jīng)元是否有興奮或抑制效應(yīng)，并很好地確定了該效應(yīng)的范圍［6］。由于激活函數(shù)判定方法的閾值選取是主觀的，可能會影響仿真工作的準(zhǔn)確性，所以本課題組在之前工作的基礎(chǔ)上，采用了一種新的方法，即將神經(jīng)元模型耦合到有限元模型中，進(jìn)一步確定經(jīng)硬腦膜刺激時，電極對視皮層的作用范圍。首先采用基于有限元法來分析不同刺激模式下腦組織內(nèi)的電場分布，然后將電位加載到在NEURON軟件中構(gòu)建的神經(jīng)元計算模型上，通過NEURON軟件計算不同位置神經(jīng)元是否發(fā)生動作電位來判定電極的作用范圍及不同刺激模式的作用效果。

1 方法

1.1 腦組織內(nèi)電場分布分析

選用COMSOLMultiphysics 3.5a（COMSOL Inc.，Burlington MA，USA）建立有限元仿真模型［5，7］，模型結(jié)構(gòu)見圖 1 中（a）和（c）。模型中采用3×3電極陣列進(jìn)行仿真研究，其結(jié)構(gòu)如圖1（b）所示。它是一種表面電極陣列，每個電極為等直徑的圓柱體，內(nèi)部設(shè)有通孔?；捉朴诰匦危笮∨c初級視皮層區(qū)域外表面積相當(dāng)，具體參數(shù)見表1。金屬導(dǎo)線位于基底內(nèi)，向外延伸至焊點處（未畫出）。由于是對單側(cè)大腦半球進(jìn)行刺激，故模型一側(cè)邊界取為大腦鐮與非刺激側(cè)大腦半球的交界面。模型總體尺寸為15 mm×10 mm×15.2 mm，具體尺寸參數(shù)信息見文獻(xiàn)［5，7］。

由于施加的是低頻刺激（小于10 kHz），此時可忽略生物組織的介電常數(shù)和電導(dǎo)率在不同空間位置的差異，將其視為各向同性均勻分布［5］，模型各部分電學(xué)參數(shù)見表2。對模型進(jìn)行邊界條件設(shè)定以及網(wǎng)格的劃分后對模型進(jìn)行求解，以得到經(jīng)硬腦膜刺激后的皮層內(nèi)電場分布。

圖1 經(jīng)硬腦膜對貓視皮層電刺激仿真模型。（a）三維仿真模型；（b）刺激電極陣列；（c）經(jīng)刺激電極軸線的二維截面Fig.1 Simulation modelofepiduralelectrical stimulation applied on the cat's visual cortex.（a）3D view of the model；（b）stimulating electrodes array；（c）2D cross section through the axis of one electrode

表1 刺激電極陣列模型參數(shù)Tab.1 Dimensions of electrodes array model

1.2 神經(jīng)計算軟件原理

NEURON軟件基本原理是通過求解線纜方程［8］，即式（1）得到神經(jīng)元在不同時刻的膜電壓。Iinj表示刺激電極流入胞內(nèi)的電流，cm和 rm分別表示膜電容和膜電阻，Vrest代表細(xì)胞的靜息電位，Vm為神經(jīng)元的跨膜電位。其基本方法就是根據(jù)神經(jīng)元生理結(jié)構(gòu)的差異將其分成若干節(jié)段，每個節(jié)段的生理結(jié)構(gòu)電學(xué)性質(zhì)完全相同，可等效為直徑不變的圓柱體，形成概念模型［7－8］，隨后根據(jù)計每個節(jié)段計算精度要求將其分成若干隔室，每個隔室用等效電路代替，形成計算模型［7－8］，并進(jìn)行計算，如圖 2所示。獲得理論結(jié)果后，與實際生理模型進(jìn)行比較，并對理論模型不斷進(jìn)行改進(jìn)，從而得到更精確的計算結(jié)果。其中隔室分的越多則計算精度就越高，但相應(yīng)的計算量越大。

表2 模型電學(xué)參數(shù)Tab.2 Electric parameters in simulation model

圖2 NEURON 計算過程［13］Fig.2 Computating process of NEURON

1.3 神經(jīng)計算方法

首先在NEURON中建立神經(jīng)元模型。為了使仿真結(jié)果為動物實驗提供理論指導(dǎo)，本研究以貓初級視皮層第III層內(nèi)的錐體細(xì)胞為研究對象，該模型由 Mainen根據(jù) Mainen和 Sejnowski的報道［9］所建立，在此基礎(chǔ)上，增加了軸突結(jié)構(gòu)如圖3（b）所示，除樹突、軸突、胞體外，還包括軸丘、起始段等完整的神經(jīng)元解剖結(jié)構(gòu)。其次，在有限元模型（圖1）體坐標(biāo)中建立一個11×11的二維正方形點陣，相鄰兩點間距0.2 mm。該點陣位于中心電極正下方，平行于電極基底，貓初級視皮層第III層內(nèi)的錐體細(xì)胞對應(yīng)位置，且矩陣中心點對應(yīng)中心電極正下方神經(jīng)元，矩陣邊緣包含了兩個近端電極下方皮層的神經(jīng)元活動信息，如圖3（c）所示。

點陣中每一點代表一個神經(jīng)元，若某一點位置上的神經(jīng)元爆發(fā)動作電位，就可認(rèn)為以該神經(jīng)元為中心的0.01 mm2區(qū)域內(nèi)的神經(jīng)元均發(fā)生興奮。由于建立的腦模型尺寸遠(yuǎn)小于貓腦直徑（約5 cm），為簡化問題，可認(rèn)為在該區(qū)域中的神經(jīng)元都為此種神經(jīng)元。最后通過NEURON中的play函數(shù)將傅里葉-有限元方法計算出的電位加載到相應(yīng)位置的神經(jīng)元模型上，并計算該神經(jīng)元是否發(fā)生興奮，從而得到電極的作用范圍。

圖3 神經(jīng)元模型與皮層關(guān)系示意圖。（a）有限元模型中的特征截面；（b）貓視皮層第III內(nèi)錐體細(xì)胞神經(jīng)元模型；（c）神經(jīng)元點陣分別在截面Ⅰ、Ⅱ、Ⅲ中的位置Fig.3 Relationship between neuron model and cortex.（a）characteristic crosssectionsin Fourier finite element model；（b）model of the pyramidal cell in layer III in cat's visual cortex；（c）the position of neuron dot array in cross sectionⅠ，Ⅱ，Ⅲ

1.4 刺激參數(shù)的選擇

實驗選取雙向脈沖、指數(shù)脈沖和正弦波3種刺激形式，來研究經(jīng)硬腦膜電刺激時，刺激形式的改變對視皮層內(nèi)神經(jīng)元作用范圍的影響。3種刺激形式的波形和主要參數(shù)見表3。3種刺激波形由Matlab自帶數(shù)學(xué)函數(shù)產(chǎn)生，通過設(shè)置波形的脈寬、幅值和頻率來改變刺激模式。

2 結(jié)果

2.1 不同波形刺激時，皮層內(nèi)神經(jīng)元的響應(yīng)范圍

實驗選取3種波形的基準(zhǔn)刺激幅值為0.8mA，刺激頻率為200Hz，其中雙向方波脈沖的正向基準(zhǔn)脈寬為0.2 ms。皮層內(nèi)電場分布及神經(jīng)元點陣中的響應(yīng)區(qū)域結(jié)果見圖4。

表3 3種刺激波形參數(shù)Tab.3 Parameters of three stimulus waveforms

圖4 3種波形刺激結(jié)果（黑色點表示未響應(yīng)的神經(jīng)元區(qū)域，白色點表示爆發(fā)動作電位的神經(jīng)元區(qū)域，箭頭指向大腦鐮方向）。（a）截面Ⅲ中的神經(jīng)元點陣；（b）雙脈沖刺激；（c）指數(shù)脈沖刺激；（d）正弦波刺激。Fig.4 The results of three stimulus waveforms（black dot represents the unexciting neuron，the white dot represents the excitatory neuron，the arrow indicates the cerebral falx direction.）.（a）neuron dot array in cross section Ⅲ；（b）bipulse stimulus；（c）exponential pulse stimulus；（d）sinusoidal pulse stimulus.

從圖4可看出，3種波形都引起發(fā)了神經(jīng)細(xì)胞的興奮，不同的刺激波形引發(fā)的神經(jīng)元區(qū)域興奮數(shù)量不等。在所有的刺激方式下，靠近大腦鐮一側(cè)的神經(jīng)元區(qū)域興奮的數(shù)量都多于遠(yuǎn)離大腦鐮方向的神經(jīng)元區(qū)域數(shù)量，且響應(yīng)平面區(qū)域呈一個以刺激電極點為中心的非對稱不規(guī)則圓形。

神經(jīng)元的興奮與其膜外電位有關(guān)，本研究在有限元模型中第Ⅰ平面上神經(jīng)元放置區(qū)域做了一條特征截線L，如圖5（a）所示，箭頭指向大腦鐮方向，截線L的電位曲線如圖5（b）所示。圖5（b）中電極中心位置為零，可以看出在對稱位置上電壓不等，近大腦鐮一側(cè)電壓略高于遠(yuǎn)離大腦鐮一側(cè)的電壓，在下面的結(jié)果中都存在這種現(xiàn)象。

圖5 特征截線L的電位分布。（a）特征曲線L的位置；（b）曲線L上的電位Fig.5 Potential distribution in characteristic line L.（a）position of line L；（b）electric potential on line L.

2.2 刺激幅值的變化對皮層內(nèi)神經(jīng)元響應(yīng)范圍的影響

刺激幅值變化也會引起皮層內(nèi)電場分布改變，因此也會引起皮層內(nèi)的神經(jīng)元區(qū)域興奮數(shù)量變化，圖6（a）顯示了在不同幅值作用下神經(jīng)元興奮數(shù)量的變化情況。三種波形的刺激頻率為200Hz，其中雙向方波脈沖的正向脈寬為0.2 ms。

從刺激結(jié)果來看神經(jīng)元響應(yīng)的區(qū)域都隨著刺激幅值的增加而增大，用指數(shù)波進(jìn)行刺激時，幅值減小到0.5 mA皮層依然有9個神經(jīng)元區(qū)域爆發(fā)動作電位，隨著刺激幅值的增大，矩陣內(nèi)神經(jīng)元區(qū)域響應(yīng)個數(shù)快速增多，幅值增加到0.8 mA時近大腦鐮一側(cè)神經(jīng)元響應(yīng)區(qū)域到達(dá)點陣的邊緣（如圖4（c）），0.9 mA時響應(yīng)區(qū)域已經(jīng)超出了點陣。雙向脈沖刺激時，隨著刺激幅值的增大矩陣內(nèi)神經(jīng)元區(qū)域響應(yīng)個數(shù)隨之增加，但只有刺激幅值高于0.7 mA時神經(jīng)元才有響應(yīng)。正弦波在同等幅值作用下響應(yīng)的神經(jīng)元范圍要遠(yuǎn)遠(yuǎn)大于前兩種刺激，幅值增大到0.9 mA時矩陣內(nèi)神經(jīng)元區(qū)域全部發(fā)生興奮。

2.3 脈寬以及頻率改變對皮層內(nèi)神經(jīng)元的影響

脈寬也是方波脈沖的重要調(diào)節(jié)參數(shù)，分別用幅值為0.8 mA、正向脈寬為0.2、0.3和0.4 ms的雙向脈沖對模型進(jìn)行刺激，刺激頻率為200Hz，所得到的神經(jīng)元區(qū)域響應(yīng)結(jié)果見6（b）。用幅值為0.8 mA、頻率分別為100、200和300Hz的指數(shù)脈沖、正弦波和雙脈沖刺激目標(biāo)皮層，其中雙向方波脈沖脈寬為0.2 ms，結(jié)果見圖6（c）。

圖6 改變刺激參數(shù)對神經(jīng)元區(qū)域興奮數(shù)量的影響。（a）改變刺激幅值；（b）雙脈沖刺激下脈寬變化；（c）改變刺激頻率Fig.6 Varies of the number of excitatory neuron regions with different stimulus.（a）different stimulation amplitude；（b）different duration under bipulse stimulus；（b）different frequency of three stimulus waveforms

從圖6（b）可以看出在刺激幅值和頻率不變時，增加正向方波的脈寬皮層內(nèi)響應(yīng)的神經(jīng)元區(qū)域也隨著增大，但是增加或減小雙向方波脈沖的頻率對刺激范圍沒有影響，如圖6（c）所示，說明用雙向方波脈沖刺激時，在正向方波占空比較小的情況下，視皮層內(nèi)神經(jīng)元對刺激頻率變化不敏感，而對脈寬的變化非常敏感，脈寬由0.2 ms增加到0.3 ms時，響應(yīng)的神經(jīng)元區(qū)域個數(shù)從7個增加到47個，由0.3 ms增加到0.4 ms時響應(yīng)區(qū)域個數(shù)增加到76。在幅值不變條件下，隨著刺激頻率的增加指數(shù)脈沖和正弦波引發(fā)的神經(jīng)元區(qū)域響應(yīng)個數(shù)都有減小的趨勢，當(dāng)頻率增大時，刺激周期減小，正向刺激電荷量相應(yīng)減少，從而引起皮層內(nèi)神經(jīng)元區(qū)域響應(yīng)數(shù)量減少。因此在考慮非方波刺激形式對視皮層內(nèi)神經(jīng)元響應(yīng)區(qū)域數(shù)量影響時，頻率也是比較重要的因素。

3 討論與總結(jié)

由于錐體細(xì)胞占到皮層內(nèi)細(xì)胞總量的80％，且為易興奮細(xì)胞［10］，因此本研究只選用貓腦初級視皮層內(nèi)第III層內(nèi)的錐體細(xì)胞作為分析對象，將其排列成間距為0.2 mm的二維點陣耦合到有限元模型中，間距設(shè)置越小，點陣中的神經(jīng)元個數(shù)就越多，所描述的響應(yīng)范圍越精確，但計算量也隨之增大，綜合考慮間距定為0.02 mm較合適。與之前的仿真方法［5］相比，本研究將選定的神經(jīng)元耦合到傅里葉-有限元模型中（圖3），將計算出的電位直接加載到各神經(jīng)元上即可得到電極的作用范圍，減少了主觀選取閾值的不確定性，可以清楚、直觀地描述出神經(jīng)元響應(yīng)區(qū)域的細(xì)微差別（如圖4），因此該方法有一定的優(yōu)越性。

研究選取雙脈沖、指數(shù)脈沖和正弦波3種刺激波形，從實驗結(jié)果看，3種波形都可引起皮層內(nèi)神經(jīng)元的興奮，且興奮區(qū)域都呈一個非對稱的近圓形（圖4），近大腦鐮一側(cè)的神經(jīng)元區(qū)域興奮的數(shù)量多于遠(yuǎn)大腦鐮側(cè)的興奮數(shù)量。近大腦鐮一側(cè)皮層被腦脊液包圍（圖1），而腦脊液電導(dǎo)率遠(yuǎn)大于其他組織電導(dǎo)率［7］，且其電導(dǎo)率的細(xì)微變化都會較大程度地影響模型電位分布［7］，所以該解剖結(jié)構(gòu)的差異使得近大腦鐮側(cè)的神經(jīng)元膜外電位略高于遠(yuǎn)大腦鐮側(cè)電位（圖5），從而引起了興奮區(qū)域的不對稱性。

在減小3種波形刺激幅值時，神經(jīng)元區(qū)域興奮數(shù)量也隨之減少，這與文獻(xiàn)［5］的結(jié)論相符。但本研究仿真中也得出了與文獻(xiàn)［5］不一致的結(jié)果：一方面，在一定范圍內(nèi)增大或減小雙向方波的刺激頻率并沒有使刺激范圍發(fā)生變化（如圖6（c）所示），對于雙向脈沖而言，在單個周期內(nèi)正向脈沖刺激完成后，組織內(nèi)剩余電荷立刻被反向電流抵消，從而消除了每個周期組織內(nèi)電荷的累積［12］，而改變刺激的重復(fù)頻率只改變了刺激的時間間隔，并沒有改變輸入的電流量，因此不能引起更多神經(jīng)元的興奮；另外，在一定范圍內(nèi)減小指數(shù)脈沖重復(fù)頻率時，電極的刺激范圍有增大趨勢，對指數(shù)波而言，減小重復(fù)頻率相當(dāng)于增加了單個周期內(nèi)流入皮層內(nèi)的電流量［11］，從而使響應(yīng)范圍增大。

本研究只關(guān)注電極在不同刺激模式下所產(chǎn)生的電場對皮層內(nèi)神經(jīng)元的作用范圍，因此是將單個神經(jīng)元視為獨立個體，形成一個平面點陣來討論的，并沒有構(gòu)建神經(jīng)元網(wǎng)絡(luò)來討論神經(jīng)元之間的影響，且只針對視皮層內(nèi)某一平面的神經(jīng)元活動做了定量分析，僅討論了平行于電極的某個二維平面的神經(jīng)元響應(yīng)范圍，沒有考慮刺激深度的問題，因此對于兩種軟件的結(jié)合使用還需要更深一步的優(yōu)化，從而得到更準(zhǔn)確的結(jié)果。

［1］郭冰冰，侯文生，陰正勤，等.一種經(jīng)硬腦膜的多通道視覺誘發(fā)電位實驗測試系統(tǒng)的設(shè)計［J］.儀器儀表學(xué)報，2009，30（6，增刊）：204－207.

［2］牛金海，劉易非，任秋實，等.視覺假體：帶給盲人光明的電子設(shè)備［J］.中國科學(xué)（E輯：信息科學(xué)），2007，37（10）：1354－1362.

［3］顏威利，徐桂芝，著.生物醫(yī)學(xué)電磁場數(shù)值分析［M］.北京：機(jī)械工業(yè)出版社，2006.

［4］Wongsarnpigoon A，GrillWM.Computationalmodelingof epidural cortical stimulation［J］.J Neural Eng，2008，5：446 －454.

［5］侯文生，章毅，鄭小林，等.基于有限元法的硬腦膜外視皮層電刺激仿真研究［J］.中國生物醫(yī)學(xué)工程學(xué)報，2010，29（4）：557－563.

［6］McIntyre CC，Grill WM，Sherman DL，et al.Cellular effects of deep brain stimulation：model-based analysis of activation and inhibition［J］.J Neurophysiol，2004，91：1457 － 1469.

［7］張毅.經(jīng)硬腦膜對視皮層電刺激的數(shù)字仿真［D］.重慶：重慶大學(xué)，2010.

［8］Dominique D.The somatic shunt cable model for neurons［J］.Biophysical Journal 1984，46：645－653.

［9］Mainen ZF，Sejnowski TJ.Influence of dendritic structure on firing pattern in model neocortical neurons［J］.Nature，1996，382：363－366.

［10］壽天德.視覺信息處理的腦機(jī)制［M］.上海：上?？萍冀逃霭嫔?，1999.

［11］Sahin M，Tie Yanmei.Non-rectangular waveforms for neural stimulation with practical electrodes［J］.J Neural Eng，2007，4：227－223.

［12］Lilly JC，Hughes JR.Brief，noninjurious electric waveform for stimulation of the brain［J］.Science，1955，121：468 －469.

［13］Carnevale NT，Hines ML.The NEURON Book［M］.New York：Cambridge University Press，2006.