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

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

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

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

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

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

1 體循環(huán)模型

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

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

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

圖2 左心室壓力波形特征

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

表1 體循環(huán)模型中的生理類(lèi)參數(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ù)靈敏度計(jì)算

2.1 靈敏度指標(biāo)

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

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

Y=f(X) .

(1)

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

(2)

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

(3)

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

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

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

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

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

(4)

(5)

(6)

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

(7)

(8)

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

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

(9)

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

(10)

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

(11)

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

(12)

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

(13)

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

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

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

(14)

(15)

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

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

由于asPCE算法中，求解多項(xiàng)式系數(shù)需要用到線(xiàn)性回歸方法，因此需要提前進(jìn)行模型樣本的運(yùn)算.在通用的多項(xiàng)式混沌展開(kāi)(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ù)來(lái)形成更精細(xì)的單位間隔，并且可以生成輸入?yún)?shù)的基于均勻分布的隨機(jī)序列.

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

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

圖3 asPCE算法流程與靈敏度計(jì)算

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

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

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

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

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

4 模型輸出LVPW對(duì)比及靈敏度計(jì)算方法對(duì)比

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

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

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

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

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

5 結(jié) 語(yǔ)

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