基于晶格玻爾茲曼方法研究不同出口壓力條件下淋巴管內氮含量的變化及影響*

赫軼男 張乾毅 韋華建 施娟

1) (桂林電子科技大學材料科學與工程學院, 桂林 541004)2) (桂林電子科技大學信息與通訊學院, 桂林 541004)(2019年12月22日收到; 2020年3月9日收到修改稿)

淋巴系統是人體內重要的防御功能系統, 具有三大免疫功能, 首先是能夠抵御細菌病毒, 使人體免于疾病的攻擊; 其次是由淋巴細胞加以輔助, 清除由新陳代謝而出的產物; 最后是由淋巴細胞來修補受損的器官與組織, 使其恢復正常的生理功能。淋巴系統沒有像血液循環系統中心臟一樣的動力泵, 淋巴液的驅動主要靠淋巴管的自主收縮來完成(肺淋巴系統是靠肺泡的運動)。淋巴管的自主收縮循環是由淋巴肌細胞內鈣離子增加產生收縮, 收縮驅動流體產生剪切力, 剪切力使淋巴內皮細胞產生一氧化氮合酶(eNOS), 一氧化氮合酶使一氧化氮增加, 一氧化氮的增加降低鈣離子使淋巴管松弛, 淋巴管松弛后流體剪切率下降, eNOS 下降,一氧化氮下降, 鈣離子增加, 淋巴肌細胞收縮, 開始新的周期。可見一氧化氮的濃度及其分布對淋巴管的收縮起關鍵作用。顯然出口壓力會影響淋巴管內流體的剪切率, 進而影響一氧化氮的濃度和淋巴管的收縮。為了研究淋巴管出口壓力對淋巴管收縮的影響, 建立了一個晶格玻爾茲曼模型, 模擬嵌入多孔組織的初始淋巴管和有兩對瓣膜的集合淋巴管, 該模型可以重現一氧化氮、鈣的相互影響以及淋巴管的自主收縮, 并研究不同出口壓力下一氧化氮的分布及其平均值.

1 引 言

淋巴系統由淋巴器官、淋巴液以及淋巴管道組成. 而淋巴管又分為毛細淋巴管、淋巴干及淋巴導管, 淋巴管作為淋巴液回歸血液循環的閉鎖管道,與人體眾多組織中的血液毛細血管平行. 淋巴液、組織液體與血漿三者之間存在著物質交換關系, 當血漿流經毛細血管時, 水和其他一些可以透過毛細血管壁的物質在毛細血管的動脈端滲出, 進入到組織細胞間隙中成為組織液, 而絕大部分的組織液在毛細血管靜脈端又重新滲入到血漿中. 少量的組織液還可以滲入到毛細淋巴管, 形成淋巴, 之后淋巴經淋巴循環由左右鎖骨下靜脈匯入血漿中. 淋巴管允許組織液從間質間隙擴散到淋巴毛細血管, 且淋巴毛細血管可以吸收小腸中的短鏈脂肪酸. 弗吉尼亞大學醫學院的Louveau 等[1]發現大腦會通過淋巴管道與免疫系統直接聯結, 由于這些淋巴管獨特的位置, 致使解剖學家這么多年未發現其位置, 這項新發現直接顛覆了幾十年來寫在教科書上的結論. 這對阿爾茨海默氏癥、自閉癥、帕金森癥、多發性硬化癥以及其他大腦疾病有著重大影響, 由此可知, 數值模擬對于研究淋巴管內協調作用有著極其重要的意義. 早期的淋巴系統網絡模擬模型, 是一種簡單明了的一維模型, 該種模擬基于流體力學的Navier-Stock(NS)方程, 由Margaris 和Black建立[2]. 而后Macdonald 等[3]對該模型進行了改進, 只模擬了小段淋巴管. 淋巴管泵送的復雜模型可以解釋振蕩行為, 并且這些行為對局部的結構以及壓力非常敏感.

