一種左心室壓力個性化估計模型參數(shù)子集選擇方法

柳 軍， 郝麗玲， 何光宇， 徐禮勝,

(1. 東北大學(xué) 醫(yī)學(xué)與生物信息工程學(xué)院， 遼寧 沈陽 110169； 2. 中國醫(yī)科大學(xué) 生物醫(yī)學(xué)工程系， 遼寧 沈陽 110122；3. 沈陽東軟智能醫(yī)療科技研究院有限公司， 遼寧 沈陽 110167)

心臟中左心室壓力波形(left ventricle pressure waveform, LVPW)直接反映了左心室的血流動力學(xué)特征，在臨床壓力容積環(huán)指標(biāo)中是不可或缺的[1]，對于主動脈瓣狹窄等臨床病癥有重要的診斷和治療意義[2].

目前LVPW的測量主要是有創(chuàng)測量和無創(chuàng)模型估計兩種方法.有創(chuàng)測量是在臨床中通過導(dǎo)管介入手術(shù)實(shí)現(xiàn)，但是病人需要承受一定的創(chuàng)傷和風(fēng)險，并且醫(yī)療成本較高.無創(chuàng)測量的方法主要是通過各種特定背景下的心血管模型估計[3-4]，同時可以適當(dāng)結(jié)合患者的外周無創(chuàng)測量信息.模型估計方法由于無需手術(shù)、患者沒有創(chuàng)傷成為目前的一個研究方向[5].而針對每個不同病人，模型參數(shù)個性化成為近年來的研究熱點(diǎn)[6].

通常一個模型中有很多參數(shù)，在臨床應(yīng)用中不可能對模型中所有參數(shù)進(jìn)行個性化取值，因?yàn)橛械哪Ｐ蛥?shù)可能無法測量或者臨床中的醫(yī)療成本過高，所以應(yīng)當(dāng)挑選出模型中對LVPW有較大影響的參數(shù)進(jìn)行個性化取值.靈敏度分析(sensitivity analysis, SA)方法近年來逐漸應(yīng)用于心血管模型參數(shù)個性化選擇中[7].

本文以一個體循環(huán)模型估計LVPW作為實(shí)例，采用自適應(yīng)稀疏多項式混沌展開(adaptive sparse polynomial chaos expansion,asPCE)算法構(gòu)建出原始模型的元模型(meta model)，采用索貝爾法(Sobol method)進(jìn)行靈敏度指標(biāo)計算，最后遴選出對LVPW有較大影響的模型參數(shù)作為參數(shù)子集.

1 體循環(huán)模型

本文以一個體循環(huán)模型作為估計LVPW的例子，在最近的研究中也有類似的模型用于估計LVPW[8]，模型的具體原理公式的推導(dǎo)在文獻(xiàn)[9]中有詳細(xì)的介紹，本文不再贅述.總體結(jié)構(gòu)上模型由一個左心室單纖維模型作為驅(qū)動，外周血管由Windkessel模型構(gòu)成，主動脈瓣和二尖瓣由二極管模型模擬.體循環(huán)模型構(gòu)成原理圖如圖1所示.

圖1 體循環(huán)模型構(gòu)成原理圖

模型的輸出為LVPW以及4個波形特征，在本文中4個波形特征分別定義為一個心動周期內(nèi)LVPW的平均壓pMB，收縮壓pSB，從二尖瓣關(guān)閉時刻到收縮壓pSB的時間間隔tr，從收縮壓pSB時刻到二尖瓣打開時刻的時間間隔td，具體波形特征在LVPW中的含義如圖2所示.

圖2 左心室壓力波形特征

本文體循環(huán)模型中涉及的模型參數(shù)分為三類：生理類參數(shù)，Windkessel外周模型參數(shù)，單纖維模型心肌纖維參數(shù).表1～表3分別給出模型中各個參數(shù)的單位、描述及其范圍.

表1 體循環(huán)模型中的生理類參數(shù)Table 1 Physiological parameters of systemic circulation model

表2 體循環(huán)模型中的外周參數(shù)Table 2 Peripheral parameters of systemic circulation model

表3 體循環(huán)模型中的心肌纖維參數(shù)Table 3 Myocardial fiber parameters of systemic circulation model

2 模型參數(shù)靈敏度計算

