基于能量概率與微孔力模型的脈沖電場對細胞電穿孔動態(tài)過程的仿真分析

姚陳果 呂彥鵬 趙亞軍 董守龍 劉紅梅

(輸配電裝備及系統(tǒng)安全與新技術(shù)國家重點實驗室(重慶大學) 重慶 400030)

基于能量概率與微孔力模型的脈沖電場對細胞電穿孔動態(tài)過程的仿真分析

姚陳果 呂彥鵬 趙亞軍 董守龍 劉紅梅

(輸配電裝備及系統(tǒng)安全與新技術(shù)國家重點實驗室(重慶大學) 重慶 400030)

采用Matlab軟件基于跨膜電壓、微孔半徑以及孔的受力隨時間和空間變化的分析來模擬細胞在1kV/cm、100μs脈沖期間電穿孔動態(tài)變化過程。采用的模型相比傳統(tǒng)模型有兩點明顯的改進：第一，利用微孔能量的概率模型判斷微孔半徑變化的初始時間，當脈沖電場注入的能量大于微孔的最大能量的概率超過0.99時，微孔半徑開始變化；第二，建立微孔受力模型來描述微孔的發(fā)展和分布，提出新的參數(shù)-孔沖量來決定微孔的大小。結(jié)果表明，除了極點不發(fā)生穿孔外，細胞膜都發(fā)生電穿孔，細胞膜大部分區(qū)域的微孔的受力接近于0，孔徑基本區(qū)域穩(wěn)定，不發(fā)生變化；但靠近極點附近的小部分區(qū)域孔徑受力為正值，如果延長脈沖時間，孔徑仍然在發(fā)展。此模型可以用來判斷微孔的形成與預測微孔的發(fā)展。

電穿孔 跨膜電壓 能量概率 微孔力

0 引言

脈沖功率技術(shù)已經(jīng)在生物電磁領(lǐng)域中得到了廣泛的應用[1-3]。當細胞處于外加脈沖電場時，細胞膜上會發(fā)生電穿孔。電穿孔是一種脈沖作用下細胞膜上脂質(zhì)分子發(fā)生重排后產(chǎn)生微孔的物理現(xiàn)象。這些微孔一般只存在幾秒鐘或幾分鐘，微孔的出現(xiàn)會極大的增加細胞膜的電導率和通透性[4],為生物活性分子和離子提供短暫的傳輸通道使其進入細胞內(nèi)部。目前電穿孔可以應用于基因傳輸[5，6]、電化學治療[7，8]及腫瘤治療[9，10]。電穿孔現(xiàn)象也在生物電磁領(lǐng)域中受到了諸多學者的青睞，也正在逐漸被醫(yī)學領(lǐng)域接受。

電穿孔的出現(xiàn)與跨膜電壓(TransmembranePotential，TMP)存在一定的關(guān)系，跨膜電壓超過一定的閾值(200～1 000mV)時，細胞膜被擊穿，膜上會出現(xiàn)短暫的微孔。脈沖撤去后，如果微孔會慢慢恢復，并且在無外界刺激的前提下，微孔恢復后不會對細胞的生理活動造成影響，這種現(xiàn)象叫做可逆電穿孔。這種可逆電穿孔效應可以使DNA等生物活性分子進入細胞中從而影響細胞的功能[11-13]。如果脈沖參數(shù)繼續(xù)增強，產(chǎn)生的微孔的尺寸很大且脈沖后微孔不能自行恢復，則將會對細胞膜產(chǎn)生不可逆的破壞，會造成細胞的死亡。這種現(xiàn)象叫做不可逆電穿孔[14,15]。由于不可逆電穿孔對微型器官和組織的非熱殺傷效應，重慶大學姚陳果等首次將不可逆電穿孔應用到治療腫瘤中，發(fā)現(xiàn)腫瘤能夠完全被消融[16]。R.Neal等在腎組織的消融實驗中取得了成功[17]，Y.Mandel等也利用不可逆電穿孔技術(shù)成功的消融了眼睛摘除后的原發(fā)性葡萄膜黑色素瘤[18]。

近幾年來，為了解電穿孔的動態(tài)發(fā)展過程，眾多學者建立了模型與仿真。重慶大學米彥等通過建立細胞等效電路來模擬跨膜電壓的頻率響應，用于反應頻率對電穿孔的影響[19]。姚陳果等通過建立分子動力學仿真來模擬電穿孔過程[20]。杜克大學W.Krassowska等利用Matlab建立單個細胞的模型來模擬跨膜電壓與孔密度的變化，用于反應電穿孔的時域發(fā)展過程[21]。

