量化電導率模型參數(shù)多樣性導致的脈沖電場消融預測的不確定性*

莊杰 韓瑞 季振宇 石富坤4)?

1) (中國科學技術(shù)大學生物醫(yī)學工程學院(蘇州),生命科學與醫(yī)學部,蘇州 215000)

2) (中國科學院蘇州生物醫(yī)學工程技術(shù)研究所,蘇州 215163)

3) (空軍軍醫(yī)大學軍事生物醫(yī)學工程學系,西安 710032)

4) (濟南國科醫(yī)工科技發(fā)展有限公司,濟南 250101)

脈沖電場消融(pulsed field ablation,PFA)以其非熱、消融邊界清晰等優(yōu)勢成為腫瘤和組織消融的關(guān)注重點.然而,組織生物物理特性和脈沖參數(shù)的多樣性導致PFA 效果難以預測.本研究考慮PFA 消融閾值和兩個常用電導率模型參數(shù)的不確定性,評估其對消融深度的影響.從模型復雜度、參數(shù)敏感度分布、模型魯棒性給出了模型選擇的指標.結(jié)果表明單因素決定的模型具有更強的魯棒性.消融閾值和組織初始電導率依次是影響消融深度評估的兩個最敏感參數(shù).該策略可用于多種數(shù)學物理和仿真模型,以增強模型可信度和簡化模型.

1 引言

細胞的生物物理特性的異質(zhì)性導致電場分布發(fā)生改變,是影響電場治療效率的關(guān)鍵因素[1-3].為了避免欠治療或過度治療,近期越來越多的工作聚焦在根據(jù)特定病人的特異性開發(fā)基于電場分布的手術(shù)規(guī)劃系統(tǒng)[4-6].但是生物多樣性導致的模型參數(shù)的不確定性如何影響高壓電場分布還缺乏研究.

強電場作用導致細胞膜通透性發(fā)生變化[7,8],即電穿孔,從而導致細胞和組織電導率發(fā)生與電場相關(guān)的變化[1].已有研究表明電導率等電學性質(zhì)的變化一方面可以反映生物體生理病理特征的變化[9,10],另一方面也會影響治療過程中電場分布的改變[1],并影響對治療電場的預測和實際治療效果.Sel 等[2,11]在腫瘤消融模型中發(fā)現(xiàn)基于Sigmoid 函數(shù)描述腫瘤電導率隨電場的變化,消融模型預測的消融面積與實驗結(jié)果較固定電導率更相符.Garcia 等[3]基于光滑Heaviside 函數(shù)研究脈沖電場腫瘤消融的治療規(guī)劃,發(fā)現(xiàn)組織電導率和消融閾值是影響治療效果預測的最重要因素.Neal 等[5]發(fā)現(xiàn)基于Sigmoid的Gompertz 函數(shù)可以很好地擬合電穿孔導致的非線性電導率變化,并發(fā)現(xiàn)電導率對電場分布的影響比溫度更大.Zhao 和Davalos[12]在基于光滑Heaviside 函數(shù)的電導率模型中增加脈沖數(shù)作為自變量,使得電導率模型更符合實驗數(shù)據(jù).姚陳果等[13]采用Gompertz 函數(shù)作為組織電穿孔下動態(tài)電導率模型研究了其對消融效果的影響,發(fā)現(xiàn)動態(tài)電導率導致的消融區(qū)域比靜態(tài)電導率更大.Perera-Bel 等[4]基于電穿孔導致組織電導率非線性變化開發(fā)了用于電穿孔治療的手術(shù)規(guī)劃系統(tǒng),可以有效提高手術(shù)治療效果.相關(guān)報道證明了電導率模型的選擇和模型參數(shù)對電場分布和治療效果預測的影響一直是相關(guān)研究的重點.然而現(xiàn)有研究沒有對比不同電導率模型對脈沖電場分布的影響,也沒有考慮模型各參數(shù)的變化對脈沖電場消融效果的影響.

不確定性量化(uncertainty quantification,UQ)目的在于描述輸入的不確定性如何傳導和影響模型輸出.而敏感度分析(sensitivity analysis,SA)用于衡量各個參數(shù)影響模型輸出的敏感度.UQ 一般包含輸出平均值、方差和預測空間,可以提高模型預測準確性,其已廣泛用于地震、氣候預報等領(lǐng)域[14].SA 根據(jù)各模型參數(shù)變化對模型輸出的方差的貢獻可以比較不同參數(shù)間的敏感度,從而幫助降低模型復雜度、優(yōu)化實驗設(shè)計.參數(shù)敏感度分析已經(jīng)被廣泛應(yīng)用于簡化電化學模型[15-17].

