不可逆電穿孔多因素動態(tài)電導率模型的研究與分析

唐 瀟 姚陳果 鄭 爽,2 劉紅梅 董守龍

（1 輸配電裝備及系統(tǒng)安全與新技術國家重點實驗室（重慶大學） 重慶 400030 2.國網(wǎng)重慶市電力公司永川供電分公司 永川 404100）

0 引言

近年來，不可逆電穿孔（Irreversible Electroration, IRE）腫瘤消融技術以其微創(chuàng)、非熱等生物電學效應引起學者的廣泛關注[1-6]，其在腦[7-9]、肝臟[10-13]、胰腺[14]、前列腺[15-16]、腎臟[17]等多種器官局部腫瘤的成功消融，驗證了其應用于臨床治療的安全性與有效性，表現(xiàn)出良好的發(fā)展前景。

準確的術前治療計劃是實現(xiàn)腫瘤安全有效消融的關鍵，也是進一步推進不可逆電穿孔腫瘤消融技術臨床應用所需解決的關鍵難題之一[18-19]。所謂術前治療計劃，即根據(jù)組織電導率模型，結(jié)合患者腫瘤組織幾何及電學特性，模擬計算脈沖電場作用下組織內(nèi)部電場分布及腫瘤消融療效，以確定最佳的治療策略（包括電壓、脈寬、重復頻率、個數(shù)等脈沖參數(shù)及電極位置、暴露長度、施加順序等電極參數(shù)）?？梢?，術前治療計劃的制定主要是基于數(shù)值模擬手段，其所應用的電導率模型的準確性與治療是否過量或欠量直接相關[20]。

在早期研究中，為了降低治療的復雜程度以及數(shù)值模擬的計算量，研究人員往往忽略脈沖作用下組織電導率發(fā)展的動態(tài)過程，考慮組織的電導率是不隨脈沖電場強度、個數(shù)等因素變化的恒定值，即靜態(tài)電導率模型[21]。然而，基于此靜態(tài)電導率模型，忽略了變化的電導率對電場分布的影響，往往使得模擬值比實際消融面積偏小，即可能導致過量治療，損傷周邊正常組織[20]。隨后，研究人員開始考慮脈沖電場作用下組織電導率分布的變化。2005年，D.Sel等基于靜態(tài)電導率模型，考慮脈沖電場強度對電導率的影響，通過對實驗數(shù)據(jù)進行擬合，提出隨電場強度呈S形變化的動態(tài)電導率曲線，并通過實驗電流與仿真電流、實驗消融面積與仿真消融面積的良好匹配證明該曲線的合理性[22]。2010年，P.A.Garcia等采用二階導數(shù)平滑曲線對電導率隨電場強度的變化進行擬合，并考慮溫度對電導率的影響，實現(xiàn)了隨電場強度及溫度變化的動態(tài)電導率模型的建立[23]。然而，從與脈沖作用無關的靜態(tài)電導率模型到現(xiàn)如今發(fā)展較為成熟的隨脈沖電場強度及溫度變化的傳統(tǒng)電導率模型，電導率均采用第一個脈沖作用下電壓、電流數(shù)據(jù)進行求解，既忽略了脈沖個數(shù)的累積效應，也忽略脈沖期間組織特性的動態(tài)發(fā)展過程[24]。2014年，E.M.Dunki- Jacobs等在臨床研究中發(fā)現(xiàn)，隨著脈沖個數(shù)的增加，組織電穿孔程度及范圍逐漸增加并趨于穩(wěn)定[25]；近期，本團隊劉紅梅等利用低壓測量脈沖獲取了馬鈴薯組織在不同電場強度、不同脈沖個數(shù)下電導率的變化，揭示了脈沖作用前后組織電特性變化規(guī)律，也證實了脈沖個數(shù)增加會導致組織電導率的增大，并建議目前的電導率模型應進一步考慮脈沖個數(shù)的影響[26]。因此，基于現(xiàn)有的傳統(tǒng)電導率模型，進一步建立考慮脈沖數(shù)目對組織消融程度及范圍有所影響的動態(tài)電導率模型，對進一步提升數(shù)值模型精確性具有重要意義。

