原俊青,駱 源

(浙江工業(yè)大學(xué) 理學(xué)院,浙江 杭州 310023)

顱內(nèi)動(dòng)脈瘤破裂會(huì)釋放大量血液,導(dǎo)致蛛網(wǎng)膜下腔出血,造成極高的發(fā)病率和死亡率[1-2]。然而治療顱內(nèi)動(dòng)脈瘤是一個(gè)復(fù)雜的過(guò)程,臨床治療對(duì)病人可能會(huì)產(chǎn)生負(fù)面影響,因此亟需建立一個(gè)可評(píng)估顱內(nèi)動(dòng)脈瘤破裂風(fēng)險(xiǎn)的數(shù)學(xué)模型。研究表明:顱內(nèi)動(dòng)脈瘤破裂與形態(tài)學(xué)參數(shù)有關(guān)[3-6],如動(dòng)脈瘤尺寸、瘤高和縱橫比等形態(tài)學(xué)參數(shù)。計(jì)算流體動(dòng)力學(xué)(CFD)研究討論了血流動(dòng)力學(xué)因素在腦動(dòng)脈瘤發(fā)病成形和破裂中的作用。壁剪應(yīng)力(WSS)、振蕩剪切指數(shù)(OSI)和相對(duì)滯留時(shí)間(RRT)被認(rèn)為是動(dòng)脈瘤破裂風(fēng)險(xiǎn)的指標(biāo)[7-10]。因此,有必要尋找一種聯(lián)合血流動(dòng)力學(xué)和形態(tài)學(xué)的方法以進(jìn)一步研究顱內(nèi)動(dòng)脈瘤的破裂風(fēng)險(xiǎn)可能性。

1 材料與方法

1.1 研究目標(biāo)

為了確定可評(píng)估動(dòng)脈瘤破裂風(fēng)險(xiǎn)的形態(tài)學(xué)參數(shù)和血流動(dòng)力學(xué)參數(shù),對(duì)190名患者的198例特異性側(cè)腦動(dòng)脈瘤模型進(jìn)行了CFD模擬,并獲得形態(tài)學(xué)特征。動(dòng)脈瘤位于基底動(dòng)脈、頸內(nèi)動(dòng)脈、大腦中動(dòng)脈、后交通動(dòng)脈和小腦后下動(dòng)脈等不同的腦區(qū)。所有患者動(dòng)脈瘤的幾何結(jié)構(gòu)都是通過(guò)三維旋轉(zhuǎn)血管造影圖像數(shù)據(jù)重建的,分別從醫(yī)學(xué)圖像和CFD模擬中計(jì)算8個(gè)形態(tài)學(xué)參數(shù)和6個(gè)血流動(dòng)力學(xué)參數(shù),并對(duì)參數(shù)進(jìn)行統(tǒng)計(jì)分析,評(píng)估每個(gè)參數(shù)的統(tǒng)計(jì)意義。

1.2 計(jì)算方法和條件設(shè)置

從Asenjo神經(jīng)外科研究所獲取一定數(shù)量的患者特異性側(cè)腦動(dòng)脈瘤圖像,并借助飛利浦Integris Allura三維旋轉(zhuǎn)血管儀顯示立體影像。根據(jù)Valencia等[11]所描述的程序,利用三維仿真模擬軟件Ansys Fluent重建96個(gè)破裂的動(dòng)脈瘤和102個(gè)未破裂的動(dòng)脈瘤。對(duì)不同病人的動(dòng)脈瘤模型進(jìn)行編號(hào),記作病例N。利用非結(jié)構(gòu)化四面體網(wǎng)格重建動(dòng)脈瘤模型,網(wǎng)格密度設(shè)置為每立方毫米1 500個(gè)單元,設(shè)置細(xì)胞的最大邊緣尺寸為0.2 mm。因此,含有最大單元數(shù)目為2 265 413的動(dòng)脈瘤模型的體積為1 521 mm3,含有最小單元數(shù)目322 215的動(dòng)脈瘤模型的體積為218 mm3。

