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

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

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

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

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

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

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

1 引言

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

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

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

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

2 理論和方法

2.1 電導(dǎo)率模型

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

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

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

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

2.2 不確定性量化

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

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

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

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

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

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

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

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

2.4 消融深度評(píng)估

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

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

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

3 結(jié)果和討論

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

圖2 所示為σH和σG的平均值、標(biāo)準(zhǔn)差和90%預(yù)測區(qū)間.IRE 作用導(dǎo)致的兩個(gè)模型的指數(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%達(dá)到0.434 S/m.σH和σG的標(biāo)準(zhǔn)差和90%預(yù)測空間展現(xiàn)出很不相同的趨勢.σH的90%預(yù)測空間相比于σG要小很多,這表明σH的輸出結(jié)果受參數(shù)不確定性影響較小.

圖2 H 模型和G 模型的電導(dǎo)率輸出的平均值、標(biāo)準(zhǔn)差和90%預(yù)測空間Fig.2.Mean value,standard deviation,and 90% prediction space of conductivity output of H and G models.

3.2 兩種電導(dǎo)率模型的敏感度分析

圖3 所示為兩種電導(dǎo)率模型的一階Sobol 敏感參數(shù)和其平均值.從圖3 可以看出,參數(shù)敏感度與電場強(qiáng)度相關(guān).對(duì)于H 模型,初始電導(dǎo)率σ0的一階Sobol 指數(shù)在各電場強(qiáng)度下都顯示出較大值(>0.3),并且其平均值也最大,是H 模型最敏感的參數(shù).其次E1的不確定性也對(duì)H 模型有較大影響.與之相對(duì)的,G 模型的各個(gè)參數(shù)敏感度分布較為分散.Em在其中較為敏感.這種模型輸出由多個(gè)低敏感參數(shù)共同決定的情況會(huì)導(dǎo)致輸入?yún)?shù)微小變化產(chǎn)生較大的輸出不確定性.在一項(xiàng)針對(duì)不同因素影響脈沖消融的仿真研究中,結(jié)果發(fā)現(xiàn)單個(gè)電極位置的微小變化對(duì)腫瘤消融區(qū)域的影響較微小,但是多個(gè)電極位置的偏差會(huì)導(dǎo)致消融區(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的不確定性輸出由多個(gè)敏感度相近的參數(shù)共同決定,會(huì)更容易受各種條件的影響.這個(gè)結(jié)果為特定應(yīng)用領(lǐng)域選擇模型提供了一種標(biāo)準(zhǔn),對(duì)于PFA 的評(píng)估和規(guī)劃來說,作為公認(rèn)影響脈沖電場分布的σ0如果可以被準(zhǔn)確測量,則使用H 模型;反之,如果σ0無法準(zhǔn)確測量,而其他參數(shù)可以較為準(zhǔn)確地確定,則G 模型更適用.圖3 的模型參數(shù)敏感度分布和排序?yàn)榛趯?shí)驗(yàn)數(shù)據(jù)的模型參數(shù)擬合求解提供了降低計(jì)算消耗的方法.在實(shí)際求解過程中,可以根據(jù)體系的先驗(yàn)知識(shí)固定低敏感參數(shù)的值,再求解高敏感參數(shù).

3.3 組織電導(dǎo)率對(duì)消融深度的影響評(píng)估

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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