本文基于已有研究成果，在建立隨脈沖電場強度及溫度變化的傳統(tǒng)動態(tài)電導率模型的基礎上，考慮脈沖個數(shù)對組織電導率的影響，建立電場、熱場與脈沖數(shù)目的多因素動態(tài)電導率模型，并通過針電極模型下仿真電流與實驗實測電流、仿真溫升與實測溫升相比較，驗證本文所建立多因素動態(tài)電導率模型的合理性。通過本文多因素動態(tài)電導率模型的建立，有助于進一步理解組織電穿孔的動態(tài)發(fā)展過程，為后續(xù)的術前治療計劃的制定奠定基礎。

1 材料與方法

馬鈴薯屬于茄類植物組織，是評估趨勢的合理模型，常用于脈沖作用下的阻抗分析、消融面積測定的研究，是研究不可逆電穿孔動態(tài)過程的常用組織[27]。因此，本文以馬鈴薯塊莖組織作為實驗對象，為避免不均勻電場引起電導率分布不均勻，采用平板電極施加脈沖。實驗前，樣本均處理成直徑為15mm、厚度為5mm的圓柱體；實驗時，兩個半徑為20mm的圓形金片電極以三明治形式，對馬鈴薯樣品進行脈沖刺激，每組參數(shù)至少重復3次實驗。

實驗所使用的脈沖發(fā)生器由課題組自主研制，用于產(chǎn)生脈寬為 100μs，重復頻率 1Hz的傳統(tǒng) IRE脈沖[29]。實驗中，采用電流傳感器（Pearson Electronics Inc.PaloAlto, USA）與高壓探頭（Teledyne LeCroy, PE5kV, USA）測量高壓脈沖施加電壓與響應電流，采用美國LumaSense Technologies公司生產(chǎn)的FOTLAB.KIT型光纖溫度傳感器測量組織中心位置的溫度，并利用示波器（WavePro 760Zi-A,Teledyne LeCroy Inc., New York, USA）對實驗中電壓、電流及溫度數(shù)據(jù)進行采集，采用頻率500MHz。實驗系統(tǒng)及電壓波形如圖1所示。

圖1 實驗系統(tǒng)及電壓波形Fig.1 Experimental system diagram and voltage waveform diagram

實驗中，所施加的脈沖參數(shù)為脈沖脈寬100μs，脈沖場強從100～1 000V/cm，以100V/cm為梯度增加，脈沖施加個數(shù)均為90，重復頻率均為1Hz。每一個實驗樣品的組織初始溫度均保持為16～17℃，以減小組織間溫度差異引起的初始電導率的顯著不同。

實驗中每組電壓、電流數(shù)據(jù)均取自脈沖尾部500個采用點數(shù)據(jù)的平均值。每組數(shù)據(jù)至少重復測量3次，以盡量減小組織異質(zhì)性帶來的誤差。

2 模型的建立

2.1 傳統(tǒng)電導率模型的建立

隨脈沖電場強度及溫度變化的電導率模型已被廣泛運用于數(shù)值模型的仿真擬合，根據(jù)文獻[30]，隨脈沖電場強度及組織溫度變化的電導率曲線為

其中

式中，E為電場強度；σ為電導率；σ0為初始電導率；T0為初始溫度；T為組織溫度；α為反映溫度升高對電導率的影響系數(shù)，設置其值為 2%[31]；A為電導率增長因子，用于描述電場作用下電導率的增加；Edel為電導率從初始電導率發(fā)展到穩(wěn)定值這一變化區(qū)間的中間值，而Erange則表示整個過渡區(qū)間的變化范圍的1/2；E0與E1分別為電導率開始變化與穩(wěn)定對應的電場強度。flc2hs函數(shù)如圖2所示。

圖2 flc2hs函數(shù)Fig.2 flc2hs function

在隨電場強度及溫度變化的電導率模型中，所測量的組織電導率主要由三部分組成：①脈沖未作用前的初始電導率σ0，表征組織固有介電特征；②脈沖作用過程中發(fā)生電穿孔引起的電導率的增加，主要與脈沖電場強度有關；③由脈沖作用產(chǎn)生焦耳熱引起的電導率增加，主要與溫度變化有關。