使用PISO算法[12]進(jìn)行計(jì)算流體動(dòng)力學(xué)模擬。血液流速和壓力分別用二階迎風(fēng)插值和二階插值,并利用二階隱式公式離散時(shí)間,時(shí)間步長(zhǎng)設(shè)為Δt=0.000 5 s;動(dòng)量和連續(xù)性方程的殘差保持在10-3以下;實(shí)驗(yàn)中模擬的兩個(gè)心臟周期為0.85 s,并取第2個(gè)周期的結(jié)果;血液被認(rèn)為是不可壓縮的卡森流體[13-14],密度為1 065 kg/mm3,以層流形式流動(dòng)。

圖1 作為初始條件的血液流速示意圖Fig.1 Schematic diagram of blood flow rate used as the initial condition

通過(guò)以上方法可以在所有動(dòng)脈瘤的母血管入口處設(shè)置相同的血液流入平均速度,避免因血液流速初始條件不一致對(duì)實(shí)驗(yàn)造成影響。

模擬時(shí)對(duì)所有血管和動(dòng)脈瘤設(shè)置無(wú)滑移邊界條件,并在所有血管出口處設(shè)置在80～120 mmHg(1 mmHg=0.133 kPa)之間振蕩的壓力脈沖。壓力脈沖通過(guò)三元Windkessel電路法計(jì)算得出,結(jié)果如圖2所示。

圖2 作為邊界條件的壓力脈沖波形示意圖Fig.2 Schematic diagram of pressure waveform used as the boundary condition

1.3 形態(tài)學(xué)參數(shù)

從形態(tài)學(xué)角度進(jìn)行動(dòng)脈瘤的破裂風(fēng)險(xiǎn)分析,首先定義8個(gè)形態(tài)學(xué)參數(shù)來(lái)表征每個(gè)動(dòng)脈瘤:瓶頸因子(Bottleneck factor,BNF)、縱橫比(Aspect ratio,AR)、瘤體長(zhǎng)度/載瘤動(dòng)脈直徑(Size ratio,SR)、非球度系數(shù)(Nonsphericity index,NSI)、波動(dòng)系數(shù)(Undulation index,UI)、動(dòng)脈瘤角度(Aneurysm angle,αA)、流動(dòng)角度(Flow angle,αF)和血管角度(Vessel angle,αV)。

瓶頸因子BNF定義為瘤體寬度與頸寬的比值;縱橫比AR定義為動(dòng)脈瘤的瘤高與頸寬的比值;SR定義為動(dòng)脈瘤最大垂直高度與載流動(dòng)脈直徑的比值。

非球度系數(shù)NSI定義為

(1)

式中:V為動(dòng)脈瘤的體積;S為動(dòng)脈瘤的表面積。

波動(dòng)系數(shù)UI定義為

(2)

式中Vch為動(dòng)脈瘤的凸包體積,即完全包圍動(dòng)脈瘤的最小體積,其上所有的點(diǎn)都為凸集上的點(diǎn)[16]。

動(dòng)脈瘤角度定義為動(dòng)脈瘤頸部與動(dòng)脈瘤最大高度之間的夾角;流動(dòng)角度定義為動(dòng)脈瘤最大高度與母體血管之間的角度;血管角度定義為母血管與動(dòng)脈瘤頸部平面之間的角度。圖3更直接地顯示了動(dòng)脈瘤形態(tài)學(xué)參數(shù)的細(xì)節(jié)。

圖3 角度αA,αF,αV定義示意圖Fig.3 Definition of angles αA,αF,αV

1.4 血流動(dòng)力學(xué)參數(shù)

為了從血流動(dòng)力學(xué)角度進(jìn)行動(dòng)脈瘤的破裂風(fēng)險(xiǎn)分析,在每個(gè)動(dòng)脈瘤上計(jì)算6個(gè)血流動(dòng)力學(xué)參數(shù):舒張壁面剪應(yīng)力(DWSS)、收縮壁面剪應(yīng)力(SWSS)、時(shí)間平均壁剪應(yīng)力(TAWSS)、相對(duì)滯留時(shí)間(RRT)、振蕩剪切指數(shù)(OSI)和動(dòng)脈瘤形成指數(shù)(AFI)。

