基于自適應(yīng)細化網(wǎng)格的初期血栓形成計算模型

邊朔，馬朝青,2

1. 煙臺大學(xué) 計算機與控制工程學(xué)院，山東 煙臺 264005

2. 高端海洋工程裝備智能技術(shù)煙臺市重點實驗室，山東 煙臺 264005

血栓是活體血管系統(tǒng)中凝血物質(zhì)和血細胞形成的凝塊，會造成組織缺血性壞死，從而影響器官功能，甚至危及生命。數(shù)據(jù)表明，急性心肌梗塞、缺血性中風(fēng)、肺栓塞和彌散性血管內(nèi)凝血等血栓性疾病導(dǎo)致的死亡人數(shù)逐年增高[1]。血栓形成是臨床各科疾病中常見的一種病理過程，可見于心房室腔、動脈、靜脈和微血管。血栓的形成機制已通過多年的研究[2-4]獲得了基本認識。19 世紀中期，Virchow 提出了血栓形成病理生理機制的三要素：血流緩慢、血管破損和血液高凝狀態(tài)[5]。在20 世紀80 年代之前，針對血栓的研究主要依靠在活體內(nèi)和活體外進行的動物實驗以及實驗室儀器實驗。通過研究者的實驗確定了導(dǎo)致血栓形成的大多數(shù)復(fù)雜生物反應(yīng)級聯(lián)[6-8]。隨著對血栓形成機制研究的不斷地深入，一系列描述凝血過程及其子過程的數(shù)學(xué)模型被提出。隨著計算機技術(shù)的迅速發(fā)展，國內(nèi)外研究者對血栓形成機制的建模更傾向使用計算機模擬實驗研究血栓形成的機理，這一類研究的核心是建立血栓形成的計算模型。

國內(nèi)外研究者根據(jù)模型中血栓形成物質(zhì)的描述方式將計算模型分為3 類：連續(xù)模型，離散粒子模型和混合模型。連續(xù)模型[9-13]中血管、血流、血細胞以及凝血因子等物質(zhì)均描述為連續(xù)介質(zhì)，模型以偏微分方程的組合為主。與連續(xù)模型相反，離散粒子模型[14-20]將所有物質(zhì)描述為細微粒子，模型以粒子移動和相互作用為主?；旌夏Ｐ蚚21-23]中的凝血物質(zhì)和血小板等血細胞描述為離散粒子，血液和血管描述為連續(xù)物質(zhì)。因此，相比單純的連續(xù)模型和離散模型，在計算復(fù)雜度和跟蹤粒子行為方面具有更好的性能[24-25]。為了能夠描述較為復(fù)雜的血栓形態(tài)和獲得更平滑的血栓表面，與之前混合模型中用粒子簇對血栓進行建模不同，Ma 等[26]的研究中將血栓描述為具有可移動自由表面的連續(xù)體，并基于水平集方法建立模型。

雖然水平集方法對追蹤移動過程中拓撲結(jié)構(gòu)發(fā)生改變的輪廓非常有效，但對模型空間高一維的描述和偏微分方程的求解增加了模型的計算復(fù)雜度。文獻[27]提出一種用于研究動態(tài)變化過程中界面演化的算法，該算法基于水平集方法，使用窄帶和快速行進的方法應(yīng)用于動態(tài)自適應(yīng)網(wǎng)格，并結(jié)合了四叉樹等數(shù)據(jù)結(jié)構(gòu)。這種方法將主要的計算精力集中放在重要的區(qū)域上，節(jié)約大量的計算資源。文獻[28]提出了一種基于三維自適應(yīng)網(wǎng)格生成方法的局部細化技術(shù)，根據(jù)幾何特征有效地進行局部細化，網(wǎng)格細化質(zhì)量和靈活性得到顯著提高。

本文對血栓混合模型的建模網(wǎng)格進行了優(yōu)化，將窄帶法與局部細化相結(jié)合，提出針對可移動界面的自適應(yīng)細化網(wǎng)格，在保證計算準確率的基礎(chǔ)上減少數(shù)據(jù)量。同時，在血栓生長模型中引入正則化水平集演化（distance regularized level set evolution，DRLSE）方法[29]，有效降低水平集方程求解過程中的計算量。利用該模型對以血小板為主的血栓的形成過程進行更高效的仿真，分析血流和血管形狀對初期血栓形成的影響。