研究表明, 一氧化氮 (NO) 作為生物活性分子通過cGMP/cAMP 及其依賴的蛋白激酶(PKA/PKG)作于KATP, 以周期性的變化參與生理狀態下淋巴管的收縮、舒張以及張力調節[4,5]. 哈佛醫學院Kunert 等[6]和Baish 等[7]闡明了淋巴管泵送行為是如何通過動力學反饋在整個淋巴系統中進行協調的. NO 與Ca2+濃度建立了自我調節的震蕩反饋環路, 維持正常的淋巴功能. 然而, 這些研究在進行模擬時, 腔內瓣膜是在逆向流動期間以數學的方式插入不可滲透的壁而施加的. 雖然這種方法能夠正確重現在體內觀察到的行為, 但它簡化了瓣膜性能, 從而忽略了瓣膜材料結構與力學性能以及產生的NO. 實際研究發現瓣膜是NO 的主要來源.

2 晶格玻爾茲曼方法及其基本理論

2.1 基本理論

在計算流體動力學方面, 晶格玻爾茲曼方法(LBM)是一個非常好的工具. 迄今為止, LBM 已經成功地應用到了多種復雜流體系統, 如多孔介質流[8]、多相流[9]、淋巴流[10]等.

LBM 最顯著的特征是它將波爾茲曼方程離散化, 與傳統的計算流體動力學方程不同的是,LBM 形式簡單, 不需要像宏觀方法那樣, 為了滿足連續性方程而在每一個時間步上都求解NS 方程, LBM 還可以將復雜的邊界條件引入. LBM 的基本思想是構建簡化的動力學模型, 遵從微觀物理過程的本質, 使宏觀特性服從宏觀方程.

使用單粒子分布函數fi代替布爾變量,LBM 演化方程為

其中fi(x,t) 是在x位置,t時刻, 具有ei速度的粒子分布函數.?(fi) 是碰撞因子, 表示碰撞對于fi的影響.

引入單弛豫時間近似[11], 碰撞項被簡化為

求出粒子分布函數后, 遵循質量與動量守恒,格點上的質量和流體的流速可以用下式計算:

其中ρ是格點上流體的質量,u是格點上流體的流速.

二維九速(D2Q9)模型如圖1 所示. 利用Chapman-Enskog 多尺度展開[12], 得到局域平衡分布函數:

其中ωi有下列形式:

圖1 D2 Q9 晶格玻爾茲曼模型的微觀速度Fig. 1. Microscopic velocity of D2 Q9 lattice Boltzmann model.

2.2 反向彈回邊界條件

在做流體力學方面的研究時, 流體相關的規律與性質需要由邊界條件來加以控制. 而LBM 的基本變量是分布函數, 故LBM 在構造邊界條件時,需要結合流體的宏觀條件. LBM 的所有邊界條件都是基于反彈邊界條件修改而來.

反向彈回邊界條件是一種無滑移邊界條件, 如圖2 所示, 該模型被簡單地分為兩層, 白點層代表固體, 在上面的黑色點是距離邊界層最近的流體點. 我們假設流體粒子從A,B,C三個位置沿著速度方向1, 3, 5 流動到邊界上位置D, 當分布函數到達后, 在下一時間步流體又從位置D沿著相反的速度方向2, 4, 6 流回到流體層, 此過程滿足質量守恒定律, 且邊界處的平均速度是0.

圖2 反向彈回示意圖Fig. 2. Bounceback.

部分反彈邊界條件[13]可以用來改變流體流動的阻力. 調整2, 4, 6 的分布函數中從1, 3, 5 反向彈回和前方流體點流過來的2, 4, 6 方向分布函數的比例, 即可以使邊界在沒有流動和自由流動之間改變.

3 數值模擬結果與討論

3.1 淋巴管模型