上述模型從不同方面描述了細胞中的電穿孔動態(tài)過程。電穿孔不僅與外加電場、空間和時間相關(guān)，還與細胞的半徑、類型以及形狀相關(guān)。上述分子動力學建模時，分析了電穿孔的形成及其變化過程，但沒有定量地分析孔徑的變化。電路模型是對整個細胞為對象進行建立的，并沒有顯示出細胞膜上的跨膜電壓等參數(shù)隨細胞膜空間位置改變的變化過程。上述杜克大學的研究[21]沒有考慮電穿孔中孔徑變化的影響，而在后續(xù)的研究[22]中，他們考慮了微孔半徑隨細胞膜空間位置及其時間的變化。但W.Krassowska認為只要在細胞膜某區(qū)域內(nèi)出現(xiàn)了第一個微孔，微孔半徑就開始發(fā)生變化。本文認為微孔的形成與能量相關(guān)，微孔想要擴大，脈沖電場注入的能量就必須大于微孔自身的能量。本文所建模型中將詳細介紹微孔能量與微孔能夠發(fā)生變化的關(guān)系。文獻[22]中，在特定脈沖參數(shù)下，微孔在脈沖期間形成并擴大。一旦脈沖撤去后，微孔全部恢復到原始狀態(tài)，因此W.Krassowska等模擬的是可逆電穿孔的動態(tài)發(fā)展過程。該模型描述了微孔參數(shù)的變化過程，但沒有詳細介紹微孔變化的原因。本文認為，微孔的發(fā)展與其受力相關(guān)，因此本文通過建立微孔受力模型來定量的描述微孔的變化過程，并通過定義新的參數(shù)——孔沖量來決定微孔半徑的大小。本文所建模型對細胞實驗具有一定的參考意義。

1 解析模型

1.1 數(shù)值模型

本文利用Matlab軟件建立半徑為10μm的單個細胞模型。細胞處于導電介質(zhì)中，脈沖參數(shù)選為1kV/cm、100μs[23]。細胞外液電勢φe和細胞內(nèi)電勢φi服從拉普拉斯方程

(1)

(2)

細胞模型如圖1所示。

圖1 球形細胞二維模型，電場方向從右向左Fig.1 Model of a spherical single cell，The electric field E is oriented from right to left

1.2 跨膜電勢

當細胞發(fā)生電穿孔時，微孔可以為電流提供傳輸通道。根據(jù)電流密度連續(xù)性原理，細胞外液與細胞膜交界面上的電流密度和細胞膜上的電流密度相等，因此可以根據(jù)電流連續(xù)性原理建立跨膜電壓隨時間和空間的變化[24，25]。

neJe=nmJm(t)=niJi

(3)

Je=λeEe

(4)

Vm=φe-φi

(5)

(6)

式中，Je、ne分別為細胞外液與細胞膜交界面的電流密度和電場法向；Jm、nm分別為細胞膜的電流密度和電場法向；Ji、ni分別為細胞膜與細胞內(nèi)液交界面的電流密度和電場法向；Ee為細胞外液與細胞膜交界面的電場法向分量；φe為細胞外液的電勢；φi為細胞內(nèi)液的電勢。

式(6)等號右邊，第1項為離子電流密度，第2項為電容電流密度，第3項為流過微孔的電流密度，其他符號定義和參數(shù)數(shù)值見表1。

表1 模型參數(shù)

Tab.1Modelparameters

符號定義數(shù)值R細胞半徑10μmh細胞膜半徑5nmλe細胞外液電導率0.14S/mλi細胞質(zhì)電導率0.3S/mλm0細胞膜初始電導率5×10-7S/mεe細胞外液介電常數(shù)80εi細胞質(zhì)介電常數(shù)80εm細胞膜介電常數(shù)5ω0孔的能量屏障系數(shù)2.65F法拉第常數(shù)9.65×104C·mol-1R氣體常數(shù)8314K-1mol-1T絕對溫度295Kq微孔生成率2.46N0初始孔密度1.5×109m-2α微孔生成系數(shù)1×109m-2·s-1N孔的相對長度0.15[21]r*親水孔的最小孔徑0.51×10-9mrmVm=0時的最小孔能量0.8×10-9mD孔半徑的擴散系數(shù)5×10-14mC空間排斥能系數(shù)9.65×10-15Jγ脂質(zhì)-孔能量系數(shù)1.8×10-11J·m-1σ0無孔時的脂質(zhì)張力系數(shù)1×10-3J·m-2σ'脂質(zhì)-水交界的張力系數(shù)2×10-2J·m-2A脂質(zhì)總面積1.26×10-9m-2FmaxVm=1V時的最大電場力0.70×10-9Nrh式(20)中的常數(shù)0.97×10-9mrt式(20)中的常數(shù)0.31×10-9mrw單個水分子的半徑0.967×10-10mη懸浮液的擴散速率0.9×10-3Pν單個水分子的摩爾體積18×10-6m3/molNA阿伏伽德羅常數(shù)6.022×1023k玻爾茲曼常數(shù)1.38×10-23

微孔的電流密度公式[26]為

JEP(t)=iEPN(t)

(7)

式中，iEP為單個微孔的電流；N(t)為微孔密度。在不區(qū)分通過微孔的離子前提下，根據(jù)能斯特-普朗克方程的模型，iEP可以表示為[26]

(8)

式中

(9)

(10)

(11)

(12)

在電穿孔過程中，細胞膜電導率隨時間的變化可表示為[27]

(13)