血液流動(dòng)會(huì)對(duì)血管內(nèi)壁產(chǎn)生剪切應(yīng)力WSS,記作τw,表示在心臟周期持續(xù)時(shí)間T的給定時(shí)間步長(zhǎng)上的WSS向量。收縮壁面剪應(yīng)力和舒張壁面剪應(yīng)力分別為在收縮時(shí)期的最大值和舒張時(shí)期的最小值的壁面剪應(yīng)力,其計(jì)算公式為

(3)

(4)

進(jìn)一步在動(dòng)脈瘤區(qū)域?qū)WSS和DWSS分別取平均值,并且分別利用母血管上的SWSS和DWSS指標(biāo)數(shù)值對(duì)動(dòng)脈瘤區(qū)域的SWSS平均值和DWSS平均值進(jìn)行歸一化。

時(shí)間平均壁剪應(yīng)力為WSS在一個(gè)心臟跳動(dòng)周期T上的平均值,其計(jì)算公式為

(5)

同樣在動(dòng)脈瘤區(qū)域上取得TAWSS平均值,并且利用母血管上的TAWSS指標(biāo)數(shù)值對(duì)其進(jìn)行歸一化。

振蕩剪切指數(shù)是一個(gè)無(wú)量綱參數(shù),計(jì)算在心臟周期內(nèi)WSS矢量方向的振蕩,將平均WSS向量與相同大小向量的均值進(jìn)行比較,其計(jì)算公式[17]為

(6)

在動(dòng)脈瘤區(qū)域取OSI平均值。

相對(duì)滯留時(shí)間為血液中的粒子顆粒在血管壁面附近的滯留時(shí)間,其計(jì)算公式[18]為

(7)

動(dòng)脈瘤形成指數(shù)量化了瞬時(shí)WSS向量和時(shí)間平均WSS向量之間的角度變化,即將瞬時(shí)WSS的局部方向與TAWSS的方向進(jìn)行比較,并忽略其大小,其計(jì)算公式[19]為

(8)

在整個(gè)動(dòng)脈瘤區(qū)域上取動(dòng)脈瘤形成指數(shù)平均值。

1.5 統(tǒng)計(jì)分析

通過(guò)Jarque-Bera檢驗(yàn)法,利用JB統(tǒng)計(jì)量確定各個(gè)參數(shù)是否符合正態(tài)分布,即

(9)

通過(guò)非參數(shù)Mann-WhitneyU檢驗(yàn)來(lái)評(píng)估動(dòng)脈瘤破裂組和未破裂組數(shù)據(jù)是否存在差異,將兩組樣本數(shù)據(jù)混合后進(jìn)行排列,并按1～m順序進(jìn)行編號(hào),計(jì)算兩組樣本的序號(hào)和(T1,T2),利用序號(hào)和(T1,T2)以及兩組樣本的數(shù)量來(lái)計(jì)算統(tǒng)計(jì)量U,即

(10)

式中:n1,n2分別為兩組樣本的病例數(shù)量。若統(tǒng)計(jì)量U大于臨界值Ua,兩組數(shù)據(jù)無(wú)差異;反之,兩組數(shù)據(jù)存在差異。

對(duì)于符合正態(tài)分布的參數(shù),應(yīng)用t檢驗(yàn)(Student’st-test)來(lái)評(píng)估正態(tài)分布的數(shù)據(jù),統(tǒng)計(jì)量t為

(11)

基于上述統(tǒng)計(jì)量,約定當(dāng)P<0.005時(shí),參數(shù)具有統(tǒng)計(jì)學(xué)上的顯著意義。此外,對(duì)各參數(shù)進(jìn)行接收機(jī)工作特性(ROC)分析,并獲得各參數(shù)的最佳閾值。通過(guò)計(jì)算ROC曲線下的面積(AUC)來(lái)衡量ROC曲線的有效性,面積越接近1,模型準(zhǔn)確性越高。