本文采用的傳統(tǒng)脈沖的脈寬為 100μs，因此未考慮瞬態(tài)電容電流，脈沖脈寬尾部電壓、電流數(shù)據(jù)均能達到穩(wěn)定，則電導率可表示為

式中，ρ為電阻率；L為組織厚度；S為樣品面積；R為電阻；U為脈沖電壓；I為組織負載電流。在低壓脈沖作用下，電導率變化小，本文將低壓25V時測量得到的平均電導率0.03S/m作為初始電導率[27]。

脈沖作用過程中，通過溫度傳感器實測溫度變化。在電極電壓為250V及500V脈沖作用下，實驗組織溫度變化如圖3所示。在250V脈沖作用時，組織溫度變化不明顯，90個脈沖作用后溫升在1℃左右，而在 500V電壓作用時，溫度隨著脈沖個數(shù)逐漸增大且變化趨勢明顯，最大溫升可達到 4℃左右。由此可知，在高電場強度脈沖作用下組織溫升明顯，對電導率影響也將更大，并且隨脈沖個數(shù)增加，溫度對電導率作用也將增大，由此在電導率模型中考慮溫度的影響是必要的。

圖3 實驗組織溫度變化Fig.3 Experimental tissue temperature change

基于已有模型及實驗參數(shù)，通過Matlab平臺，采用最小二乘法對式（1）擬合，可得到不同脈沖個數(shù)作用下隨電場強度及溫度變化的電導率曲線。圖4為消融效果與電導率曲線擬合。

不同脈沖數(shù)目下電導率擬合曲線與實驗數(shù)據(jù)點擬合程度良好，脈沖個數(shù)N=1～90均建立隨電場強度與溫度的變化規(guī)律，各脈沖數(shù)目下的電導率曲線擬合度R2均在 0.9以上。從圖4c可知，隨脈沖個數(shù)增加，引起電導率變化的電場強度降低，即脈沖數(shù)目的增加，將減小引起電導率變化所需的電場強度。

圖4 消融效果與電導率曲線擬合Fig.4 Ablation effect and conductivity curves fitting

當前廣泛采用的傳統(tǒng)電導率模型均是采用第一個脈沖作用時電壓、電流數(shù)據(jù)進行電導率曲線的擬合，以此作為整個脈沖過程中組織的電導率特性。在本文中，傳統(tǒng)電導率模型根據(jù)式（1），可得到該模型中的參數(shù)A=12.55，Edel=334.353，E0=92.014，即傳統(tǒng)單因素模型可表示為

2.2 多因素動態(tài)電導率曲線

在建立隨脈沖電場強度及溫度變化的電導率曲線基礎上，將脈沖數(shù)目的影響考慮其中。文中通過分別觀測量化動態(tài)電導率模型中變量A、Edel及E0隨脈沖數(shù)目的變化趨勢，建立變量隨脈沖數(shù)目的關系。

在對參數(shù)進行擬合時，需滿足兩個條件：①滿足數(shù)據(jù)點隨脈沖數(shù)目變化的趨勢，與已有數(shù)據(jù)點相貼近；②在脈沖個數(shù)趨于0或無窮時，均滿足理論要求?；谶@兩個條件，對于變量A的變化規(guī)律可進行擬合，即

參數(shù)A隨脈沖個數(shù)變化的擬合曲線如圖5所示。擬合度R2=0.933 2，由A的變化趨勢可以看出，當脈沖數(shù)目N趨近于0時，A值趨近于0，電導率σ(E)=σ0，此時，未施加脈沖，組織未發(fā)生電穿孔。隨著脈沖的施加，組織穿孔程度增加，電導率增加，且存在快速增加與平緩增加變化趨勢，由此可看出，在脈沖個數(shù)較小時，引起的穿孔變化較大，隨著個數(shù)增加，穿孔程度達到飽和，引起的電導率的變化趨于平穩(wěn)狀態(tài)。當個數(shù)N趨近于無窮時，已達到最大穿孔程度，電導率達到穩(wěn)定狀態(tài)，A(N)=15.158，此時電導率穩(wěn)定在 16.158 4σ0。擬合函數(shù)規(guī)律與理論分析可緊密聯(lián)系，證明該函數(shù)存在合理性。