1.3 孔密度

電穿孔過程中理論上微孔可以親水孔的形式存在于細胞膜上來傳導電流的?？缒る妷簳S著外加電場變化，間接影響微孔密度的變化，孔密度的變化公式為[28，29]

(14)

式(14)中，初始孔密度可以設(shè)置為1.5×109m-2。式(14)中孔密度不僅能夠表示微孔密度的增加，也可以表示微孔密度的減小。在脈沖期間，孔密度會增加多個數(shù)量級；但脈沖撤去后，孔密度則以極小的速度減小。

1.4 孔徑的發(fā)展

未施加電場前，細胞膜上所有的微孔都是疏水孔，如圖2[30]所示。孔徑的大小由其屏障能量控制。當細胞處于高脈沖場強時，如果r>r*，大部分疏水孔將會自發(fā)的轉(zhuǎn)為親水孔，親水孔能夠迅速發(fā)展到最小能量對應的孔徑rm。因此本文假定初始的孔徑為rm。微孔的能量公式為[30]

(15)

式中，第1項為電場提供的能量;第2項為空間的排斥力提供的能量;第3項為線張力提供的能量;第4項為表面張力提供的能量。Ed表示為施加電場前微孔所能存儲的最大屏障能量，rd表示Ed所對應的孔徑。當外加電場作用于細胞時，Ed逐漸減小，且對應的孔徑也在減小。如果Ed對應的孔徑下降到rm，則Ed接近于0，表示孔徑不受屏障能量的束縛，孔徑會隨時間的變化而變化。但如果電場提供的能量較小時，微孔處于穩(wěn)定狀態(tài)。本文可以通過式(16)～式(18)來判斷孔徑是否發(fā)生變化[31]。

(16)

(17)

(18)

P(Vm，t)為微孔由疏水孔轉(zhuǎn)為親水孔且孔徑開始發(fā)展的概率，此時產(chǎn)生的隨機孔可認為是親水孔。如果P(Vm，t)>0.99，則微孔會擴大；否則微孔會處于一個穩(wěn)定的狀態(tài)。微孔擴大的公式可表示為[32]

(19)

圖2 Vm=0時，疏水孔與親水孔孔徑與微孔能量曲線圖Fig.2 Under the transmembrane potential Vm=0,energy of hydrophobic (rr*) pores as a function of pore radius

1.5 微孔受力模型

微孔的發(fā)展受到力的約束，如圖3[33]所示，當細胞處于強電場下，對于單個孔而言，在細胞膜上會出現(xiàn)電場力、線張力、表面張力和空間排斥力。

F(rj，Vm，σeff)為單個微孔的受力，即微孔力，其表達式為[32]

圖3 Vm=0時，單個微孔的綜合受力Fig 3 Comprehensive force on a single pore at Vm=0

2πγ+2πσeffrjr≥rm

(20)

式中，第1項為電場力;第2項為空間排斥力;第3項為線張力;第4項為表面張力。有效張力系數(shù)σeff的表達式為[32]

(21)

式中，A為細胞的表面積；Ap為細胞膜的上所有孔的面積，其表達式為

(22)

當細胞處于脈沖場強中，脈沖開始時，微孔受力平衡，表明微孔處于穩(wěn)定狀態(tài)。隨著時間的發(fā)展微孔受力變?yōu)檎５捎赑(Vm，t)<0.99，所以微孔不會變化，但當P(Vm，t)>0.99，且微孔受力仍然為正，微孔開始變化。如果脈沖撤去后，微孔受力仍然為正，表示微孔繼續(xù)擴大；如果為負，表示微孔會縮小到原始狀態(tài)。

1.6 仿真計算

用Matlab軟件對式(1)～式(4)、式(6)、式(13)～式(16)、式(19)、式(20)進行數(shù)值計算，來模擬電穿孔動態(tài)過程。細胞處于導電介質(zhì)中，細胞膜被劃分為200個子區(qū)域，θ的范圍在0°～360°，則每個區(qū)域內(nèi)Δθ=1.8°。仿真計算步驟如下：

1)設(shè)置跨膜電壓初始值為0.001V(接近于0)，孔密度為1.5×109m-2，電導率為5×10-7S·m-1，初始孔徑為0.8nm。

2)利用式(1)～式(4)、式(6)、式(13)、式(14)來更新跨膜電壓、孔密度、電導率。

3)用式(16)來更新P(Vm，t)，如果P(Vm，t)≥0.99，微孔半徑開始變化，孔徑變化的可用式(19)來更新；否則微孔的狀態(tài)不發(fā)生變化。

4)計算微孔每個區(qū)域的受力。脈沖期間，如果子區(qū)域內(nèi)微孔受力趨近于0，孔徑不會變化，表示此區(qū)域達到平衡，即電穿孔過程結(jié)束；否則回到步驟2)用已更新的跨膜電壓、孔密度和孔徑作為初始值繼續(xù)迭代直至脈沖結(jié)束。

2 數(shù)值結(jié)果

基于以前對細胞實驗和仿真的學習[22]，本文選用1kV/cm、100μs脈沖參數(shù)作為激勵源。

