基于晶格玻爾茲曼方法的多段淋巴管模型*

張乾毅 韋華健 李華兵

(桂林電子科技大學(xué)材料科學(xué)與工程學(xué)院,桂林 541004)

淋巴系統(tǒng)對人體的免疫及細胞的內(nèi)環(huán)境穩(wěn)態(tài)都有著重要的作用.與血液循環(huán)系統(tǒng)相似,淋巴系統(tǒng)也是遍布全身的管道系統(tǒng),主要由淋巴液、淋巴管和淋巴器官構(gòu)成.淋巴管的自發(fā)收縮驅(qū)動管內(nèi)淋巴液的流動.淋巴管的自發(fā)收縮-舒張機制由Ca2+和NO 濃度的振蕩反饋決定,NO 在管內(nèi)的分布對淋巴管的收縮循環(huán)起到重要作用.因淋巴液流動而作用在淋巴瓣膜上的剪切力是瓣膜產(chǎn)生NO 的主要原因.在真實系統(tǒng)中,某段淋巴管中的NO 分布會受到與其連接的淋巴管的影響,特別是上游的連接片段.通過晶格玻爾茲曼方法,建立了1 個具有瓣膜結(jié)構(gòu)的多段淋巴管模型,再現(xiàn)了淋巴管內(nèi)Ca2+和NO 的反饋機制,瓣膜變化和淋巴液流動情況.該模型中存在3 種淋巴管,分別是初始淋巴管、中間淋巴管和出口淋巴管.淋巴管的段數(shù)可以通過修改計算參數(shù)無限擴充.本文計算的段數(shù)為3—5 段,每段淋巴管中有兩對瓣膜.通過模型研究了多段淋巴管中NO 濃度分布、壓力分布、NO 平均濃度變化,以及3 段管模型中各管的流量隨時間變化情況.

1 引言

淋巴瓣膜出現(xiàn)在集合淋巴管中,是毛細淋巴管向集合淋巴管過渡的標(biāo)志,淋巴瓣膜將集合淋巴管分成了若干小段,起到防止淋巴液回流的作用[1].與血液循環(huán)系統(tǒng)中,將心臟作為單一泵源并通過沒有激勵源的血管進行循環(huán)的方式不同,淋巴管表現(xiàn)為一組串聯(lián)排列的分布式泵[2],每段被瓣膜分隔開的獨立淋巴管段都具有自發(fā)收縮驅(qū)動淋巴液流動的能力[3-6].集合淋巴管可以根據(jù)外部壓力條件,自行調(diào)整其收縮行為以維持體液穩(wěn)態(tài),若存在外部的流體壓力梯度可以驅(qū)動管內(nèi)淋巴液流動,集合淋巴管可最小化收縮,否則便主動收縮驅(qū)動淋巴液流動[7].

淋巴管的自主收縮機制為,當(dāng)淋巴管中Ca2+濃度超過淋巴肌細胞的閾值濃度時,淋巴肌細胞中會形成蛋白質(zhì)間的交叉連接,從而引起自主收縮.淋巴管的收縮導(dǎo)致淋巴液流動并產(chǎn)生剪切力,淋巴管內(nèi)皮細胞響應(yīng)動態(tài)流體剪切力并激活其中的一氧化氮合成酶(eNOS)產(chǎn)生并釋放一氧化氮(NO)[8],引起淋巴管舒張.由此可見,是NO 和Ca2+濃度的周期性變化引起了淋巴管的收縮-舒張過程.

相較于動脈有二十多個模型可供參考,淋巴系統(tǒng)的數(shù)值模擬還處在初級階段.最早的淋巴系統(tǒng)網(wǎng)絡(luò)模擬是由Reddy 和Black[9]根據(jù)Navier-Stokes方程設(shè)計的簡化一維模型.后來,Bertram 等[10,11]用集總方法模擬了小段淋巴管,發(fā)現(xiàn)瓣膜阻力隨著施加在瓣膜上的壓力梯度的變化而變化,但其瓣膜是虛擬的,只起到單向作用.研究表明,瓣膜是產(chǎn)生一氧化氮的重要來源.在前期研究中,建立了存在真實瓣膜的集合淋巴管模型,但其存在著明顯局限性,只有1 個單獨的淋巴管段.而在真實的淋巴系統(tǒng)中,淋巴管段之間首尾連接,某段淋巴管中的NO 的分布會受到其上下游相鄰管段的影響,本文在之前模型的基礎(chǔ)上,建立了1 個可以調(diào)整淋巴管段數(shù)目的多段淋巴管模型.