本文針對近5 年廣受關(guān)注的脈沖電場消融[18,19],構(gòu)建三維肺靜脈隔離(pulmonary vein isolation,PVI)模型,研究電導率模型參數(shù)和消融閾值不確定性對脈沖電場分布和消融深度的影響.針對電導率模型的參數(shù)敏感度分析,研究發(fā)現(xiàn)最敏感參數(shù)的個數(shù)與模型魯棒性相關(guān).并且可以根據(jù)模型復雜度、參數(shù)敏感度分布、模型魯棒性多個角度選擇模型.在PVI 模型上的研究表明,組織電導率和消融閾值的準確測量是提升脈沖電場消融預測的關(guān)鍵因素,尤其消融閾值是最大敏感參數(shù).

2 理論和方法

2.1 電導率模型

本文采用由Gompertz 曲線(G 模型)和一個光滑Heaviside 函數(shù)(H 模型)描述的兩種電導率模型[2,6].G 模型的電導率σG(E) 可以寫成[20]:

其中,σ0是處理前組織的初始電導率,σm是電穿孔后生物組織電導率達到最大的值.電場E是一個代表施加電場大小的矩陣.另外,指數(shù)A=(Em+Ei)/2,B=(Em-Ei)/2,其中Em和Ei分別代表不可逆和可逆電穿孔閾值;指數(shù)D是由實驗決定的數(shù)值;T0和T分別代表治療前后的溫度,α是溫度系數(shù).

在H 模型中,模型電導率σH(E),如(2)式所示采用一個連續(xù)二階可導的光滑Heaviside 函數(shù)flc2hs 表示:

式中,Emid=(E1+E0)/2,Emid-range=(E1-E0)/2,其中E0和E1分別為達標觸發(fā)電穿孔和達到飽和的電場閾值;參數(shù)W是一個與脈寬正比例相關(guān)的因子[12].光滑Heaviside 函數(shù)flc2hs 的具體形式定義形式可參考文獻[6].

2.2 不確定性量化

不確定性量化用于量化模型輸入不確定性是如何傳遞到模型輸出的.模型輸出的不確定性一般采用平均值、方差和90%預測空間表示.90%預測空間表示有5%的輸出值在這個空間以上或以下.UQ 的計算采用擬蒙特卡羅(Quasi-Monte Carlo)方法.相比于傳統(tǒng)蒙特卡羅方法,這種方法通過采用一個低差異序列選擇樣本從而提高計算速度,減少模型評價需要的數(shù)目[21].

在腫瘤和組織上的電導率測試表明其基本遵從一個正態(tài)分布[22].為了計算方便,其他參數(shù)也采用正態(tài)分布.如表1 所示,H 模型電導率的平均值采用Zhao 等[6]針對5 μs 脈寬的研究,其不確定性大概是10%平均值.而G 模型電導率均值采用(1)式擬合添加了噪聲的H 模型.在UQ 分析中,各個參數(shù)的不確定性設(shè)定為平均值的10%.

表1 不同模型的參數(shù)Table 1. Parameters for different models.

2.3 參數(shù)敏感度分析

常用基于方差的Sobol 敏感度分析方法進行敏感度分析[23].Sobol 方法利用一個敏感度因子表示一個參數(shù)的不確定性對模型輸出的方差貢獻,是一種全局方法[24].Sobol 因子在0—1 之間,越接近1,表示越敏感.此處,采用一階敏感度因子Θi描述參數(shù)的敏感度:

其中,E(O|Pi) 代表當參數(shù)Pi固定時,模型輸出的期望值.Δ(O) 是輸出值的方差,衡量輸出值在平均值Σ附近的變化.

利用Sobol 因子的平均敏感度可以更加明顯地對比不同參數(shù)的敏感度大小.平均敏感度如下所示:

P-i代表除了參數(shù)Pi以外的其他不確定參數(shù).

2.4 消融深度評估

本文通過使用有限元法求解一個靜電場近似方程獲得電場分布.雖然(1)式和(2)式也具有電穿孔閾值的定義,但是為了方便與其他研究對比,將400 V/cm 作為消融閾值確定消融深度[25,26].圖1所示為1/6 的三維PVI 消融模型.模型由12 個緊貼肺靜脈口的電極環(huán)陣列組成,肺靜脈口近似為喇叭形.這種電極形式類似于臨床上使用的電極布局[27].每個電極材料定義為金(Au),長4 mm,相鄰兩個電極間距3 mm,電極截面直徑1.4 mm.PVI 內(nèi)充滿血液.PVI 電導率采用(3)式或(4)式定義.

圖1 1/6 三維PVI 脈沖電場消融模型Fig.1.One-sixth three-dimensional PVI pulsed electric field ablation model.

為了研究仿真模型消融深度的不確定性輸出,初始電導率σ0的均值設(shè)定為0.230 S/m,標準差為均值的10%或者20%.為了研究消融閾值對消融深度的影響,本文消融閾值Et設(shè)定為一個從360—440 V/cm 的均勻分布[28].

