邊朔，馬朝青,2

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖1 窄帶示意

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

2 初期血栓計(jì)算模型

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

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

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

2.1 血漿模塊

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

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

2.2 vWF 模塊

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

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

當(dāng)血小板所在位置剪切速率大于閾值γact=1 000-1時(shí)，血小板被激活，繼而發(fā)生黏附和聚集反應(yīng)；當(dāng)剪切速率大于閾值γno-ag=10 000-1時(shí)，僅產(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)，僅在血漿拖拽力作用下隨血流移動，血漿拖拽力計(jì)算方式為

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

2.4 血栓模塊

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

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

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

3 仿真實(shí)驗(yàn)結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

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

3.2 仿真實(shí)驗(yàn)效率比較與分析

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

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

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

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

4 結(jié)論

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

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

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

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

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