2 晶格玻爾茲曼方法

晶格玻爾茲曼方法(lattice Boltzmann method,LBM)是1 種介觀尺度的模擬方法,在流體力學(xué)方面的模擬中被廣泛應(yīng)用,如多孔介質(zhì)流[12]、血液流[13]、多相流[14]等.LBM 以微觀的流體粒子分布函數(shù)為基礎(chǔ),用分布函數(shù)對速度集求各階矩得到相關(guān)的宏觀量,避開了傳統(tǒng)數(shù)值方法直接求解宏觀變量所帶來的復(fù)雜非線性項,大大簡化了計算的難度.在流體碰撞和流動的過程中,粒子總體滿足質(zhì)量和動量守恒,粒子之間的相互作用由局部分布函數(shù)決定.

Mcnamara 和Zanetti[15]用單粒子分布函數(shù)fi取代了格子氣細胞自動機中的布爾變量,其LBM演化方程為

其中i表示微觀速度的索引,對于D2Q9 模型(見圖1),i=0,1,2,···,8;fi(x,t) 是在x位置、t時刻,具有ei速度的粒子的分布函數(shù),δt為時間步長,δx為格子長度;Ω(fi)是碰撞因子,表示碰撞對fi的影響.Chen 等[16]和Qian 等[17]提出用單弛豫時間來替代碰撞因子項.繼而,LBM 方程可以寫為

圖1 D2Q9 模型Fig.1.D2Q9 model.

其中,ρ是宏觀密度,u是格點上流體的流速.本文采用LBM 中的D2Q9 模型[18]進行計算,其平衡分布函數(shù)具體定義如下式:

式中,c為基準(zhǔn)速度,通常情況下取1.ωi為各方向的權(quán)重系數(shù),當(dāng)i=0 時,ωi=4/9 ;當(dāng)i=1,2,3,4 時,ωi=1/9;當(dāng)i=5,6,7,8 時,ωi=1/36.ei的形式為

邊界條件的選擇對于LBM 的計算效率、計算精度和穩(wěn)定性都起到非常重要的作用.因為在迭代過程中,經(jīng)過格子的碰撞和流動遷徙過程之后,流場內(nèi)格點的分布函數(shù)已經(jīng)更新,而邊界上仍存在未更新的格點分布函數(shù).因而需要用確定的邊界條件,并采用對應(yīng)的邊界處理格式對未更新的格點分布函數(shù)進行確定,才能繼續(xù)進行下一步的迭代.

壓力邊界條件用來處理存在壓差條件下,流體的流動以及邊界滑移時流體的流動.如圖2 所示,在左入口處,經(jīng)過流動遷徙后指向流場內(nèi)部的3 個分布函數(shù)f1,f5和f8為未知量.相同的,在右出口處f3,f6和f7也為未知量.需要擬定好出入口密度值,并假設(shè)出入口處速度的垂直分量uy=0,并進一步求出出入口速度水平分量ux和未知的分布函數(shù).根據(jù)(3)式、(4)式可以得出以下關(guān)系[13]:

圖2 壓力邊界條件示意圖Fig.2.Schematic diagram of pressure boundary condition.

用LBM 對實際物理體系進行模擬時,必須對模型中的物理量進行量綱轉(zhuǎn)換,即實際物理量與離散格子量之間的轉(zhuǎn)換.在本文的計算中,用L表示長度量綱,T表示演化的時間量綱,G表示質(zhì)量量綱.在各量綱符號后加上“′”代表格子上的量綱.將顆粒的直徑D,流體密度ρ和黏滯系數(shù)v作為基準(zhǔn)物理量,則:.根據(jù)這些基準(zhǔn)物理量,可繼續(xù)求得其他物理量.

3 淋巴管模型