圖5 參數(shù)A隨脈沖個數(shù)變化的擬合曲線Fig.5 The fitting curves of parameter A with the number of pulses

對于變量Edel可擬合為

參數(shù)Edel隨脈沖個數(shù)變化的擬合曲線如圖6所示。擬合度R2=0.997 2，由Edel變化趨勢可以看出，當N趨近于0時，Edel(N)趨近于無窮，在未施加脈沖時電導率不發(fā)生變化；隨著脈沖個數(shù)的增加，組織電導率變化更快，由此Edel值減小；當N趨近于無窮時，Edel(N)趨近于穩(wěn)定值154.287。

圖6 參數(shù)Edel隨脈沖個數(shù)變化的擬合曲線Fig.6 The fitting curves of the parameter Edel with the number of pulses

對于參數(shù)E0，可描述為

參數(shù)E0隨脈沖個數(shù)變化的擬合曲線如圖7所示，擬合度R2=0.941。E0為電導率開始發(fā)生變化對應的最低電場強度，在未施加脈沖時，該值趨近于無窮，隨著脈沖數(shù)目增加，引起電導率變化所需的電場強度降低，并逐漸趨近于0。

圖7 參數(shù)E0隨脈沖個數(shù)變化的擬合曲線Fig.7 The fitting curves of the parameter E0 with the number of pulses

綜上所述，結(jié)合式（4）～式（6）可得，隨脈沖電場強度、溫度及脈沖數(shù)目變化的多因素動態(tài)電導率模型為

其中

該電導率模型在原有考慮脈沖電場強度及溫度影響的基礎上加入脈沖數(shù)目，在溫度為初始溫度時，隨電場強度及脈沖數(shù)目變化的電導率曲面如圖8所示。

圖8 隨脈沖電場強度及個數(shù)變化的電導率曲面Fig.8 Conductivity surface with pulse field strength and number

3 電導率模型的驗證

為驗證所建立的電導率模型的準確性，本文考慮針電極施加脈沖作用于組織誘導其消融的情境下，其對應的仿真模擬的響應電流與實驗實測響應電流相比較來進行驗證。實驗中仍然采用馬鈴薯組織作為實驗對象，并將其處理為半徑 23mm、厚度5mm的圓柱體。電壓分別為500V、700V、1 000V，脈寬為100μs的脈沖電場施加半徑0.5mm、針間距10mm、暴露長度5mm的醫(yī)用不銹鋼針電極。實驗電流通過圖1所示的電流傳感器進行采集。

為獲取模擬電流值，本文借用COMSOL Multyphsics建立1∶1物理尺寸的有限元仿真模型，其幾何模型如圖9所示。

圖9 針電極作用仿真幾何模型Fig.9 Needle electrode simulation model

采用電流模塊及生物傳熱模塊求解其在電場作用下的電熱場分布，具體求解方程如下。

組織內(nèi)的電場分布通過Laplace方程得到

式中，φ為電動勢。

兩電極電動勢分別為

求解域外部邊界條件設置為電絕緣。模型中組織電導率采用本文所建立的模型（見式（7）），針電極采用不銹鋼介電特性。

對于組織內(nèi)的傳熱，采用Pennes生物傳熱方程計算組織內(nèi)部溫升為

式中，k為組織熱導率，k=0.570 2W/(m·K)；T為組織溫度；ρl為組織密度，ρl=1 008kg/m3；cp為組織比熱容，cp=3 161.36J/(kg·K)[1]；此外，根據(jù)實測溫度，傳熱模型中初始溫度設置為16℃。為電流流過組織所產(chǎn)生的焦耳熱，d為一個脈沖的高平時間，τ為脈沖周期，由此將單個脈沖的焦耳熱均分在一個脈沖周期上，簡化計算，且對溫升影響較小。

仿真電流通過對針電極表面法向電流分量進行面積分來求取，即

式中，nJ為某一電極表面的法向電流密度；Is為計算得到的電流。

在不同電壓等級及不同樣品厚度條件下進行仿真與實驗。不同電壓等級樣本在90個脈沖作用下的仿真電流與實驗電流對比如圖10所示。