由于多元線性回歸不能很好地衡量動(dòng)脈瘤破裂風(fēng)險(xiǎn),故對(duì)所有具有統(tǒng)計(jì)學(xué)意義的參數(shù)采取多元logistic回歸,構(gòu)建最相關(guān)的參數(shù)與動(dòng)脈瘤的破裂風(fēng)險(xiǎn)大小之間的關(guān)系模型。將該關(guān)系模型描述為動(dòng)脈瘤破裂概率(p)與不破裂概率(1-p)的比值,稱為優(yōu)勢(shì)比,具體為

(12)

式中:p為動(dòng)脈瘤破裂風(fēng)險(xiǎn)大小;ai(i=1,2,…,j)為待估計(jì)的回歸參數(shù)。整個(gè)過(guò)程都使用Matlab軟件,通過(guò)極大似然估計(jì)法計(jì)算上述關(guān)系模型。

2 結(jié)果與分析

2.1 形態(tài)學(xué)和血流動(dòng)力學(xué)特征

動(dòng)脈瘤未破裂組和破裂組的參數(shù)均值和標(biāo)準(zhǔn)差見(jiàn)表1。從形態(tài)學(xué)角度來(lái)看,除了參數(shù)NSI和αA的差異較小,其余的參數(shù)在破裂組中均有較高的值。參數(shù)NSI,UI和αA的P值都明顯大于0.005,沒(méi)有統(tǒng)計(jì)學(xué)意義,對(duì)比以往研究成果[20]可知:該結(jié)果可能是由于顱內(nèi)動(dòng)脈瘤數(shù)據(jù)集較小所致。由表1可知研究結(jié)果中只有SR是有統(tǒng)計(jì)學(xué)意義的。至于參數(shù)AR,BNF,αF和αV,其P值較接近0.005,需要利用更多的動(dòng)脈瘤數(shù)據(jù)來(lái)評(píng)估其統(tǒng)計(jì)意義。

表1 顱內(nèi)動(dòng)脈瘤破裂組和未破裂組的形態(tài)學(xué)和血流動(dòng)力學(xué)參數(shù)Table 1 Morphological and hemodynamic parameters for the ruptured and unruptured groups

表1也顯示了破裂組和未破裂組的血流動(dòng)力學(xué)參數(shù)的統(tǒng)計(jì)分析結(jié)果。結(jié)果表明:動(dòng)脈瘤未破裂組的血流動(dòng)力學(xué)參數(shù)DWSS,SWSS和TAWSS均高于破裂組的對(duì)應(yīng)值,這與以往研究[9,21-22]結(jié)果相呼應(yīng)。此外,動(dòng)脈瘤破裂組DWSS和SWSS的比值和未破裂組DWSS和SWSS的比值較為相近,這意味著動(dòng)脈瘤腔與母血管內(nèi)的剪應(yīng)力類似,但是給出這一結(jié)論需要進(jìn)行更多醫(yī)學(xué)上的考慮?？傮w上未破裂組的WSS高于破裂組,而破裂組的RRT和OSI高于未破裂組,表明血液中的顆粒物質(zhì)傾向于在瘤壁面附近停留更長(zhǎng)的時(shí)間,并且在心臟周期中剪應(yīng)力變換方向更加頻繁。對(duì)于參數(shù)AFI,兩組的差異很小,因此難以利用其評(píng)估動(dòng)脈瘤破裂風(fēng)險(xiǎn)的概率。DWSS,SWSS和RRT具有統(tǒng)計(jì)學(xué)意義,與已有的研究結(jié)果[9,23]吻合較好。然而,其他血流動(dòng)力學(xué)參數(shù)也應(yīng)具有統(tǒng)計(jì)學(xué)意義[9,21],如OSI和TAWSS,這可能需要用到更多的病例數(shù)據(jù)進(jìn)行論證支撐[22]。用ROC分析計(jì)算了所有形態(tài)學(xué)和血流動(dòng)力學(xué)參數(shù)的最佳閾值。閾值和AUC見(jiàn)表2,由表2可知:一方面,參數(shù)NSI,UI和αA的AUC最小,而其P值在這些參數(shù)中相對(duì)較高;另一方面,SR的AUC最高(0.798),P值最低(0.002),意味著該參數(shù)是評(píng)估動(dòng)脈瘤破裂風(fēng)險(xiǎn)的良好指標(biāo),這一結(jié)論與過(guò)去的研究結(jié)果相符合,無(wú)論是三維結(jié)構(gòu)還是二維結(jié)構(gòu)的動(dòng)脈瘤,瘤體長(zhǎng)度與載瘤體動(dòng)脈的直徑比(SR)都是一個(gè)重要的評(píng)估因子[23-24]。AUC排名第2的參數(shù)是RRT(0.759),其次是DWSS和SWSS,兩者的AUC相近,分別為0.752和0.736。