在真實的淋巴系統(tǒng)中,淋巴管段間彼此連接.某段淋巴管段的收縮情況及其內(nèi)部NO 和Ca2+的分布情況都受到其他連接的淋巴管段的影響,特別是其上游連接管段.建立的多段淋巴管模型包含有3 種淋巴管段,總段數(shù)可以通過修改MPI 計算參數(shù)進行調(diào)整.每種淋巴管段中都含有兩個腔內(nèi)瓣膜和淋巴管壁.如圖3(a)所示為初始淋巴管,其直徑為0.01 cm,虛線部分表示具有可供組織液自由流入的多孔入口.在多孔入口處設(shè)定恒定的反彈比,只允許初始淋巴管中15%的淋巴液滲透回組織.圖3(b)為中間淋巴管.圖3(c)為出口淋巴管,其尾端表示0.0296 cm 長的靜態(tài)淋巴管壁,用來降低出口處淋巴管壁的影響.淋巴管段的外部為組織,假設(shè)組織處不會流動,只有當(dāng)淋巴管膨脹,淋巴管壁的位置越過組織時,組織才會轉(zhuǎn)化成流體.當(dāng)組織轉(zhuǎn)化為流體時,新產(chǎn)生的流體將由其相鄰的流體計算得出.通過LBM 中的D2Q9 模型來模擬淋巴管內(nèi)外部的流體.淋巴管壁、淋巴瓣膜通過曲線邊界條件實現(xiàn),流體作用力通過壓力張量積分法計算.組織處、初始淋巴管的入口和出口淋巴管的靜態(tài)部分通過反向彈回邊界條件實現(xiàn).最外層的黑色加粗壁由平衡分布函數(shù)實現(xiàn).其速度為0,密度為1.為了保護出口淋巴管中尾端的瓣膜,在距離出口0.0016 cm處引入了1 個部分反向彈回邊界,流體可以自由流出,若外部流體有從出口處流入的趨勢,該邊界會阻止85%的流體流入.若尾端的瓣膜完全關(guān)閉,部分反彈變?yōu)榉聪驈椈?

圖3 淋巴管示意圖Fig.3.Schematic diagram of lymphatic vessel.

模型中淋巴管壁和瓣膜均被離散化,只能沿著豎直方向移動,二者均受彎曲力作用,可表示為

其中,ym和yn是格子相鄰兩點y軸上的位置.KB在淋巴管中是個常數(shù),而在瓣膜中用帶有上標(biāo)v(表示瓣膜valve)的表示瓣膜的抗彎強度,此處抗彎強度具體表示為

其中,i表示當(dāng)前段數(shù),n表示瓣膜的總段數(shù).是在固定點位置處的最大值,是最小值的近似值.系數(shù)A調(diào)節(jié)瓣膜的柔軟度,瓣膜在淋巴液回流時要有足夠的強度抵擋流體的沖擊,因此瓣膜的柔性系數(shù)由淋巴管的中間線附近的自由點處向淋巴管壁固定點處逐漸增大.

彈性力FE來自淋巴管外組織處:

式中,R表示淋巴管和瓣膜的半徑,R0在淋巴管中是個常數(shù).但在瓣膜中,為了保持瓣膜的拋物線形狀,假設(shè)瓣膜在初始位置的拋物線形狀表示為

式中,yl0表示瓣膜的初始位置,x0表示瓣膜的固定點水平方向位置,正負號分別表示上下瓣膜,B用于調(diào)整初始位置是偏向于閉合或者張開.淋巴管的最小半徑為Rl.為避免瓣膜張開程度過大,用一個最大B值確定的位置為其極限位置.

因淋巴管和瓣膜都具有黏彈性,黏滯阻力可表示為

負號表示Fr總是作用于壁面沿速度v的反方向.

當(dāng)兩個瓣膜位置接近時,因為缺少流體格點而不能使用LBM 來計算流體力.在這種情況下,使用潤滑理論來計算流體力[19],又因為瓣膜只在y方向上運動,只考慮垂直方向上的潤滑力,可表示為

式中,l和h分別表示瓣膜片段的長度和瓣膜片段間的距離,va是該段瓣膜的速度.

除了上述彎曲力、彈性力和黏滯阻力外,作用于淋巴管上的力還包括淋巴肌力和流體力.淋巴管由淋巴肌細胞組成,當(dāng)大量鈣離子擴散到淋巴肌細胞細胞質(zhì)時,會觸發(fā)淋巴肌細胞收縮.淋巴肌力與Ca2+和NO 的濃度相關(guān),可表示為

式中,CCa和CNO分別表示鈣離子和一氧化氮的濃度,KM是決定淋巴肌力強度的常量.Ca2+的反應(yīng)擴散方程表示為