圖10 針電極實驗電流與仿真電流對比Fig.10 Comparison of experimental current and simulated current of needle electrode

仿真電流與實驗電流在90個脈沖作用范圍內(nèi)的誤差采用最大相對誤差、最小相對誤差及相對誤差平均值來表述，其中相對誤差及相對誤差平均值分別為

式中，fi為仿真電流值；yi為實測電流值；N為脈沖數(shù)目。由電流變化趨勢及誤差分析可看出，仿真電流與實驗電流匹配良好，變化趨勢相近，不同電壓水平所得到的實驗電流與仿真電流的相對誤差平均值均在 4%以下，證明所建立的多因素電導率模型的合理性。實驗值與仿真值的誤差比較見表1。

表1 實驗值與仿真值的誤差比較Tab.1 Comparison of experimental and simulated values

由于低壓脈沖作用下溫度變化不明顯，難以將實驗測量溫升與仿真溫升相結(jié)合，因此本文取高壓1 000V脈沖作用下的中心點溫度進行溫升說明。90個脈沖作用后，實測溫升平均值為4.46℃，而仿真溫升為 4.35℃，與實測溫度相差 2.45%，由此證明模型所考慮溫度變化與實際脈沖作用過程中的溫升作用相近，具有合理性。

4 討論

4.1 實驗及仿真模型中電導率的變化

Rodolfo等證明脈沖個數(shù)的累積引起組織電流的遞增[33]。同樣，根據(jù)本文實驗數(shù)據(jù)也可直觀觀測出隨著脈沖個數(shù)的變化，電流呈現(xiàn)增長趨勢。圖11為板電極電壓為 250V（電場強度 500V/cm）時前20個脈沖作用下電壓及電流的實驗數(shù)據(jù)?？梢钥闯?，電壓數(shù)據(jù)基本保持不變，而隨著脈沖個數(shù)的增加，電流逐漸增加。從理論上進行分析，高壓窄脈沖作用下細胞穿孔形成電流通路，在下一個脈沖作用前，組織在已有穿孔細胞上進一步形成穿孔或擴大孔徑，由此電流通路進一步增加，電流增大，即電導率增加。因此將脈沖數(shù)目的累積通過電導率模型的方式進行描述是與實際組織變化相貼近的，本文所建立隨脈沖電場強度、溫度及脈沖數(shù)目變化的電導率模型為進一步描述脈沖作用下組織電導率動態(tài)變化過程提供模型基礎，為后續(xù)動態(tài)消融效果的評估奠定仿真基礎。

圖11 電壓250V時的實驗電壓與電流Fig.11 Experimental voltage and current at 250V

本文所建立模型依據(jù)實驗電壓、電流，能夠反映出組織電導率隨脈沖個數(shù)增加的動態(tài)變化。在針電極電壓為1 000V時，仿真模型中分別取脈沖個數(shù)N=1、N=30、N=60及N=90的電導率及溫度分布如圖12所示?？梢钥闯觯瑢τ谀骋幻}沖數(shù)目下的電導率而言，針電極附近高電場強度區(qū)域電導率最大，并隨電場強度變化而變化；而對于不同個數(shù)下的電導率分布而言，同一位置處在電場強度作用下電導率隨脈沖個數(shù)的增加而增強，且在高電場強度區(qū)域，隨個數(shù)變化較為顯著，而在低電場強度區(qū)域，電導率變化較小。由此可看出，仿真模型中電導率不僅隨電場強度變化，且隨脈沖施加個數(shù)的增加逐漸增強。同時，從溫度分布可以看出，隨著脈沖個數(shù)的增加，溫升逐漸增加，體現(xiàn)溫度隨脈沖個數(shù)的積累。綜上，多因素動態(tài)電導率模型能夠體現(xiàn)隨脈沖數(shù)目增加所帶來的電導率的累積與溫度的累積，與實際電導率及溫度變化貼近。

圖12 脈沖個數(shù)為1、30、60及90時電導率及溫度分布Fig.12 Conductivity and temperature distribution when the number of pulses is 1, 30, 60 and 90

4.2 與傳統(tǒng)電導率模型的差異