2.1 靈敏度指標(biāo)

本文靈敏度分析采用Sobol法進(jìn)行計算.該方法是基于方差分解[10]的思想，用模型輸出的方差來衡量輸出不確定性，輸出的方差可以分解成偏方差的組合.這些偏方差可以表示每個輸入?yún)?shù)或者多個輸入?yún)?shù)相互影響的不確定性對輸出的影響.

假設(shè)模型輸出輸入關(guān)系為

Y=f(X) .

(1)

其中：Y是模型輸出量;X是模型輸入?yún)?shù)，假設(shè)是M維的向量，可以表示為X=(X1X2…XM).本文計算的靈敏度指標(biāo)包括主靈敏度指標(biāo)和總靈敏度指標(biāo)，主靈敏度指標(biāo)定義如下：

(2)

其中：Si是主靈敏度指標(biāo)；V表示求方差運(yùn)算；E表示求期望運(yùn)算.主靈敏度指標(biāo)提供的信息表明應(yīng)該準(zhǔn)確確定哪個輸入?yún)?shù)來獲得最大的方差降幅.總靈敏度指標(biāo)定義如下：

(3)

其中：ST,i是總靈敏度指標(biāo)；-i表示去除Xi的參數(shù)集合.總靈敏度指標(biāo)用來決定哪個參數(shù)是可以在其數(shù)值范圍內(nèi)固定而對最后的結(jié)果幾乎沒有影響.

2.2 元模型構(gòu)建和靈敏度指標(biāo)計算

2.2.1 元模型構(gòu)成及靈敏度計算

各個參數(shù)的索貝爾靈敏度指標(biāo)最初是由蒙特卡洛統(tǒng)計學(xué)方法[11]求出，由于蒙特卡洛方法需要進(jìn)行大量的模型樣本運(yùn)算，計算量巨大，后來有學(xué)者[12]提出了采用多項式混沌展開來替代原始模型，根據(jù)展開式的系數(shù)來計算索貝爾靈敏度指標(biāo).本文采用自適應(yīng)稀疏多項式混沌展開算法進(jìn)行元模型的構(gòu)建[13]，然后根據(jù)多項式系數(shù)計算主靈敏度和總靈敏度兩個靈敏度指標(biāo).

模型的輸出可以表示為輸入?yún)?shù)的一系列截斷后的有限項正交多項式的組合，本文采用勒讓德正交多項式，輸入?yún)?shù)需要轉(zhuǎn)換到[-1,+1]的范圍內(nèi)，以滿足勒讓德正交多項式的要求.模型的輸出輸入關(guān)系為

(4)

(5)

(6)

多變量互相階數(shù)和定義如下：

(7)

(8)

多項式的單變量階數(shù)和與多變量互相階數(shù)和可以由設(shè)計者設(shè)定，考慮到多項式的高階項在多項式中的影響很小，本文設(shè)定如下：

A{3,3}={α|Degree(α)≤3 and Order (α)≤3} .

(9)

本文多項式組合的系數(shù)是通過最小化模型真值和多項式計算值差的平方的方法求得.在求得多項式系數(shù)后，利用多項式的系數(shù)計算索貝爾靈敏度指標(biāo).元模型的方差通過下式求出：

(10)

其中:Hα是一個規(guī)范化因子，其具體取值取決于多項式的類型.本文采用勒讓德正交多項式，其計算公式如下：

(11)

最后可以計算得到主靈敏度指標(biāo)：

(12)

其中Ai是包括多項式中那些αi是正值而其他階次是零的集合.總靈敏度指標(biāo)的計算如下：

(13)

其中AT,i是包括多項式中所有αi的階次是正值的集合.

2.2.2 asPCE算法構(gòu)建元模型

