基于晶格玻尔兹曼方法研究不同出口压力条件下淋巴管内氮含量的变化及影响*

赫轶男 张乾毅 韦华建 施娟

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, 所有这些影响都是我们后续将会进一步进行研究.