現(xiàn)在已有的傳統(tǒng)模型考慮了脈沖電場強度及溫度對電導率的影響，基于作用的第一個脈沖的電壓、電流數(shù)據(jù)通過擬合的方式獲得A、Edel及E0的參數(shù)，以此作為整個脈沖作用的電導率曲線。在針電極模型下，將本文建立的隨脈沖電場強度、溫度及脈沖數(shù)目變化的多因素電導率曲線與傳統(tǒng)電導率模型進行對比如圖13所示。由圖13可以看出，在電壓為500V、700V及1 000V時，傳統(tǒng)模型下隨著脈沖個數(shù)變化，電流保持恒定，各電壓等級下，傳統(tǒng)電導率模型 90個脈沖仿真電流與實驗電流的誤差分別為(16.399±5.608)%、(16.853±5.815)%及(16.171±4.549)%。本文建立多因素電導率模型仿真電流隨脈沖數(shù)目的增加有逐漸增大趨勢，且已證明與實驗數(shù)據(jù)更為相近。由此可以看出，在傳統(tǒng)電導率模型基礎上，有所改進的多因素動態(tài)電導率模型實現(xiàn)了減小仿真電流與實驗電流誤差的優(yōu)勢，在傳統(tǒng)模型基礎上有所進步。

圖13 多因素模型與傳統(tǒng)模型仿真電流對比Fig.13 Comparison of multi-factor model and traditional model simulation current

在電壓等級為500V、700V及1 000V的針電極模型中，對兩針電極中心點位置（見圖8的P1）進行兩種模型的溫升分析，在三個電壓等級作用下，多因素電導率模型在90個脈沖作用后的溫升（500V、1.213℃；700V、2.618℃；1 000V、6.037℃）均高于傳統(tǒng)電導率模型仿真溫升（500V、1.213℃；700V、2.291℃；1 000V、4.821℃），由此可認為，在高壓多脈沖治療時，采用傳統(tǒng)動態(tài)電導率模型進行術前治療計劃時，可能對組織劑量估算有所偏差，將造成不必要的組織損傷。本文所建立的多因素電導率模型在原有隨脈沖電場強度及溫度變化電導率模型基礎上實現(xiàn)了時間維度的變化，且證明多因素動態(tài)電導率模型與實驗結(jié)果更一致。由此認為，該模型可用于術前治療計劃中選取更為合理的脈沖數(shù)目，解決對形態(tài)、深度及大小等因素均不相同的腫瘤組織施加固定數(shù)目脈沖的問題，實現(xiàn)更為準確且優(yōu)化的脈沖參數(shù)及電極參數(shù)的選擇。為實現(xiàn)對各脈沖數(shù)目下治療效果評估，接下來將探索在不同脈沖個數(shù)下的消融判別依據(jù)，通過多因素電導率模型與消融判據(jù)相結(jié)合，實現(xiàn)隨脈沖數(shù)目變化的消融效果的仿真驗證，建立完整的動態(tài)消融效果評估體系。

5 結(jié)論

本文通過金板電極施加高壓脈沖作用下的組織電壓、電流數(shù)據(jù)建立隨脈沖電場強度、溫度及脈沖數(shù)目變化的動態(tài)電導率模型，該模型在傳統(tǒng)電導率模型基礎上有所改進，進一步考慮脈沖數(shù)目動態(tài)變化，從一定程度上反映組織隨脈沖數(shù)目的增加電穿孔程度逐漸增加直到穩(wěn)定狀態(tài)的動態(tài)階段；并通過應用于針電極作用實現(xiàn)了仿真電流與實驗電流的良好匹配，證明了本文建立的電導率模型的正確性。同時，通過與傳統(tǒng)動態(tài)電導率模型的對比可看出，本文所建立多因素電導率模型通過將脈沖數(shù)目納入模型考慮范圍內(nèi)，使得模型更加接近實際情況，與傳統(tǒng)模型相比，與實驗數(shù)據(jù)相差更小。因此，通過本文研究，為建立更加精確的術前治療仿真模型建立了基礎，對治療計劃的脈沖參數(shù)最優(yōu)化選擇具有重大意義，有利于實現(xiàn)精準治療的目標。