在構(gòu)建元模型和求解多項式系數(shù)的過程中，本文采用asPCE的算法.該算法的優(yōu)點(diǎn)是可以減少模型運(yùn)算樣本數(shù)量.在構(gòu)建元模型的過程中，有選擇地挑選多項式組合中對輸出結(jié)果有較大影響的項，作為最終元模型多項式組合的組成部分.在算法的迭代過程中，將重復(fù)進(jìn)行兩步重要多項式的挑選操作：前向選擇和后向選擇.在前向選擇中，從零開始逐漸增加多項式的最高階次和互相階次，挑選出可以顯著提高元模型準(zhǔn)確度的部分，構(gòu)成前向選擇的多項式組合；在構(gòu)建完前向選擇的多項式組合后，再進(jìn)行后向選擇，在后向選擇中，去掉那些前向選擇中留下的似乎重要但是并沒有顯著提高元模型準(zhǔn)確度的部分.多項式的最高階次和互相階次從零開始逐漸增大到設(shè)定值，在此過程中，算法需要的模型運(yùn)算樣本點(diǎn)數(shù)量將逐漸增加，前向選擇和后向選擇的操作在不斷地迭代進(jìn)行，直到最后確定元模型所包括的多項式組合.

(14)

(15)

其中fPCEi(→X(i))表示除去樣本X(i)所構(gòu)建的元模型.

2.3 模型樣本運(yùn)算

由于asPCE算法中，求解多項式系數(shù)需要用到線性回歸方法，因此需要提前進(jìn)行模型樣本的運(yùn)算.在通用的多項式混沌展開(general polynomial chaos expansion,gPCE)方法中，至少需要NP次模型運(yùn)算樣本.而在asPCE算法中，模型運(yùn)算樣本的數(shù)量可以大大減少，本文按照gPCE方法進(jìn)行了NP=4 060次模型運(yùn)算，但是實(shí)際asPCE算法迭代中只是用到了其中一部分.

在進(jìn)行模型運(yùn)算之前，模型的輸入?yún)?shù)需要進(jìn)行隨機(jī)取值.本文采用Sobol序列準(zhǔn)隨機(jī)數(shù)生成方法[14]，該方法以2為底數(shù)來形成更精細(xì)的單位間隔，并且可以生成輸入?yún)?shù)的基于均勻分布的隨機(jī)序列.

由于模型的輸入?yún)?shù)是高維隨機(jī)序列生成，在模型運(yùn)算中，模型輸出值有很多情況是不符合實(shí)際人體的生理數(shù)據(jù)，因此，要把這些不正確的模型輸出值排除在最后元模型的構(gòu)建過程中.本文根據(jù)健康人群和各種心血管病人群的特點(diǎn)以及文獻(xiàn)[15]的相關(guān)類似設(shè)置，設(shè)置了以下的模型輸出過濾規(guī)則：

1) 模型輸出的計算是基于連續(xù)多個心動周期逐漸收斂穩(wěn)定，因此設(shè)定相鄰兩個周期的2范數(shù)閾值為0.05，當(dāng)超過15個心動周期并且相鄰2個心動周期波形的2范數(shù)小于0.05，就認(rèn)為模型輸出波形達(dá)到收斂;2) 左心室壓力波形的平均值范圍為1～10 kPa;3) 左心室壓力波形的收縮壓范圍為5～24 kPa;4) 左心室壓力波形在舒張期內(nèi)的壓力范圍為0.3～3 kPa;5) 左心室壓力波形中收縮期時長與舒張期時長的比值范圍為0.5～1.2.

圖3 asPCE算法流程與靈敏度計算

3 靈敏度指標(biāo)結(jié)果及元模型與真實(shí)模型輸出對比

本文模型有27個輸入?yún)?shù)和4個左心室壓力波形特征輸出指標(biāo)，經(jīng)過以上主靈敏度和總靈敏度指標(biāo)的計算，所有輸入?yún)?shù)的靈敏度指標(biāo)計算結(jié)果如圖4所示.從最后的主靈敏度和總靈敏度結(jié)果可以看出，4個左心室壓力波形特征分別對應(yīng)的重要的輸入?yún)?shù)有所區(qū)別.取4個波形特征所對應(yīng)的所有重要參數(shù)的并集作為最后模型輸出左心室壓力波形的重要參數(shù)子集，共計12個參數(shù)，并在表4給出了每個參數(shù)對應(yīng)的4個波形特征的主靈敏度最大值和總靈敏度最大值.