2.1 跨膜電壓

球形細胞處于外電場中會產(chǎn)生超極化和去極化。最大和最小的跨膜電壓都出現(xiàn)在細胞的赤道上。最大的跨膜電壓出現(xiàn)在超極化區(qū)域(0°～90°，270°～360°)的赤道上，最小的出現(xiàn)在去極化區(qū)域(90°～270°)的赤道上。細胞膜可以視為電容和電阻的組合。在超極化區(qū)域內(nèi)，如圖4a中的θ=0°，30°，80°，90°，假定細胞膜的電容不發(fā)生變化[21]。由于細胞膜上發(fā)生電穿孔，微孔會為離子傳輸提供通道，因此細胞膜電導率會大大增加?？缒る妷弘S著時間推移先增大后減小到穩(wěn)定值。在去極化區(qū)域內(nèi)，跨膜電壓會先減小后增加到穩(wěn)定值。設(shè)置細胞膜的靜息電位為0，如果θ1+θ2=180n，超極化區(qū)域θ1與去極化區(qū)域θ2具有對稱性。以超極化區(qū)域θ=0°為例，脈沖開始時，由于外電場的作用，細胞膜開始充電，因此跨膜電壓迅速上升到峰值1.34V。電穿孔產(chǎn)生的閾值電壓范圍在200～1 000mV，因此電穿孔發(fā)生在跨膜電壓上升期間。但由于電穿孔的產(chǎn)生，膜電導率迅速增加導致跨膜電壓下降并趨于穩(wěn)定值0.4V。每個區(qū)域的充電時間與放電時間也不同，但是跨膜電壓的變化過程是相似的。隨著時間的發(fā)展，以θ=0°為例，跨膜電壓迅速上升到峰值后緩慢下降到穩(wěn)定值，θ=30°處，充電時間和放電時間都要比θ=0°處的時間長，此區(qū)域峰值為1.32V，穩(wěn)定值為0.39V；θ=80°處，峰值為1.16V，穩(wěn)定值為0.003 7V。從赤道到極點的區(qū)域，峰值和穩(wěn)定值在不斷減小。在極點上，由于法向電場強度為0，因此跨膜電壓不發(fā)生變化。由于對稱性，上述分析也可以用到去極化區(qū)域內(nèi)，但不同的是其峰值和穩(wěn)定值為負，如圖4a所示。

隨著微孔的發(fā)展，跨膜電壓在同一時間不同區(qū)域的變化是不同的。脈沖開始時，所有區(qū)域的跨膜電壓都為0。當時間為0.83μs，跨膜電壓呈余弦變化。此時最大電壓出現(xiàn)在超極化區(qū)域的赤道上，最小電壓出現(xiàn)在去極化區(qū)域的赤道上。當跨膜電位超過一定幅值時(200～1 000mV)，細胞膜上出現(xiàn)微孔。如1.5μs，膜電導率迅速增加，赤道首先發(fā)生電壓跌落。到10μs時，電壓跌落的區(qū)域擴大，而且跨膜電壓的最大值在持續(xù)減小且出現(xiàn)的位置向極點轉(zhuǎn)移。尤其是在1.5～10μs之間表現(xiàn)的非常明顯。跨膜電壓在不同的位置其穩(wěn)定值是不同的。如圖4b所示，不同的脈沖時刻下，跨膜電壓的峰值及其峰值出現(xiàn)的位置是不同的。電導率隨著電穿孔的產(chǎn)生是不斷增大的，但由于每個區(qū)域的跨膜電壓變化不同，因此其電導率也不同，則導致每個區(qū)域的時間常數(shù)τ=C/G也不同。電穿孔首先發(fā)生在赤道，因此電導率首先增加，時間常數(shù)開始減小。離赤道越遠，電穿孔發(fā)生的時間越晚，時間常數(shù)就越大，即充電時間越長，因此跨膜電壓峰值出現(xiàn)的時間越晚。但在100μs時，靠近極點附近的區(qū)域內(nèi)，微孔半徑達到了μm級別，因此電導率非常大，并且跨膜電壓跌落非常嚴重，基本接近于0。

圖4 跨膜電壓的時間與空間的分布Fig 4 Time evolution and spatial distribution of thetransmembrane voltage Vm

2.2 微孔半徑與受力