表2 形態(tài)學(xué)和血流動(dòng)力學(xué)參數(shù)的閾值和AUC

2.2 多元統(tǒng)計(jì)分析

利用逐步回歸法在參數(shù)SR,DWSS,SWSS和RRT中選取最重要的變量。每引入一個(gè)參數(shù),就要對(duì)舊的參數(shù)逐個(gè)進(jìn)行檢驗(yàn),剔除偏回歸平方和不顯著的自變量,繼而進(jìn)行l(wèi)ogistic回歸確定動(dòng)脈瘤破裂風(fēng)險(xiǎn)評(píng)估的數(shù)學(xué)模型。分析過(guò)程中保留SR和DWSS作為相關(guān)參數(shù),軟件計(jì)算得到動(dòng)脈瘤破裂風(fēng)險(xiǎn)概率的生物數(shù)學(xué)模型,即

ODDall=e0.332SR-0.586DWSS+0.012

(12)

該模型中AUC為0.80。當(dāng)參數(shù)SR的數(shù)值增加一個(gè)單位量時(shí),動(dòng)脈瘤的破裂概率增加至原來(lái)的1.394倍;當(dāng)參數(shù)DWSS增加一個(gè)單位量時(shí),動(dòng)脈瘤的破裂概率降低至原來(lái)的0.557。因此,動(dòng)脈瘤破裂的概率對(duì)血流動(dòng)力學(xué)參數(shù)DWSS的變化更敏感,也說(shuō)明了瘤體長(zhǎng)度與載瘤體動(dòng)脈直徑的比值變大,會(huì)增大動(dòng)脈瘤破裂風(fēng)險(xiǎn)。將該模型的有效性與其他兩個(gè)模型的概率進(jìn)行比較:一個(gè)模型只包含參數(shù)SR,稱之為純形態(tài)模型;另一個(gè)模型只包含參數(shù)DWSS,稱之為純血流動(dòng)力學(xué)模型。對(duì)于第1種情況,純形態(tài)生物數(shù)學(xué)模型為

ODD形態(tài)學(xué)=e-1.542+0.387SR

(13)

該模型中AUC為0.79,當(dāng)參數(shù)SR的數(shù)值增加一個(gè)單位量時(shí),動(dòng)脈瘤破裂的概率增加至原來(lái)的1.459倍。對(duì)于第2種情況,純血流動(dòng)力學(xué)生物數(shù)學(xué)模型為

ODD動(dòng)力學(xué)=e1.243-0.634DWSS

(14)

該模型中AUC為0.75,當(dāng)參數(shù)DWSS的數(shù)值增加一個(gè)單位量時(shí),動(dòng)脈瘤破裂的概率降低至原來(lái)的0.530。比較上述3種數(shù)學(xué)模型,效率最高的模型是式(12)模型,即聯(lián)合模型,其與純形態(tài)學(xué)模型和純血流動(dòng)力學(xué)模型的差異分別為5.10%和7.58%。因此,這3種模型的ROC曲線非常相似,結(jié)果如圖4所示。