為了檢驗(yàn)元模型與真實(shí)模型的逼近程度，本文進(jìn)行了左心室壓力波形特征的對比，在輸入?yún)?shù)相同的情況下，比較從真實(shí)模型和元模型輸出的左心室壓力波形的4個特征，具體對比結(jié)果如圖5所示.其中所有相關(guān)性圖中，r=0.99,P<0.001.并且4個特征分別對應(yīng)的兩組數(shù)據(jù)沒有顯著差異(P>0.05).從上面的LVPW的波形特征相關(guān)性曲線和一致性結(jié)果圖中，可以清楚地看到元模型與真實(shí)模型具有很好的相關(guān)性和一致性.

表4 對應(yīng)左心室壓力波形特征的重要輸入?yún)?shù)Table 4 Important parameters corresponding to the features of LVPW

圖4 靈敏度分析結(jié)果

4 模型輸出LVPW對比及靈敏度計算方法對比

為了驗(yàn)證當(dāng)模型中不重要參數(shù)取值固定的左心室壓力波形情況，進(jìn)行了左心室壓力波形和波形特征的對比.把模型的參數(shù)分成兩個部分，重要參數(shù)和非重要參數(shù).重要參數(shù)取值相同(隨機(jī)取值)，非重要參數(shù)一種情況隨機(jī)取值，另外一種情況固定取值(經(jīng)驗(yàn)值).在這兩種情況下，對真實(shí)模型輸出的左心室壓力波形和波形特征進(jìn)行對比，結(jié)果如圖6和圖7所示，其中所有的相關(guān)性曲線中，r=0.99,P<0.001.LVPW和波形特征在兩種情況下沒有顯著差異(P>0.05).本文進(jìn)行了3 431次實(shí)驗(yàn)，隨機(jī)取其中一例波形對比圖如圖6所示.經(jīng)過靈敏度選擇后，重要參數(shù)子集的數(shù)量比全部參數(shù)數(shù)量顯著減少，從對比結(jié)果圖中可以看出，在重要參數(shù)取值相同，非重要參數(shù)取值固定與取值隨機(jī)兩種情況下，最后輸出的左心室壓力波形和波形特征都具有很好的一致性.經(jīng)過靈敏度選擇后，重要參數(shù)作為重點(diǎn)關(guān)注的參數(shù)子集，非重要參數(shù)取值固定對LVPW波形的估計幾乎沒有重大影響.

圖5 真實(shí)模型和元模型輸出的波形特征相關(guān)性曲線和對應(yīng)的Bland-Altman圖

圖6 不同參數(shù)取值情況下模型輸出LVPW對比

傳統(tǒng)的靈敏度計算經(jīng)常采用蒙特卡洛方法，其劣勢在于需要巨量的模型計算樣本，計算量大，特別對于參數(shù)數(shù)量較多的模型，計算量呈幾何數(shù)量級增長.采用asPCE算法構(gòu)建原始模型的元模型進(jìn)行靈敏度指標(biāo)的計算，一方面可以對不了解原始模型內(nèi)部原理的研究者提供一個簡單的多項式組合的元模型，另一方面，基于此算法進(jìn)行靈敏度指標(biāo)的計算，可以顯著降低模型計算樣本數(shù)量的需求，對模型參數(shù)較多的情況，效果更加明顯.本文模型有27個參數(shù)，多項式的最高階次為3，假定蒙特卡洛樣本數(shù)量10 000的情況下，蒙特卡洛算法，gPCE算法和本文提出的asPCE算法所需要的模型運(yùn)算樣本數(shù)量依次是290 000,4 060和2 582.

圖7 不同參數(shù)取值情況下模型輸出波形特征的相關(guān)性分析和Bland-Altman分析

5 結(jié) 語

當(dāng)臨床所用心血管模型參數(shù)很多時，LVPW估計的個性化和參數(shù)優(yōu)化計算工作將會很復(fù)雜.本文提出基于靈敏度分析的方法對模型中的參數(shù)進(jìn)行重要參數(shù)子集選擇，最終選出對LVPW有較大影響的參數(shù)子集，這樣的參數(shù)子集可為臨床患者的LVPW個性化估計提供參考，臨床中可以重點(diǎn)關(guān)注這些靈敏度指標(biāo)較高的參數(shù).參數(shù)子集中的參數(shù)數(shù)量減少，可以明顯降低模型參數(shù)優(yōu)化計算的復(fù)雜度.實(shí)驗(yàn)結(jié)果表明，基于靈敏度選擇的參數(shù)子集的LVPW估計，與模型全參數(shù)估計具有很高的一致性.