1 水平集方法中的自適應(yīng)網(wǎng)格細化

水平集是計算流體力學(xué)中常用到的一個重要概念。1988 年，Oshe 等[30]首次提出了水平集方法，有效解決了曲面演化問題的數(shù)值方法，并適用于任意維數(shù)空間。水平集方法的主要思想是將n維曲面的演化問題在n+1 維空間中解決，即將界面嵌入到一個超曲面中，并求解表示曲面演化的隱式水平集函數(shù)方程式。水平集方程主要包括3 個要素：超曲面的數(shù)據(jù)表示、控制曲面演化的一系列偏微分方程以及相應(yīng)的數(shù)值解法。

首先使用符號距離公式（signed distance function,SDF）建立表示超曲面的水平集函數(shù)：

式中：d是點x到曲面最近點的歐式距離，其符號取決于點x位于曲面的內(nèi)部還是外部，通常情況下定義內(nèi)部為負號，外部為正號。

水平集函數(shù)的演化遵守式(1)：

式中：F是曲線上各點的演化速度，方向為曲線的法線方向。已知初始水平集 φt=0，通過式(1)可以求得水平集函數(shù)在速度F的作用下經(jīng)過一個時間步長 Δt后的取值。

在迭代求水平集函數(shù)數(shù)值解的過程中涉及的偏微分方程求解，不論是使用低精度的迎風(fēng)（upwind）差分法，還是高精度的龍格庫塔（Runge-Kutta）方法，都會產(chǎn)生較大的計算量，限制水平集方法在各領(lǐng)域的應(yīng)用。

在水平集方法的各應(yīng)用領(lǐng)域中，速度F通常僅定義在零水平集上。但水平集函數(shù)是在整個計算空間上進行更新的，因此需要將F也擴展到整個空間。這使得原始水平集方法在求解大遲鈍輪廓時計算量大且效率低下?？紤]到在某一時間點上水平集函數(shù)的一系列取值中零水平集只有一個，并且被關(guān)注的也只有零水平集所代表的曲面，因此可以采用窄帶算法[31]僅求取零水平附近的函數(shù)值來降低計算量。窄帶方法首先在輪廓的內(nèi)部和外部分別設(shè)定寬度為 θ的帶狀區(qū)域，形成一個帶寬為 2θ的窄帶圍繞零水平集，如圖1(a)。在對水平集函數(shù)進行求解時，只需要將速度方程擴張至窄帶邊界，并且僅對窄帶內(nèi)部的網(wǎng)格進行水平集更新。

利用水平集方法追蹤運動界面時，界面精度和質(zhì)量損失都與網(wǎng)格的分辨率有關(guān)。如果所有網(wǎng)格都進行細化，則計算時間和存儲空間將大大增加。因此，本文采用自適應(yīng)網(wǎng)格技術(shù)對窄帶區(qū)域網(wǎng)格進行細化，對遠離界面區(qū)域的網(wǎng)格進行粗化處理，以確保界面附近的網(wǎng)格具有更高的分辨率，提高網(wǎng)格的計算精度。假設(shè)窄帶的帶寬為 2θ，如果網(wǎng)格上水平集的值滿足 |φ|≤θ，則需要進行細化（圖1(b)）。

圖1 窄帶示意

在窄帶法中，窄帶內(nèi)的網(wǎng)格被存儲壓縮為一維數(shù)組。但水平集在演化過程中，不僅窄帶會定期更新，窄帶內(nèi)網(wǎng)格中水平集數(shù)據(jù)也頻繁地被調(diào)取和更新。細化后的笛卡爾網(wǎng)格在數(shù)據(jù)存儲和查找上更為復(fù)雜，因此不能繼續(xù)使用原始的壓縮矩陣形式存儲細化數(shù)據(jù)。針對三維笛卡爾細化網(wǎng)格，八叉樹結(jié)構(gòu)更容易管理網(wǎng)格數(shù)據(jù)，提高計算和存儲效率。

2 初期血栓計算模型

