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

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

(輸配電裝備及系統(tǒng)安全與新技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室(重慶大學(xué)) 重慶 400030)

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

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

(輸配電裝備及系統(tǒng)安全與新技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室(重慶大學(xué)) 重慶 400030)

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

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

0 引言

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

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

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

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

1 解析模型

1.1 數(shù)值模型

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

(1)

(2)

細(xì)胞模型如圖1所示。

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

1.2 跨膜電勢(shì)

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

neJe=nmJm(t)=niJi

(3)

Je=λeEe

(4)

Vm=φe-φi

(5)

(6)

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

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

表1 模型參數(shù)

Tab.1Modelparameters

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

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

JEP(t)=iEPN(t)

(7)

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

(8)

式中

(9)

(10)

(11)

(12)

在電穿孔過程中，細(xì)胞膜電導(dǎo)率隨時(shí)間的變化可表示為[27]

(13)

1.3 孔密度

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

(14)

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

1.4 孔徑的發(fā)展

未施加電場(chǎng)前，細(xì)胞膜上所有的微孔都是疏水孔，如圖2[30]所示?？讖降拇笮∮善淦琳夏芰靠刂啤．?dāng)細(xì)胞處于高脈沖場(chǎng)強(qiáng)時(shí)，如果r>r*，大部分疏水孔將會(huì)自發(fā)的轉(zhuǎn)為親水孔，親水孔能夠迅速發(fā)展到最小能量對(duì)應(yīng)的孔徑rm。因此本文假定初始的孔徑為rm。微孔的能量公式為[30]

(15)

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

(16)

(17)

(18)

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

(19)

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

1.5 微孔受力模型

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

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

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

2πγ+2πσeffrjr≥rm

(20)

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

(21)

式中，A為細(xì)胞的表面積；Ap為細(xì)胞膜的上所有孔的面積，其表達(dá)式為

(22)

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

1.6 仿真計(jì)算

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

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

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

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

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

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

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

2.1 跨膜電壓

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

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

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

2.2 微孔半徑與受力

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

(23)

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

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

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

3 討論

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

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

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

文獻(xiàn)[16]中，脈沖處理過后的細(xì)胞在掃描電子顯微鏡下觀察發(fā)現(xiàn)了微米級(jí)別的微孔。文獻(xiàn)[38]用人工膜代替細(xì)胞膜進(jìn)行脈沖處理后，也在人工膜上發(fā)現(xiàn)了微米級(jí)別的微孔。而在本文仿真中，極點(diǎn)附近出現(xiàn)了微米級(jí)的微孔，與實(shí)驗(yàn)?zāi)軌蚝芎玫奈呛稀?/p>

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

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

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

4 結(jié)論

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

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

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

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

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

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

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

[3] 米彥，張晏源，儲(chǔ)貽道，等.基于非平衡Blumlein型多層微帶傳輸線的高壓納秒脈沖發(fā)生器[J].電工技術(shù)學(xué)報(bào)，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] 郭飛，李成祥，唐賢倫，等.沖激輻射天線實(shí)現(xiàn)皮秒脈沖電場(chǎng)在人體大腦模型中聚焦的研究[J].電工技術(shù)學(xué)報(bào)，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]米彥，姚陳果，李成祥，等.基于場(chǎng)-路復(fù)合模型的細(xì)胞內(nèi)外膜跨膜電位時(shí)頻特性[J].電工技術(shù)學(xué)報(bào).2011，26(2)：14-21.MiYan，YaoChenguo，LiChengxiang，etal.Time-frequencycharacteristicsoftransmenbranepotentialsoncellularinnerandoutermembranesbasedondielectric-circuitcompoundmodel[J].TransactionsofChinaElectrotechnicalSociety，2011，26(2)：14-21.

[20]姚陳果，王劍飛，李成祥，等.脈沖電場(chǎng)誘導(dǎo)細(xì)胞脂雙層膜電穿孔的仿真分析[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

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

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

TM151；TM836；Q64

姚陳果 男，1975年生，教授，博士生導(dǎo)師，研究方向?yàn)殡姎庠O(shè)備在線監(jiān)測(cè)與故障診斷技術(shù)、生物醫(yī)學(xué)中的電工新技術(shù)及高電壓新技術(shù)。

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

呂彥鵬 男，1991年生，博士研究生，研究方向?yàn)槊}沖功率技術(shù)及其生物醫(yī)學(xué)應(yīng)用。

E-mail：lvyanpeng@cqu.edu.cn