微孔的變化主要是受到作用在微孔的力來控制的。以圖5a中的θ=0°為例，在脈沖起始階段，微孔受力由負值變?yōu)檎担悄芰扛怕蔖(Vm，t)<0.99，所以孔徑不發(fā)生變化；一旦概率P(Vm，t)>0.99且微孔受力為正，微孔開始發(fā)生變化。微孔受力與跨膜電壓和微孔密度相關(guān)。在充電期間，跨膜電壓迅速增大，電場力也迅速的增大，因此微孔受力迅速達到峰值，導致微孔半徑持續(xù)擴大。由于電穿孔的出現(xiàn)，跨膜電壓隨后會迅速下降，因此微孔受力會趨向于0，所以孔徑緩慢的增長并趨于穩(wěn)定。由于孔力與孔徑與跨膜電壓的平方相關(guān)，因此θ=0°和θ=180°，θ=30°和θ=150°，θ=80°和θ=100°的曲線分別重合。在不同的子區(qū)域中，由于微孔的受力有所不同，因此微孔半徑的發(fā)展過程也是不同的。如圖5b所示，θ=0°處孔徑為7.189nm，θ=30°時孔徑為8.050nm。在圖5a中，需要注意的是：不同位置上受力曲線與時間軸所包圍的面積是不同的；且θ=30°處所包圍的面積大于θ=0°時，所以半徑也比θ=0°時大。通過此仿真，可以定義一個新的變量來決定微孔的尺寸，表達式為

(23)

圖5 微孔受力與半徑隨時間的變化Fig.5 Time evolution of pore force and pore radius

此變量可稱為孔沖量，即微孔受力在微孔上對時間的積分?？讻_量隨時間的變化曲線如圖6所示?；诳讻_量的定義，如果孔沖量具有極限值，即隨時間的變化趨于穩(wěn)定。則孔徑也趨于穩(wěn)定，如果極限值越大則孔徑也越大。θ=30°處的孔沖量比θ=0°處大，因此θ=30°處的孔徑較大。但在θ=80°處，微孔受力并不趨于0，而是正值并持續(xù)增大，因此孔沖量并不收斂，所以孔徑繼續(xù)增大。由于對稱性，所以在超極化上的分析與去極化的分析是一樣的。

圖6 孔沖量M隨時間的變化Fig.6 Time evolution of pore momentum M

圖7 孔徑的空間分布Fig.7 Spatial distribution of pore radius

圖8 微孔受力的空間分布Fig.8 Spatial distribution of pore force

孔徑在不同的區(qū)域內(nèi)具有不同的分布，如圖7和圖8所示。t=2μs時，孔徑從赤道到極點逐漸減小，此時最大值為8.427nm。微孔受力在所有的位置上都是正值，所以孔徑在整個細胞膜上將會不斷的增加。隨時間的變化，孔徑最大值出現(xiàn)的位置由赤道向極點轉(zhuǎn)移。在t=50μs時刻，微孔受力在赤道上是負值，因此微孔開始收縮。在t=100μs時刻，微孔受力為負值的區(qū)域擴大，表明微孔收縮的區(qū)域也在擴大。雖然微孔受力為負，但是極小。在50μs時刻，θ=0°處的微孔受力為-3.200 5×10-13N，半徑為7.911 7nm；在100μs時刻，微孔受力為-1.688 9×10-12N，半徑為7.190nm，孔徑變化很小，因此可視微孔受力為負值的區(qū)域處于平衡狀態(tài)。但是由于不同子區(qū)域的孔沖量不同，即使區(qū)域穩(wěn)定，每個區(qū)域孔徑的穩(wěn)定值是不同的。脈沖結(jié)束后，大部分區(qū)域都穩(wěn)定，最小的孔徑為7.189 2nm。但是在極點附近，仍有一小部分區(qū)域的受力為正值，因此孔徑仍然在發(fā)展。最大的孔徑為1.461 5μm。如果繼續(xù)加大仿真時間，那么此區(qū)域的孔徑仍然繼續(xù)增大。

3 討論

本文研究了電脈沖場作用于單細胞時電穿孔的動態(tài)發(fā)展過程。仿真結(jié)果較好地驗證了模型的正確性，同時也可通過與國、內(nèi)外相關(guān)文獻對比，證明本文模型和計算結(jié)果的正確性。

文獻[34]中，脈沖期間測得細胞赤道的穩(wěn)定跨膜電壓為±0.5V，與本文中測得細胞赤道跨膜電壓±0.4V相對比，誤差在允許范圍之內(nèi)。另一方面，整個細胞膜跨膜電壓的空間變化與本文的仿真相吻合(如圖4a所示)，即脈沖起始階段赤道的跨膜電壓最大。隨著脈沖時間變化，跨膜電壓的最大值出現(xiàn)在穿孔區(qū)域和未穿孔區(qū)域的邊界上，與本文結(jié)果中跨膜電壓的變化具有一致性。

文獻[35]測出細胞膜上的部分微孔半徑約為20～120nm，可以使DNA分子自由的進入到細胞內(nèi)部(DNA分子直徑在5～9nm)。本文結(jié)果中脈沖結(jié)束后產(chǎn)生的最小微孔半徑為7.189 2nm，在極點附近產(chǎn)生的微孔半徑超過了100nm，證明本文模型中細胞膜上產(chǎn)生的微孔也能夠使DNA分子輕易地進入細胞內(nèi)部。本文結(jié)果中產(chǎn)生的孔徑也與ChangandReese等人的實驗結(jié)果[35]較為吻合。除此之外，國外相關(guān)研究學者利用不同直徑的熒光染料分子作為探針證明了脈沖作用下在細胞膜上能夠產(chǎn)生不同尺寸的微孔[36，37]，實驗表明小分子染料可大量進入細胞[36]，熒光強度較高；半徑較大的染料分子僅可少量進入細胞，熒光強度較弱[37]；而半徑過大的染料分子則不能進入細胞[36]，細胞內(nèi)部無熒光。因此證明了本文仿真結(jié)果中脈沖作用下能夠在細胞膜上產(chǎn)生不同尺寸的微孔的正確性。