我們建立的淋巴管模型是由兩個瓣膜限定的具有嵌入組織的單個淋巴管段, 如圖3 所示. 在完全松弛狀態下, 淋巴管的直徑為100 μm, 即2R0=100 μm. 組織液可以從左側具有固定幾何形狀的多孔管道進入, 而這段管道代表了初始淋巴毛細血管[14]. 為了模擬淋巴毛細血管中的主要瓣膜, 該模型的實線區域是不可滲透的, 而虛線區域是半滲透的, 使用部分反彈邊界條件實現[13]. 我們使用恒定的反彈比ξ=0.85 , 這表示當壓力條件有利于在此方向上的流動時, 有15%的流體可以滲透回組織.滲透部分和不可滲透部分的面積比例為1∶1, 而淋巴毛細血管瓣膜可以使液體流入, 但禁止回流. 淋巴管的右端表示出口, 將與下游淋巴管或淋巴結連接. 在這部分, 我們引入長為234 μm 的固定端, 用來提供結構支撐, 以減少右側邊界的影響.

圖3 淋巴管段示意圖Fig. 3. Lymphatic section.

我們使用LBM 來計算流體的流動、剪切力、組織和淋巴管中的壓力[11]. 上下淋巴管壁被分解為200 段, 瓣膜分解為32 段, 每段只沿y方向運動淋巴管內壁上的內皮細胞可以產生NO[15], 且隨著淋巴管內壁上剪切應力的增加而增加. 產生的NO 會在流體和組織中對流和擴散:

其中CNO是 NO 的濃度,DNO是NO 的擴散系數,h是化學反應速率常數.h?t決定了化學反應時間,?t是格子時間, 隨著h的增加, 化學反應速度增加.和是NO 的衰減率和生成系數. 我們認為, NO 的產量與切線方向的應力成正比[16].υl表示淋巴管內表面的切線方向的流體速度,xn表示沿著法線方向的梯度. 在淋巴管內, 瓣膜也可以產生NO.

淋巴管的收縮與松弛由淋巴肌細胞的Ca2+濃度決定, Ca2+可以進入、離開淋巴肌細胞的細胞質,也可以通過連接點到達相鄰細胞, Ca2+的反應擴散方程如下:

其中CCa是Ca2+的濃度,DCa是Ca2+從一個細胞擴散到鄰近細胞的擴散系數,是Ca2+的衰減率,是Ca2+的產生速率,是淋巴管擴張產生Ca2+的非線性項[17],R和分別是淋巴管的局域半徑和非線性項的參考半徑,Rl是管子的最小半徑,δ↑是一個非對稱Kronecker函數, 如果CCa從低于閾值Cth增加到Cth, 則該函數被設定為1, 其他情況被設置為0.

3.2 淋巴管壁所受的力以及瓣膜

有五個力施加在淋巴管壁上: 流體力F、淋巴肌力FM、彈性力FE、彎曲力FB、黏性阻力Fr.

1)流體力F可以通過壓力張量積分求得:

其中δij是Kronecker 函數且i=x,j=y,ρ是格點上流體的質量,u是格點上流體的流速, 這里我們使用到面積 ds最近的流體點的分布函數來代替fi.n是面積 ds上的單位矢量,us是面積 ds的速度.

2)淋巴肌力FM與上文提及的Ca2+濃度和NO 的濃度有關:

其中kM是決定力的常量系數.

3)彈性力FE來源于淋巴管組織:

為了限制收縮幅度和模擬組織機械阻力, 如果R

4)彎曲力FB我們假設淋巴管的左右兩側是固定, 淋巴管的各段可以沿y方向移動:

其中i表示某一離散片段.

5)黏性阻力Fr在進行計算時需要考慮到淋巴管的黏彈性, 且黏性阻力作用于壁面速度ν相反的方 向上:

計算中用兩塊可變形的拋物面來模擬淋巴瓣膜, 如圖4 所示.

圖4 靜止狀態下淋巴管瓣膜Fig. 4. The lymphatic valves at rest.