式中,DCa是Ca2+的擴散系數(shù),是Ca2+的衰減速率,基于NO 通過肌球蛋白輕鏈磷酸酶作用于平滑肌細胞可以減少收縮力的產(chǎn)生,Ca2+的總衰減速率需乘上 (1+KCa,NOCNO)[5].是Ca2+的生成速率,11 次項是淋巴管拉伸產(chǎn)生Ca2+.δ ↑是非對稱Kroneckerδ函數(shù),當(dāng)CCa從低于閾值Cth增加到Cth時,該函數(shù)設(shè)定為1,其他情況設(shè)定為0;當(dāng)CCa達到閾值Cth時,它觸發(fā)由電壓門控和鈣誘導(dǎo)的鈣通道介導(dǎo)動作電位[20,21];表示當(dāng)CCa＞Cth時CCa的急劇增大.R和RCa分別表示淋巴管的局部半徑和非線性項參考半徑,Rl是淋巴管的極限最小半徑.一氧化氮的產(chǎn)生和擴散可以表示為

式中,CNO表示NO 的濃度,DNO表示NO 的擴散系數(shù),分別表示NO 的衰減速率和生成速率.

4 計算結(jié)果與分析

本次計算代碼使用C++編寫,計算代碼在搭載了8 個GPU 的服務(wù)器上執(zhí)行.GPU 型號為NVID IA Quadro GP100,每個GPU 含有3584 個CUDA核心,處理雙精度浮點數(shù)的能力為5.2 TFLOP/s,搭配16 GB HBM2 顯存,理論帶寬高達717 GB/s.編譯環(huán)境為CUDA10.0.每個GPU 用來計算一段淋巴管,淋巴管數(shù)分別為3,4,5 段.MPI 用來交換淋巴管段邊界間的數(shù)據(jù).通過量綱轉(zhuǎn)換將實際淋巴管物理參數(shù)與格子物理量進行轉(zhuǎn)換,其詳細參數(shù)參考文獻[22].

圖4 表示多段管中的NO 濃度分布和壓力分布情況,右側(cè)顏色標(biāo)尺自下而上的顏色變化表示管內(nèi)NO 濃度和壓力強度的增大.圖4(a)、圖4(c)和圖4(e)分別表示淋巴管的段數(shù)為3,4,5 段時NO 的濃度分布情況,瓣膜是產(chǎn)生NO 的主要結(jié)構(gòu),瓣膜附近的NO 濃度要明顯高于其他部分.隨著瓣膜數(shù)量的增加,NO 在整個管腔內(nèi)的分布更為復(fù)雜,Ca2+的分布變得越來越隨機,進而導(dǎo)致Ca2+的濃度達不到淋巴肌細胞收縮的閾值,淋巴管整體的收縮情況也變得混亂.此外,瓣膜附近NO 的濃度還受淋巴管當(dāng)前收縮階段的影響.當(dāng)淋巴管收縮時,瓣膜附近處NO 的濃度將沿著整個淋巴管,自上游向下游逐漸升高.反之,若淋巴管舒張時,瓣膜附近處NO 的濃度將自上游向下游逐漸降低.圖4(b)、圖4(d)和圖4(f)分別表示淋巴管的段數(shù)為3,4,5 段時壓力的分布情況,淋巴管外紅色點表示組織.如圖4(d)所示,四段管中的前兩段淋巴管管內(nèi)顏色為紅色,表示這部分壓力非常大,應(yīng)處于收縮的初始階段;后兩段淋巴管中壓力呈現(xiàn)出逐漸遞減的趨勢,證明了淋巴管的自發(fā)收縮從上游向下游逐漸過渡的過程.四段管整體處在收縮周期的初始階段,此時淋巴液流動速度較小,作用在瓣膜和淋巴管壁上的剪切力也較小,因而管內(nèi)部NO 的濃度整體偏低,故圖4(c)相較于圖4(a)和圖4(e)沒有出現(xiàn)NO濃度較高的紅色區(qū)域.圖4(b)和圖4(f)中,各管內(nèi)壓力強度偏低,且管間壓力分布也較為均勻.

圖4 淋巴管的段數(shù)不同時,(a),(c),(e) NO 濃度和(b),(d),(f)壓力分布(顏色標(biāo)尺頂部紅色部分最大壓強為3.045 × 103 Pa,底部藍色部分最小壓強為3.010 × 103 Pa) (a),(b) 3 段;(c),(d) 4 段;(e),(f) 5 段Fig.4.(a),(c),(e) NO concentration and (b),(d),(f) pressure distribution for different segment lymphatic vessel (The maximum pressure of the red part at the top of the color scale is 3.045 × 103 Pa,and the minimum pressure of the blue part at the bottom is 3.010 × 103 Pa):(a),(b) Three segments;(c),(d) four segments;(e),(f) five segments.