文獻[16]中，脈沖處理過后的細胞在掃描電子顯微鏡下觀察發(fā)現(xiàn)了微米級別的微孔。文獻[38]用人工膜代替細胞膜進行脈沖處理后，也在人工膜上發(fā)現(xiàn)了微米級別的微孔。而在本文仿真中，極點附近出現(xiàn)了微米級的微孔，與實驗能夠很好的吻合。

W.Krassowska等[39]在理論分析中提出孔密度大的地方微孔半徑反而較小，即兩者呈反比例關(guān)系。如圖9所示，W.Krassowska的理念與本文模型結(jié)果相吻合，驗證了本文模型計算結(jié)果的正確性。

圖9 脈沖后(100 μs)微孔密度與微孔半徑的關(guān)系Fig.9 Relationship between pore densityand pore radius after the pulse (100 μs)

上述文獻證明了本文模型與結(jié)果的正確性，但脈沖作用于細胞的機理仍需進一步探究。下一步工作可利用微流體裝置設(shè)計完整的實驗方案從而更進一步探討脈沖對細胞的作用機理。

4 結(jié)論

本文基于跨膜電壓、孔密度、孔徑以及微孔受力等變量隨時間和空間的變化的分析，并通過加入能量概率模型與微孔受力模型來模擬電穿孔的動態(tài)發(fā)展過程。得到如下結(jié)論：

1)脈沖期間，除細胞膜極點位置上，其他區(qū)域都發(fā)生電穿孔。

2)微孔受力趨近于0表明微孔處于穩(wěn)定狀態(tài)，但不同區(qū)域的穩(wěn)定值不同?？讖椒€(wěn)定值的大小可由孔沖量決定。但靠近極點的區(qū)域微孔受力仍然為正值，表明此區(qū)域仍處于不平衡狀態(tài)；此外，孔沖量也不能夠收斂，表明孔徑仍處于發(fā)展狀態(tài)。

3)增加仿真時間，如果微孔受力為仍在為正，那么孔徑將會繼續(xù)增大。

本文模型旨在模擬細胞膜上電穿孔的形成與發(fā)展過程，而電穿孔的發(fā)展可能會導致細胞膜的破裂繼而導致細胞死亡。本文模型可以對細胞實驗具有一定的參數(shù)指導作用，尤其是針對殺死腫瘤細胞的實驗，不同的脈沖參數(shù)所產(chǎn)生的微孔大小及其分布都是不同的，因此可以根據(jù)本文模型的仿真分析來選擇細胞實驗脈沖參數(shù)。

[1]SersaG，CemazarM，MiklavcicD，etal.Electrochemotherapyoftumors[J].RadiologyandOncology，2006，40(3)：163-174.

[2] 熊蘭，馬龍，胡國輝，等.具有負載普適性的高壓雙極性方波脈沖源研制[J].電工技術(shù)學報，2015，30(12)：51-60.XiongLan，MaLong，HuGuohui，etal.Anewlyhigh-voltagesquarebipolarpulsegeneratorforvariousloads[J].TransactionsofChinaElectrotechnicalSociety,2015，30(12)：51-60.

[3] 米彥，張晏源，儲貽道，等.基于非平衡Blumlein型多層微帶傳輸線的高壓納秒脈沖發(fā)生器[J].電工技術(shù)學報，2015，30(11)：100-109.MiYan，ZhangYanyuan，ChuYidao，etal.High-voltagenanosecondpulsegeneratorbasedonnon-balancedblumleintypemultilayeredmicrostriptransmissionline[J].TransactionsofChinaElectrotechnicalSociety,2015，30(11)：100-109.

[4]AbidorIG，ArakelyanVB，ChernomordikLV，etal.ElectricalbreakdownofbilayerlipidmembranesI.Themainexperimentalfactsandtheirqualitativediscussion[J].BioelectrochemandBioenergetics，1979，6(1)：37-52.

[5]MatthewsKE，DevSB，ToneguzzoF,etal.Electroporationforgenetherapy[J].InMethodsinMolecularBiology：AnimalCellElectroporationandElectrofusionProtocols，1995，48(5)：273-280.

[6]NishiT，YoshizatoK，YamashiroS，etal.High-efficiencyinvivogenetransferusingintraarterialplasmidDNAinjectionfollowinginvivoelectroporation[J].CancerResearch，1996，56(5)：1050-1055.

[7]SukharevS，KlenchinV，SerovS，etal.ElectroporationandelectrophoreticDNAtransferintocells.TheeffectofDNAinteractionwithelectropores[J].BiophysicalJournal，1992，63(5)：1320-1327.