本文針對初期血栓形成過程建立particle-continuum 混合模型。初期血栓常見于動脈或血流速度較快的靜脈，是在血栓形成初期由大量血小板和少量纖維蛋白組成的小節(jié)狀或贅生物狀凝塊，也稱為血小板血栓或白色血栓。因此，本模型包含的主要元素有擬設(shè)有破損內(nèi)壁的血管、血管性血友病因子（von Willebrand Factor，vWF）、血漿、血小板和初期血栓。在混合模型中，血管、vWF、血漿和血栓屬于連續(xù)介質(zhì)（continuum），為了能夠觀察血小板在血栓形成過程中的運動行為，利用離散顆粒（particle）模擬血小板。圖2 是模型的三維示意圖，模型中的血管被定義為單入口單出口的通道，血漿充滿整條血管并從入口流向出口。在血管中段預(yù)設(shè)小區(qū)域的血管內(nèi)膜破損和血管狹窄。破損后暴露出的皮下膠原纖維引發(fā)凝血機制，由血漿從血管入口帶入的血小板會在vFW 的作用下黏附在破損內(nèi)膜處，并引發(fā)更多的血小板聚集，逐步生長為初期血栓。

圖2 帶有破損血管壁的初期血栓生長三維模型示意

根據(jù)模型中包含的元素，計算模型被分為4 個子模塊：血漿模型、vWF 模型、血小板模型和血栓模型。血漿模塊對血管中的血漿流動進行建模，該模塊為整個模型提供血液動力學(xué)環(huán)境。vWF 模塊對vWF 因子的對血小板的影響進行判斷，其釋放和擴散情況受血漿模塊提供的速度場影響，該模塊中的vWF 因子將誘導(dǎo)血小板激活。血小板模塊對血小板的移動進行建模，血小板的移動受血漿的速度場和vWF 的共同影響。血栓模塊對血栓的生長進行建模，組成血栓的血小板由血小板模塊提供。

2.1 血漿模塊

在該模塊中，本文假設(shè)血漿為不可壓縮的粘性牛頓液體。將引入Navier-Stokes（N-S）方程[17]作為血漿的速度函數(shù)：

通過對N-S 方程的分步求解可獲得血管內(nèi)的血漿速度場，并將其同時用于vWF 模塊、血小板模塊和血栓模塊。

2.2 vWF 模塊

vWF 是參與初期血栓形成的主要凝血因子之一，主要由內(nèi)皮細胞合成分泌后存在于血漿中，內(nèi)皮細胞在受刺激或損傷后，血漿內(nèi)局部vWF 濃度會提高，導(dǎo)致血栓形成。vWF 在初期血栓形成中主要參與血小板的黏附和聚集過程。首先，vWF與血小板膜糖蛋白GPIⅠb-Ⅸ復(fù)合物結(jié)合，引導(dǎo)血小板降低移動速度并黏附于血管壁損傷部位形成血小板覆蓋；然后與GPⅡb-Ⅲa 結(jié)合，誘導(dǎo)血小板聚集，進一步擴大血小板血栓的體積。vWF 對血小板的作用受血液剪切速率影響，其與GPIⅠb-Ⅸ的結(jié)合在剪切速率大于1 000-1時發(fā)生、與GPⅡb-Ⅲa的作用在剪切速率1 000-1～10 000-1發(fā)生[31]。

在本模塊中，根據(jù)局部的剪切率確定vWF 因子是否對血小板產(chǎn)生作用以及產(chǎn)生何種作用。使用血漿模塊計算得出的血漿流速u可計算血管內(nèi)血流的剪切速率：

當血小板所在位置剪切速率大于閾值γact=1 000-1時，血小板被激活，繼而發(fā)生黏附和聚集反應(yīng)；當剪切速率大于閾值γno-ag=10 000-1時，僅產(chǎn)生黏附作用。由于vWF 因子的濃度在破損內(nèi)膜處最高，設(shè)定血小板所受vWF 的影響程度隨血小板與破損內(nèi)膜距離的增大線性遞減。

2.3 血小板模塊

根據(jù)初期血栓形成過程中血小板在vWF 影響下產(chǎn)生的黏附和聚集作用，定義2 種血小板狀態(tài)：靜息狀態(tài)和激活狀態(tài)；以及3 種作用于血小板上的力：血漿拖拽力、黏附誘導(dǎo)力和聚集誘導(dǎo)力。