淋巴瓣膜在靜止狀態下偏向出口位置. 因此,我們將瓣膜的靜態形狀設置為拋物線, 方程為:

其中y0是瓣膜靜止或極限位置,H是靜止狀態下淋巴管的位置,x0是瓣膜的錨點, 負正號分別代表上下瓣膜,A起到控制瓣膜張角的作用.

其中i表示瓣膜的某一離散片段.

計算彈性力時, 還需考慮當瓣膜的兩個薄膜極為接近時, 需要乘以一個11 冪項來增加

其中yc代表淋巴管的中心線位置.

當瓣膜達到極限位置ylp時, 也同樣需要增加

當兩個瓣膜相互接近從而導致缺失流體格點,無法使用LBM 計算流體力時, 可使用潤滑力來計算流體力[18].

3.3 計算參數

我們采用型號為Nvidia Quadro GP100 的專業級顯卡, 該顯卡具有較高的性價比, 擁有著16 G 的顯卡內存, 以及3584個CUDA 核心, 一個專業級顯卡所具備的CUDA 核心直接反映著計算速度的快慢.

在晶格玻爾茲曼方法的計算中, 需使用到量綱轉換, 即實際物理量和相應的離散格子量之間的轉換. 這里選擇τ=0.75 ,u= 0.01 cm2/s,D′=25和D= 0.01 cm, 其中D′和D分別代表著格子上淋巴管直徑與實際上的淋巴管直徑. 因此, 實際時間和空間尺度分別為T= 1.33 × 10–6s,L=0.0004 cm. 在本次計算中, 使用到的計算參數如表1 和表2 所列.

表1 Ca2+與NO 的化學參數Table 1. Chemical parameters of Ca2+ and NO.

表2 淋巴管與瓣膜參數Table 2. Parameters of Lymphatic and valve.

進行模擬時, 因為淋巴管壁和瓣膜的質量不同, 所以我們需要在格子上使用不同的密度. 瓣膜的密度取1 格子單位, 淋巴管壁密度取80 格子單位. 所有的流體節點初始密度設定為1 格子單位,初使速度設定為0, NO 的濃度設定為0. 通過施加速度為0 的平衡分布, 邊界處的密度保持為常數ρin, 出口處的密度通過施加恒定壓力邊界條件也保持恒定值ρout[18]. 邊界上的NO 濃度保持恒定值0, 且NO 可以通過流體、瓣膜結構和淋巴管壁擴散. 在我們的模擬中, 保持ρin= 1. 由于淋巴液被視為水, 因此壓力單位為P= (L/T)2g·cm–3=9.045 × 104g·cm?1·s–2, 而 入 口 壓 力 保 持 在3.015×104g·cm?1·s?2.

3.4 淋巴管出入口壓力差對NO 濃度的影響

我們的模擬計算表明淋巴管收縮取決于NO 濃度的變化, 而NO 濃度變化情況取決于流體剪切率的變化. 因此從組織到淋巴管的壓差可改變流體速度和泵送狀態. 為了了解淋巴管對可能遇到的組織壓力變化的反應, 我們設計了壓力階躍變化的方案. 即保持入口密度不變ρin=1 g/cm3, 通過改變出口處密度, 進而改變壓力差.

計算參數如表3 和表4 所列.

圖5 是時間t= 2.296 s時管內的NO 分布.注意壓差不同淋巴管的收縮周期不同[10], 所以同時刻淋巴管所處的狀態不同. 如圖5 所示, DP=30.15 g·cm–1·s–2, 此時淋巴管開始收縮, 左瓣膜關閉, 右瓣膜打開, 由此產生壁面剪切率變化影響到NO 的產量. 瓣膜的運動進一步增加了瓣葉表面上的剪切應力, 在瓣膜附近產生較多的NO. 由于回流的原因, 一些NO 從左瓣膜流出, 高濃度的NO 從打開的右側瓣膜排出, 從右側瓣膜周圍的NO 的凸型狀可以看出在瓣膜附近NO 濃度較高,Glenn 等[19]已在小白鼠的淋巴管中已經觀察到此現象.