[8]PuciharG，MirL，MiklavcˇicD.Theeffectofpulserepetitionfrequencyontheuptakeintoelectropermeabilizedcellsinvitrowithpossibleapplicationsinelectrochemotherapy[J].Bioelectrochemistry，2002，57(2)：167-172.

[9] 郭飛，李成祥，唐賢倫，等.沖激輻射天線實現(xiàn)皮秒脈沖電場在人體大腦模型中聚焦的研究[J].電工技術(shù)學報，2016，31(3)：195-202.GuoFei，LiChengxiang，TangXianlun，etal.Focusingofpicosecondpulsedelectricfieldsinhumanbrainmodelwithimpulseradiatingantenna[J].TransactionsofChinaElectrotechnicalSociety.2016，31(3)：195-202.

[10]SharabiS，LastD，GuezD，etal.Dynamiceffectsofpointsourceelectroporationontheratbraintissue[J].Bioelectrochemistry，2014，99(10)：30-39.

[11]DevB，RabussayD，WideraG，etal.Medicalapplicationsofelectroporation[J].IEEETransactionsonPlasmaScience，2000，28(1)：206-223.

[12]PotterH.Electroporationinbiology：methods，applications，andinstrumentation[J].AnalyticalBiochemistry，1988，174(2)：361-373.

[13]WeaverC，ChizmadzhevA.Theoryofelectroporation：areview[J].BioelectrochemandBioenerg，1996，41(2)：135-160.

[14]ZimmermannU，F(xiàn)riedU，MussauerH，etal.Electromanipulationofmammaliancells：fundamentalsandapplication[J].IEEETransactionsonPlasmaScience，2000，28(1)：72-82.

[15]EsserT，SmithC，GowrishankarR，etal.Towardssolidtumortreatmentbyirreversibleelectroporation：intrinsicredistributionoffieldsandcurrentsintissue[J].TechnologyinCancerResearchandTreatment，2007，6(4)：261-273.

[16]YaoC，SunC，MiY，etal.Experimentalstudiesonkillingandinhibitingeffectsofsteeppulsedelectricfield(SPEF)totargetcancercellandsolidtumor[J].IEEETransactionsonPlasmaScience，2004，32(4)：1626-1633.

[17]NealR，GarciaA，KavnoudiasH，etal.Invivoirreversibleelectroporationkidneyablation：experimentallycorrelatednumericalmodels[J].IEEETransactionsonBiomedicalEngineering，2015，62(2)：561-569.

[18]MandelY，LauferS，BelkinM，etal.Irreversibleelectroporationofhumanprimaryuvealmelanomainenucleatedeyes[J].PLoSONE，2013，8(9)：1-7.

[19]米彥，姚陳果，李成祥，等.基于場-路復合模型的細胞內(nèi)外膜跨膜電位時頻特性[J].電工技術(shù)學報.2011，26(2)：14-21.MiYan，YaoChenguo，LiChengxiang，etal.Time-frequencycharacteristicsoftransmenbranepotentialsoncellularinnerandoutermembranesbasedondielectric-circuitcompoundmodel[J].TransactionsofChinaElectrotechnicalSociety，2011，26(2)：14-21.

[20]姚陳果，王劍飛，李成祥，等.脈沖電場誘導細胞脂雙層膜電穿孔的仿真分析[J].高電壓技術(shù)，2010，36(2)：423-427.YaoChenguo，WangJianfei，LiChengxiang，etal.Simulationanalysisofelectroporationinducedbyintensenanosecondpulseelectricfieldincellmembranebilayers[J].HighVoltageEngineering,2010，36(2)：423-427.

[21]DeBruinK，KrassowskaW.Modelingelectroporationinasinglecell.i.effectsoffieldstrengthandrestpotential[J].BiophysicalJournal，1999，77(3)：1213-1224.

[22]KrassowskaW，F(xiàn)ilevP.Modelingelectroporationinasinglecell[J].BiophysicalJournal，2007，92(2)：404-417.

[23]JiangC，DavalosR，BischofC.Areviewofbasictoclinicalstudiesofirreversibleelectroporationtherapy[J].IEEETransactionsonBiomedicalEngineeing，2014，62(1)：4-20.

[24]PuciharG，MiklavcˇicˇD，KotnikT.Atime-dependentnumericalmodeloftransmembranevoltageinducementandelectroporationofirregularlyshapedcells[J].IEEETransactionsonBiomedicalEngineering，2009，56(5)：1491-1501.

[25]HibinoM，ShigemoriM，ItohH，etal.Membraneconductanceofanelectroporatedcellanalyzedbysubmicrosecondimagingoftransmembranepotential[J].BiophysicalJournal，1991，59(1)：209-220.

[26]BarnettA.Thecurrent-voltagerelationofanaqueousporeinalipidbilayermembrane[J].BiochimicaBiophysicalActa，1990，1025(1)：10-14.

[27]PuciharG，KotnikT，ValicˇB，etal.Numericaldeterminationofinducedtransmembranevoltageinducedonirregularlyshapedcells[J].AnnalsofBiomedicalEngineering，2006，34(4)：642-652.