正常狀態(tài)下，血小板為靜息狀態(tài)，僅在血漿拖拽力作用下隨血流移動，血漿拖拽力計算方式為

式中：Fad_max、Fag_max分別為黏附力和聚集力的最大值；Dad、Dag分別為血小板中心與破損內(nèi)膜和血栓表面的最近距離；nad、nag分別為從血小板指向破損內(nèi)膜和血栓表面最近點的法向量，用于確定受力方向?？紤]vWF 濃度隨與破損區(qū)域的距離增大而減小，誘導(dǎo)力的大小和血小板與破損內(nèi)膜或血栓表面的距離成反比。

2.4 血栓模塊

血小板在黏附和聚集的過程會被激活并融合到血栓中。由于模型中血小板屬于顆粒，血栓屬于連續(xù)介質(zhì)，融合過程是從離散顆粒到連續(xù)介質(zhì)的轉(zhuǎn)化的過程。血栓模塊主要利用水平集方法，完成這一過程中的轉(zhuǎn)化，實現(xiàn)血栓的生長。同時，為了提高水平集方法的計算效率，水平集的定義及水平集方程的求解引入了自適應(yīng)細化網(wǎng)格方法。

首先，與血液、血管壁和血栓均使用基于符號距離公式的水平集函數(shù) φ描述:

發(fā)生黏附或聚集的血小板融入血栓時，血栓的表面將向血漿一側(cè)擴張。本文定義血栓表面的變化分為2 個步驟，首先按血小板形狀擴張表面，然后用小幅度的曲率速度平滑血小板與原表面的連接邊緣。為了在擴張過程使水平集保持符號距離公式的形式，本文引入DRLSE。在第一步的擴張中，使用該方法代替Fast Matching Method 高效地初始化水平集；第二步的平滑化中加入正則項，在水平集方程求解過程中省略多次再初始化步驟。DRLSE 的水平集演化方程為

3 仿真實驗結(jié)果與分析

3.1 狹窄血管內(nèi)的仿真實驗

血管狹窄是導(dǎo)致血栓形成的重要原因之一，本文應(yīng)用提出的初期血栓形成模型中建立具有局部狹窄的動脈血管3D 模型，對初期動脈血栓的形成過程進行了仿真，并分別對狹窄程度和破損位置對血栓的生長速度和幾何形態(tài)的影響進行分析。

仿真實驗中的血管設(shè)定是直徑為 4 mm的頸動脈，血流速度為1 m/s。實驗分別在管狀直血管和有局部狹窄的血管中進行，血管狹窄設(shè)定為動物實驗中通常使用的針尖壓迫血管產(chǎn)生的局部狹窄。為了模擬這一過程，在對血管進行3D 形態(tài)建模時，利用二維高斯方程獲得血管壁的峰狀突起，通過調(diào)整高斯方程的系數(shù)控制突起的高度和范圍。圖3 為實驗中使用的2 種狹窄血管，最狹窄處分別占血管半徑25%和50%，黃色為血管突起區(qū)域。

圖3 仿針使用的2 種狹窄血管

首先，設(shè)定全部突起部分存在血管內(nèi)壁破損，圖4 為狹窄程度為25%時血栓形成過程仿真圖，黃色區(qū)域為血管破損區(qū)域及破損上形成的初期血栓。圖5(a)為血管無狹窄和25%狹窄時破損出現(xiàn)0.3 s 后血栓內(nèi)血小板的數(shù)量。結(jié)合圖4可見，無狹窄情況下狹窄出口和入口處血栓生長程度相似，狹窄使生長速度上升且出口處速度高于入口出速度。E. Westein 等[33]通過動物實驗證明了血管出現(xiàn)狹窄后的血流狀態(tài)改變會導(dǎo)致受vWF 主導(dǎo)的血栓生長速度加快，并且血小板在凸起后半段的聚集速度大于前半段，如圖5(b)所示。圖5 中，本文仿真結(jié)果（圖5(a)）與該實驗結(jié)果（圖5(b)）一致，可驗證本文模型的有效性。

圖4 25%狹窄下的血栓形成，黃色區(qū)域中突起為血栓

圖5 無狹窄和有狹窄下同一時間內(nèi)血栓內(nèi)包含血小板數(shù)量