3 結(jié)果和討論

3.1 H 模型σH 和G 模型電導率σG 的不確定性量化結(jié)果

圖2 所示為σH和σG的平均值、標準差和90%預測區(qū)間.IRE 作用導致的兩個模型的指數(shù)增長階段從800—1500 V/cm.σH的平均值從0.244 S/m增大到0.451 S/m.σG最小值0.268 S/m 出現(xiàn)在400 V/cm,在1800 V/cm 增長62%達到0.434 S/m.σH和σG的標準差和90%預測空間展現(xiàn)出很不相同的趨勢.σH的90%預測空間相比于σG要小很多,這表明σH的輸出結(jié)果受參數(shù)不確定性影響較小.

圖2 H 模型和G 模型的電導率輸出的平均值、標準差和90%預測空間Fig.2.Mean value,standard deviation,and 90% prediction space of conductivity output of H and G models.

3.2 兩種電導率模型的敏感度分析

圖3 所示為兩種電導率模型的一階Sobol 敏感參數(shù)和其平均值.從圖3 可以看出,參數(shù)敏感度與電場強度相關(guān).對于H 模型,初始電導率σ0的一階Sobol 指數(shù)在各電場強度下都顯示出較大值(>0.3),并且其平均值也最大,是H 模型最敏感的參數(shù).其次E1的不確定性也對H 模型有較大影響.與之相對的,G 模型的各個參數(shù)敏感度分布較為分散.Em在其中較為敏感.這種模型輸出由多個低敏感參數(shù)共同決定的情況會導致輸入?yún)?shù)微小變化產(chǎn)生較大的輸出不確定性.在一項針對不同因素影響脈沖消融的仿真研究中,結(jié)果發(fā)現(xiàn)單個電極位置的微小變化對腫瘤消融區(qū)域的影響較微小,但是多個電極位置的偏差會導致消融區(qū)域的明顯下降[28].而溫度參數(shù)雖然敏感度較高,但由于溫度系數(shù)α敏感度極低,可見H 模型和G 模型中溫度整體并不產(chǎn)生明顯影響.溫度不敏感與不可逆電穿孔的非熱效應(yīng)一致,都是由于脈沖電場作用不產(chǎn)生明顯溫升的結(jié)果[4].

圖3 H 模型和G 模型的各參數(shù)一階Sobol 敏感度指數(shù)和平均Sobol 指數(shù)Fig.3.First order Sobol sensitivity index and average Sobol index of each parameter of the H and G models.

從圖2 可看出,σH的魯棒性比σG更高,這是由于σG的不確定性輸出由多個敏感度相近的參數(shù)共同決定,會更容易受各種條件的影響.這個結(jié)果為特定應(yīng)用領(lǐng)域選擇模型提供了一種標準,對于PFA 的評估和規(guī)劃來說,作為公認影響脈沖電場分布的σ0如果可以被準確測量,則使用H 模型;反之,如果σ0無法準確測量,而其他參數(shù)可以較為準確地確定,則G 模型更適用.圖3 的模型參數(shù)敏感度分布和排序為基于實驗數(shù)據(jù)的模型參數(shù)擬合求解提供了降低計算消耗的方法.在實際求解過程中,可以根據(jù)體系的先驗知識固定低敏感參數(shù)的值,再求解高敏感參數(shù).

3.3 組織電導率對消融深度的影響評估

從3.2 節(jié)分析可以得到H 模型僅有一個最敏感參數(shù),即為初始電導率,具有很好的魯棒性.因此本節(jié)采用H 模型電導率來研究三維PVI 模型的消融電場分布.初始電導率σ0的不確定性設(shè)定為平均值的10%或20%,用于對比不同不確定性下電場分布的不同.20%是根據(jù)實驗確定的范圍[12].

圖4 所示為三維PVI 模型上的電場分布仿真結(jié)果,圖4(a)是完整模型的俯視圖,圖4(b)是1/6結(jié)果圖.紅色區(qū)域是電場分布區(qū)域.在電極中間,電場強度較低.這點與前期研究仿真結(jié)果類似,電極邊緣電場較強而電極中間電場較弱[29].消融深度如圖4(a)中的白線定義,從電極表面到400 V/cm處的距離.

圖4 三維PVI 消融模型的電場分布Fig.4.Electric field distribution of three-dimensional PVI ablation model.

為了研究初始電導率σ0的不確定性如何影響消融深度,本文在方差為均值的10%和20%的兩個σ0樣本空間中分別隨機選擇了200 組數(shù)據(jù)進行消融深度的仿真計算.圖5 所示為消融深度數(shù)值個數(shù)分布和累積概率,從圖5 可見消融深度基本遵循高斯分布.方差10%的情況下,91.5%的樣本消融深度為3.675 mm.200 個樣本的消融深度的標準差為0.096 mm,是均值的2.6%.相對比而言,20%的情況,61.0%的樣本具有3.676 mm 的消融深度,其方差為0.172 mm,是均值的4.7%,達到10%情況的1.8 倍.從累積概率也可以看出20%的情況下作為模型輸出的消融深度呈現(xiàn)更加分散的分布.消融深度的概率分布圖表明消融深度分布隨初始電導率變化的趨勢不是線性的.