圖4 與聯(lián)合、純形態(tài)學(xué)純血流動(dòng)力學(xué)模型 相關(guān)的ROC曲線圖Fig.4 The ROC curves associated with the combined, pure morphological and pure hemodynamic models

采用聯(lián)合模型,即形態(tài)學(xué)-血流動(dòng)力學(xué)模型,可以更準(zhǔn)確、更全面地預(yù)測(cè)腦動(dòng)脈瘤破裂的風(fēng)險(xiǎn)。需注意的是,本研究中具有統(tǒng)計(jì)學(xué)意義的參數(shù)與其他研究中獲得的參數(shù)部分不同,例如Xiang等[9]表示SR,TAWSS和OSI是相關(guān)參數(shù),而Qin等[25]指出EL和D/W是最重要的參數(shù),其中D是動(dòng)脈瘤頸部到動(dòng)脈瘤圓尖端最長(zhǎng)尺寸之間的比率,W是動(dòng)脈瘤的頂部寬度,EL是能量損失。然而,Qin等[25]只考慮位于大腦中動(dòng)脈的動(dòng)脈瘤,而Xiang等[9]研究了位于不同腦區(qū)域的動(dòng)脈瘤?？傊?現(xiàn)有研究并未得到一致的答案,將來(lái)仍需針對(duì)評(píng)估動(dòng)脈瘤破裂風(fēng)險(xiǎn)的課題進(jìn)行更深入的研究。

3 結(jié) 論

從形態(tài)學(xué)和血流動(dòng)力學(xué)角度分析顱內(nèi)動(dòng)脈瘤數(shù)據(jù),對(duì)8個(gè)形態(tài)學(xué)參數(shù)和6個(gè)血流動(dòng)力學(xué)參數(shù)進(jìn)行多元統(tǒng)計(jì)分析,確定了用于評(píng)估動(dòng)脈瘤破裂風(fēng)險(xiǎn)的重要參數(shù)。實(shí)驗(yàn)得出如下結(jié)論:1) 瘤體長(zhǎng)度/載瘤動(dòng)脈直徑(SR)和舒張壁面剪應(yīng)力(DWSS)是對(duì)動(dòng)脈瘤破裂風(fēng)險(xiǎn)影響最顯著的參數(shù),可用于評(píng)估顱內(nèi)動(dòng)脈瘤破裂的風(fēng)險(xiǎn),利用這兩個(gè)參數(shù)分別得到的形態(tài)學(xué)-血流動(dòng)力學(xué)聯(lián)合模型(AUC為0.80)、純形態(tài)模型(AUC為0.79)和純血流動(dòng)力學(xué)模型(AUC為0.75)可以判斷動(dòng)脈瘤破裂風(fēng)險(xiǎn)的大小,并將風(fēng)險(xiǎn)量化;2) 3種模型有效性相似,聯(lián)合模型評(píng)估動(dòng)脈瘤破裂風(fēng)險(xiǎn)的效率更高,聯(lián)合模型顯示了瘤體長(zhǎng)度與載瘤動(dòng)脈直徑的比值增加一單位,動(dòng)脈瘤的破裂概率提升至原來(lái)的1.394倍,舒張壁面剪應(yīng)力增加一個(gè)單位,動(dòng)脈瘤破裂風(fēng)險(xiǎn)降低至原來(lái)的0.557,說(shuō)明動(dòng)脈瘤的形態(tài)尺寸以及血液流速的變化會(huì)使顱內(nèi)動(dòng)脈瘤的破裂風(fēng)險(xiǎn)以指數(shù)增長(zhǎng)的形式變化。考慮到CFD模擬相對(duì)較快,在臨床治療中為了提高準(zhǔn)確率,可以多關(guān)注參數(shù)SR和DWSS的數(shù)值變化,利用顱內(nèi)動(dòng)脈瘤破裂風(fēng)險(xiǎn)的聯(lián)合模型輔助臨床醫(yī)生進(jìn)行判斷。