如圖5 所示, DP= 18.09 g·cm–1·s–2, 此時, 淋巴管開始松弛, 左瓣膜打開, 右瓣膜關閉. 流體從左側流入淋巴管, 與右側相比, 該區域產生更多NO. 即使淋巴管松弛后, NO 的濃度仍然較高, 當淋巴管達到舒張期峰值時, 流體速度變小, 并且產生很少的NO. NO 迅速降解或擴散, 其濃度下降.同時, 淋巴管準備進行下一次收縮. 收縮產生較高的剪切應力, 從而產生NO, 使其充滿了淋巴管.

由于瓣膜小葉之間的有效淋巴管直徑較小, 瓣膜處的剪切應力較大, 且瓣膜的兩個表面又都能產生NO, 這便是NO 的濃度較高的原因. 在左側的初始淋巴管段中, 帶有單向閥的多孔壁處會發生一些回流, 即在連通區域有近乎85% 的流體被反射回來.

表3 出口壓強高于入口壓強時正壓力差Table 3. Positive pressure when outlet pressure is higher than inlet pressure.

表4 出口壓強低于入口壓強時負壓力差Table 4. Negative pressure when outlet pressure is lower than inlet pressure.

圖5 t = 2.296 s時, NO 濃度分布圖Fig. 5. t = 3.003 s, NO concentration distribution map.

與預期結果相符, 淋巴管在壓力差為負值時,即出口壓力低于入口. 在這種情況下, 壓力也可以驅動淋巴液流動. 壓差高到一定程度時, 內皮表面上的剪切應力產生了高濃度的NO, 如圖6 所示.且高濃度的NO 會降低淋巴管中Ca2+的濃度, 使淋巴管停止收縮. 此時淋巴管的半徑小于靜止時的半徑, 這意味著淋巴管內的壓力低于外部. DP<–10 g·cm–1·s–2時, NO 達到飽和, 不再增加.

當淋巴管壓力差為正值時, 即出口壓力高于入口. 在這種情況下, 由于回流, 高壓迫使一些流體返回淋巴管, 使整段淋巴管膨脹. 同時Ca2+的濃度增加, 激發下一次收縮, 迫使淋巴液流出淋巴管.由于出口壓力增加, 淋巴管的收縮也無法提升平均流量, 因此剪切率下降造成NO 的平均值下降. 且當 ?P過高后會造成出口瓣膜無法打開, 收縮沒有周 期性.

圖6 NO 平均濃度與壓強差關系圖Fig. 6. Relationship between NO average concentration and pressure difference.

4 結 論

我們使用晶格玻爾茲曼方法模擬了淋巴管段的收縮, 同時采用了基于淋巴管結構力學和NO 的動力學模型, 最終再現自我持續循環的淋巴收縮并輸運組織液. 計算結果表明NO 的平均含量隨著淋巴管兩端壓力差的變化而變化, 在負壓強差時, 隨著壓強差絕對值的增加趨于平穩. 而在正壓強差時, 隨著壓強差絕對值的增加先降低而后趨于穩定. 在負壓力差的情況下, 由于入口壓力高于出口壓力, 只要流體達到一定的流速, 這時淋巴管內可以產生足夠的NO, 抑制了鈣離子, 淋巴管肌細胞停止收縮, 由壓差驅動流體實現收縮的最小化. 反之如果出口壓力過高, 收縮也無法令瓣膜打開, 流體無法正常輸運, 淋巴管內NO 和鈣離子就無法產生有規律的振蕩, 淋巴管收縮失去規律性. 由于在淋巴管內細胞也可以產生NO, 另外瓣膜還會存在一定的滲漏, 也會產生NO, 而且癌變細胞可以產生更高的NO, 所有這些影響都是我們后續將會進一步進行研究.