圖6 為血管狹窄程度為25%的情況下，破損出現(xiàn)在不同位置時初期血栓的形成過程。當破損區(qū)域在狹窄上游時，初期血栓中血小板多聚集于血流入口處，因為該區(qū)域具有率先影響和捕獲血小板的位置優(yōu)勢。當破損區(qū)域位于狹窄入口和出口時，狹窄促使頂端處出現(xiàn)明顯的聚集。但是當破損出現(xiàn)在下游時，血小板反而在血流出口方向聚集。

圖6 不同血管內(nèi)壁破損位置的血栓形成過程

本文在25%和50%狹窄的血管中，分別對破損區(qū)域出現(xiàn)在狹窄上游、入口、出口和下游這4 種條件下進行了10 次血栓形成仿真實驗，獲得了血栓在0.5 s 內(nèi)的生長情況并繪制圖7。如圖7所示，血管狹窄程度增大會促進血栓形成，且狹窄出口處的形成速度大于狹窄入口處，這與之前實驗結(jié)果（圖5(a)）和文獻[33]中的實驗結(jié)果（圖5(b)）一致。結(jié)合前文分析的血小板聚集區(qū)域，由于當破損在出口時，更多的血小板聚集在狹窄的頂峰區(qū)域，該情況下更容易造成血管的堵塞。

圖7 0.5 s 時不同狹窄程度下血栓在不同破損區(qū)域的生長情況

3.2 仿真實驗效率比較與分析

在仿真使用模型的血栓模塊中，結(jié)合了窄帶法的自適應(yīng)細化網(wǎng)格被使用在水平集方程求解上，以降低需要存儲和計算的數(shù)據(jù)量。圖8 是在網(wǎng)格間距為血小板半徑和2.5 倍血小板半徑下進行仿真實驗的血栓形態(tài)。按血栓模塊設(shè)計的DRLSE 方法，血栓能夠較為完整地保留血小板的球狀形態(tài)。使用粗網(wǎng)格的實驗中，由于水平集演化方程求解中的誤差，血栓中血小板形態(tài)出現(xiàn)失真和變形；在網(wǎng)格細化提高計算精確度后，血小板的形態(tài)得到了更好地保留。

圖8 不同網(wǎng)格間距下血栓的形態(tài)

在水平集方法中，窄帶法雖然可以降低算法的計算量，但均勻的局部細化仍然會導(dǎo)致數(shù)據(jù)的冗余。圖9 是無狹窄、狹窄度為25% 和50%這3 種情況下，局部均勻細化和非均勻的自適應(yīng)細化下網(wǎng)格數(shù)量隨血栓生長的曲線圖。雖然在血栓生長至一定程度后細化網(wǎng)格增長不再明顯。但隨狹窄程度增加，均勻細化下網(wǎng)格數(shù)量成倍增長。使用自適應(yīng)局部細化后，網(wǎng)格數(shù)量下降至均勻細化的2%，并且狹窄程度和血栓生長對細化網(wǎng)格數(shù)量無明顯影響。

圖9 仿真過程中不同狹窄程度下細化網(wǎng)格數(shù)變化曲線

4 結(jié)論

本研究中針對初期血栓形成機制的研究，建立了基于水平集的計算模型，并采用自適應(yīng)的局部細化網(wǎng)格和DRLSE 方法優(yōu)化了水平集函數(shù)求解過程。通過研究，得出結(jié)論如下：

1)本文提出的初期血栓形成計算模型能夠被用于實現(xiàn)狹窄血管中初期血栓形成的仿真實驗，實驗結(jié)果展示的狹窄對血栓形成的影響情況與動物實驗一致。

2)通過實驗分析可知，血管狹窄會促進血栓的形成，血小板在狹窄出口處的聚集較為明顯，血管堵塞的風(fēng)險較大。

3)結(jié)合了局部水平集和笛卡爾自適應(yīng)網(wǎng)格的自適應(yīng)局部細化網(wǎng)格可以在提高計算精確度的基礎(chǔ)上有效降低模型計算中水平集函數(shù)求解的數(shù)據(jù)量。

4)文中所提出的計算模型雖然引入了vWF在血栓形成中的作用，但仍有更多的凝血因子可以加入模型以提高模擬過程的完整性，這也是下一步研究工作的重點內(nèi)容。