圖5 表示3—5 段淋巴管中的NO 平均濃度隨時間的變化情況.當(dāng)淋巴管開始收縮時,最左側(cè)瓣膜關(guān)閉,淋巴液的流動使作用在瓣膜處的剪切力變大,NO 的平均濃度相應(yīng)升高.當(dāng)淋巴管達到舒張期峰值時,淋巴液的流速變小,產(chǎn)生很少的NO,NO的降解或擴散使其平均濃度有所下降.繼而進入下一次收縮-舒張循環(huán),NO 的平均濃度呈現(xiàn)出周期性的增減變化,并且隨著循環(huán)周期的往復(fù),NO 的最低基準(zhǔn)濃度持續(xù)升高.特別的是,3 段管中平均濃度的峰值經(jīng)過多次收縮循環(huán)后由最低變成了最高,5 段管中平均濃度的峰值變成了最低.

圖5 NO 平均濃度隨時間變化Fig.5.Average NO concentration changes with time.

圖6 為3 段管中3 種淋巴管段的流量隨時間變化圖.如圖6 所示,在初始淋巴管和中間淋巴管中都出現(xiàn)了明顯的負流量峰.但在初始淋巴管中,首先出現(xiàn)負流量峰,之后出現(xiàn)了正流量峰.這是因為從組織處流入的組織液首先經(jīng)過初始淋巴管,使得初始淋巴管最先膨脹.淋巴管壁的拉伸進一步增加了淋巴肌細胞細胞質(zhì)中鈣離子的濃度.因此,多段淋巴管的收縮是從上游的淋巴管開始,并逐漸過渡到下游的淋巴管進行收縮.初始淋巴管最先開始收縮,因為在多孔入口處允許15%的淋巴液滲透回組織,這使管道中的向右的流體阻力大于向左的阻力,導(dǎo)致總體流量向左大于向右,因而出現(xiàn)負流量峰.當(dāng)初始淋巴管左側(cè)瓣膜完全關(guān)閉時,淋巴液只能向右流動,正流量峰繼而出現(xiàn).

圖6 三段管中流量隨時間變化Fig.6.Flux in three-segment lymphatic vessel changes with time.

對于中間的集合淋巴管,因為鈣離子的擴散方向為由左向右,故中間淋巴管的收縮方向與其保持一致,正流量峰出現(xiàn).鈣離子繼續(xù)擴散到下游,使得右側(cè)下游淋巴管開始收縮,負流量峰出現(xiàn).最右側(cè)的出口淋巴管,因為淋巴液可以輕松流出,所以出口淋巴管表現(xiàn)出峰值較高的正流量峰.

5 結(jié)論

通過LBM,建立了1 個多段淋巴管模型.基于NO 對淋巴管收縮的抑制作用,以及其產(chǎn)生依賴于剪切力的事實,NO 在淋巴管內(nèi)的分布對淋巴管的整體收縮起到重要作用.在活體的淋巴管中,瓣膜處的直徑最小,因而在該區(qū)域的剪切力相應(yīng)的更大.此外,瓣膜的瓣葉與淋巴管壁間淋巴液的流動意味著NO 可以從瓣葉的兩側(cè)釋放并集中在瓣膜附近.隨著淋巴管內(nèi)瓣膜數(shù)量的增多,特別是當(dāng)淋巴液回流時,管內(nèi)整體的NO 濃度分布變得十分復(fù)雜.但NO 的平均濃度始終隨著淋巴管的收縮-舒張循環(huán)呈現(xiàn)出周期性的增減變化,并且平均濃度隨著循環(huán)周期的往復(fù)而升高.從流量隨時間變化圖中可以得出,多段淋巴管的收縮是從上游的初始淋巴管開始,并逐漸過渡到下游的集合、出口淋巴管進行收縮.當(dāng)淋巴管收縮時,瓣膜附近處NO 的濃度將沿著整個淋巴管,自上游向下游逐漸升高;當(dāng)組織液流入,淋巴管舒張時,瓣膜附近處NO 的濃度將沿著整個淋巴管,自上游向下游逐漸降低.下一步,將在淋巴管中引入可形變細胞并施加多倍重力,通過細胞可以更直接地觀察管內(nèi)流體流動情況,進而更好地研究剪切力與NO 濃度的變化情況.