圖5 改變σ0,不確定度分別為10% ((a),(c))和20% ((b),(d))時,消融深度的統(tǒng)計結(jié)果Fig.5.The statistical results of the ablation depth by changing σ0 with uncertainty of 10% ((a),(c)) and 20%((b),(d)).

在實際的PFA 治療中,由于脈寬、頻率、極性、電極布局等治療條件的不同,組織消融閾值都會有變化.為研究消融閾值的不確定性對脈沖電場消融深度的影響,在目前的三維PVI 模型中,將消融閾值Et設(shè)定為360—440 V/cm 的均勻分布,在平均值400 V/cm 的10%范圍.圖6 所示為PVI初始電導率為0.230 S/m 情況下,消融深度的分布.圖6(a)所示為消融深度的分布為3.45—3.9 mm.從圖6(b)可見,均勻分布的消融閾值情況下,消融深度遵循高斯分布,平均值為3.681 mm,標準差為0.307 mm (均值的8.3%).消融深度的正態(tài)分布是中心極限定理的結(jié)果.與圖5 的σ0導致的消融深度變化相比,可以看出消融閾值Et是影響消融深度評估的最主要因素.

圖6 消融閾值的均勻分布導致消融深度呈高斯分布Fig.6.The uniform distribution of the ablation threshold results in a Gaussian distribution of the ablation depth.

在一項經(jīng)同行評議的用于脈沖電場腫瘤消融手術(shù)規(guī)劃的仿真研究中,Kos 等[28]探討了腫瘤電導率和消融閾值變化對消融面積的影響.作者觀察到消融閾值是影響消融面積的最大因素,消融閾值增加10%會導致腫瘤消融面積降低超過10%.消融閾值的確定具有非常大的挑戰(zhàn)性,它與細胞和組織的生物物理特性以及脈沖參數(shù)都息息相關(guān)[30].這導致在實際治療過程中,往往只能根據(jù)動物實驗和臨床預實驗等估計消融閾值.本文提出的UQ 分析,通過考察消融閾值的參數(shù)空間可以增強術(shù)前對消融區(qū)域的評估,從而提升治療效果.在Kos 等[28]的研究中組織電導率對消融覆蓋范圍的影響比消融閾值小,并且只有在電導率變化超過30%才出現(xiàn)可見的影響.這可能是由于在他們的仿真模型中采用的電導率是不隨電場變化的常數(shù).這些與本文相似的結(jié)論驗證了本文結(jié)果的正確性.同時本文方法提供了關(guān)于電導率模型的全面的參數(shù)敏感度分布.并進一步在更精細的參數(shù)空間探索影響消融深度的最敏感參數(shù).

圖7 總結(jié)了UQ 和SA 用于模型分析的流程,尤其針對有多個模型可選的研究.本文提出的方法可以方便、高精度地通過模型各參數(shù)的輸入空間和模型輸出不確定性空間判斷模型的魯棒性和敏感度參數(shù),提升模型的優(yōu)化計算和預測準確率.相比于先前研究中常用的OFAT,本文方法具有參數(shù)連續(xù)性高、可全局計算等優(yōu)點.

圖7 用于多模型評估的UQ 和SA 示意圖Fig.7.Schematic diagram for UQ and SA for multi-model evaluation.

4 結(jié)論

針對脈沖電場消融相關(guān)的電導率模型和三維PVI 電場仿真模型,計算了模型不確定性輸入對模型輸出的影響.1)不同電導率模型在共有參數(shù)具有同樣比例的變化時表現(xiàn)出不同的魯棒性;2) H模型輸出不確定性主要由初始電導率決定,展現(xiàn)出更強的魯棒性,而G 模型由多個敏感度相近的參數(shù)影響其不確定性輸出;3)在消融深度評估模型中,σ0變化越大導致的消融深度變化也越大,但呈現(xiàn)非線性比例關(guān)系;4) 與電場相關(guān)的非線性電導率模型中初始電導率對消融深度的影響比固定電導率的影響更大;5) 影響消融深度的最主要參數(shù)是消融閾值.

本文采用的方法以直觀、全局、精細的方式給出了模型參數(shù)對模型輸出的影響和排序,可以有效提高模型預測準確性,降低計算消耗,挑選適用模型.這些方法具有普適性,可用于各種參數(shù)化的數(shù)學物理模型,增強對研究系統(tǒng)的理解和優(yōu)化實驗設(shè)計.