[28]NeuJ，KrassowskaJ.Asymptoticmodelofelectroporation[J].PhysicalReviewE，1999，59(3)：3471-3482.

[29]DeBruinK，KrassowskaW.Electroporationandshockinducedtransmembranepotentialinacardiacfiberduringdefibrillationstrengthshocks[J].AnnalsofBiomedicalEngineering，1998，26(4)：584-596.

[30]SmithK，NeuJ，KrassowskaW.ModelofcreationandevolutionofstableelectroporesforDNADelivery[J].BiophysicalJournal，2004，86(5)：2813-2826.

[31]TonerM，CravalhoE.Kineticsandlikelihoodofmembraneruptureduringelectroporation[J].PhysicsLettersA，1990，143(8)：409-412.

[32]NeuJ，SmithJ，KrassowskaW.Electricalenergyrequiredtoformlargeconductingpores[J].Bioelectrochemistry，2003，60(1-2)：107-114.

[33]JiangC，QinZ，BischofJ.Membrane-targetingapproachesforenhancedcancercelldestructionwithirreversibleelectroporation[J].AnnalsofBiomedicalEngineering，2014，42(1)：193-204.

[34]HibinoM，ItohH，KinositaK.Timecoursesofcellelectroporationasrevealedbysubmicrosecondimagingoftransmembranepotential[J].BiophysicalJournal,1993，64(6)：1789-1800.

[35]ChangD，ReeseT.Changesinmembranestructureinducedbyelectroporationasrevealedbyrapid-freezingelectronmicroscopy[J].BiophysicalJournal，1990，58(1)：1-12.

[36]DemiryurekY，NickaeenM，ZhengM，etal.Transport，resealing，andre-porationdynamicsoftwo-pulseelectroporation-mediatedmoleculardelivery[J].BiochimicaetBiophysicaActa(BBA)-Biomembranes，2015，1848(8)：1706-1714.

[37]PakhomovAG，GianulisE，VernierPT，etal.Multiplenanosecondelectricpulsesincreasethenumberbutnotthesizeoflong-livednanoporesinthecellmembrane[J].BiochimicaetBiophysicaActa(BBA)-Biomembranes，2015，1848(4)：958-966.

[38]TekleE，AstumianRD，F(xiàn)riaufWA，etal.AsymmetricporedistributionandlossofmembranelipidinelectroporatedDOPCvesicles[J].BiophysicalJournal，2001，81(2)：960-968.

[39]NeuJC，KrassowskaW.Modelingpostshockevolutionoflargeelectropore[J].PhysicalReviewE，2003，67(2)：1-12.

Simulation Analysis on Dynamic Process of Electroporation by the Model Based on Energy Probability and Pore Force in Cell Exposed to Pulsed Electric Field

Yao Chenguo Lü Yanpeng Zhao Yajun Dong Shoulong Liu Hongmei

(StateKeyLaboratoryofPowerTransmissionEquipment&SystemSecurityandNewTechnologyChongqingUniversityChongqing400030China)

Matlabisusedtofornumericalanalysisonthetransmembranevoltage，theradiusofthemicro-pores，poreforce，generationandevolutionofporesinasphericalcellexposedto1kV/cm、100μspulsedelectricfield(PEF).Thismodeldifferstraditionalmodelintwodistinctways.Atfirst，thisstudyusestheenergyprobabilitymodeltodecidetheinitialtimeofporeexpansion.Whentheprobabilitythattheinjectedenergyexceedstheminimumenergyofporeisgreaterthan0.99，theporecanexpand.Secondly，poreforcemodelisusedtodescribethedevelopmentanddistributionofporeinthecellmembrane，andanewparameter-poremomentumisputforwardtodecidethesizeofthepores.Resultsshowthatelectroporationhappensonmostdistributionsexceptthepolesregions.Poreforceisclosedtozeroatmostregionswhichmeanthebalanceoccursontheseregions.However，poreforcenearpolesregionsisstillpositive，soporesareexpandingifthetimeofthepulseisextended.Thismodelcanbeausefultooltodecideporeformationandpredictionofporedevelopment.

Electroporation，transmembranevoltage，energyprobability，poreforce

國家自然科學基金項目(51477022)、重慶市研究生科研創(chuàng)新項目(CYB16018)、國家創(chuàng)新研究群體基金項目(51321063)、中央高?；究蒲袠I(yè)務費項目(106112015CDJZR158804)和重慶市杰出青年基金項目(cstc2014jcyjjq90001)資助。

2015-07-05 改稿日期2015-11-12

TM151；TM836；Q64

姚陳果 男，1975年生，教授，博士生導師，研究方向為電氣設(shè)備在線監(jiān)測與故障診斷技術(shù)、生物醫(yī)學中的電工新技術(shù)及高電壓新技術(shù)。

E-mail：yaochenguo@cqu.edu.cn(通信作者)

呂彥鵬 男，1991年生，博士研究生，研究方向為脈沖功率技術(shù)及其生物醫(yī)學應用。

E-mail：lvyanpeng@cqu.edu.cn