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

张乾毅 韦华健 李华兵

(桂林电子科技大学材料科学与工程学院,桂林 541004)

淋巴系统对人体的免疫及细胞的内环境稳态都有着重要的作用.与血液循环系统相似,淋巴系统也是遍布全身的管道系统,主要由淋巴液、淋巴管和淋巴器官构成.淋巴管的自发收缩驱动管内淋巴液的流动.淋巴管的自发收缩-舒张机制由Ca2+和NO 浓度的振荡反馈决定,NO 在管内的分布对淋巴管的收缩循环起到重要作用.因淋巴液流动而作用在淋巴瓣膜上的剪切力是瓣膜产生NO 的主要原因.在真实系统中,某段淋巴管中的NO 分布会受到与其连接的淋巴管的影响,特别是上游的连接片段.通过晶格玻尔兹曼方法,建立了1 个具有瓣膜结构的多段淋巴管模型,再现了淋巴管内Ca2+和NO 的反馈机制,瓣膜变化和淋巴液流动情况.该模型中存在3 种淋巴管,分别是初始淋巴管、中间淋巴管和出口淋巴管.淋巴管的段数可以通过修改计算参数无限扩充.本文计算的段数为3—5 段,每段淋巴管中有两对瓣膜.通过模型研究了多段淋巴管中NO 浓度分布、压力分布、NO 平均浓度变化,以及3 段管模型中各管的流量随时间变化情况.

1 引言

淋巴瓣膜出现在集合淋巴管中,是毛细淋巴管向集合淋巴管过渡的标志,淋巴瓣膜将集合淋巴管分成了若干小段,起到防止淋巴液回流的作用[1].与血液循环系统中,将心脏作为单一泵源并通过没有激励源的血管进行循环的方式不同,淋巴管表现为一组串联排列的分布式泵[2],每段被瓣膜分隔开的独立淋巴管段都具有自发收缩驱动淋巴液流动的能力[3-6].集合淋巴管可以根据外部压力条件,自行调整其收缩行为以维持体液稳态,若存在外部的流体压力梯度可以驱动管内淋巴液流动,集合淋巴管可最小化收缩,否则便主动收缩驱动淋巴液流动[7].

淋巴管的自主收缩机制为,当淋巴管中Ca2+浓度超过淋巴肌细胞的阈值浓度时,淋巴肌细胞中会形成蛋白质间的交叉连接,从而引起自主收缩.淋巴管的收缩导致淋巴液流动并产生剪切力,淋巴管内皮细胞响应动态流体剪切力并激活其中的一氧化氮合成酶(eNOS)产生并释放一氧化氮(NO)[8],引起淋巴管舒张.由此可见,是NO 和Ca2+浓度的周期性变化引起了淋巴管的收缩-舒张过程.

相较于动脉有二十多个模型可供参考,淋巴系统的数值模拟还处在初级阶段.最早的淋巴系统网络模拟是由Reddy 和Black[9]根据Navier-Stokes方程设计的简化一维模型.后来,Bertram 等[10,11]用集总方法模拟了小段淋巴管,发现瓣膜阻力随着施加在瓣膜上的压力梯度的变化而变化,但其瓣膜是虚拟的,只起到单向作用.研究表明,瓣膜是产生一氧化氮的重要来源.在前期研究中,建立了存在真实瓣膜的集合淋巴管模型,但其存在着明显局限性,只有1 个单独的淋巴管段.而在真实的淋巴系统中,淋巴管段之间首尾连接,某段淋巴管中的NO 的分布会受到其上下游相邻管段的影响,本文在之前模型的基础上,建立了1 个可以调整淋巴管段数目的多段淋巴管模型.

2 晶格玻尔兹曼方法

晶格玻尔兹曼方法(lattice Boltzmann method,LBM)是1 种介观尺度的模拟方法,在流体力学方面的模拟中被广泛应用,如多孔介质流[12]、血液流[13]、多相流[14]等.LBM 以微观的流体粒子分布函数为基础,用分布函数对速度集求各阶矩得到相关的宏观量,避开了传统数值方法直接求解宏观变量所带来的复杂非线性项,大大简化了计算的难度.在流体碰撞和流动的过程中,粒子总体满足质量和动量守恒,粒子之间的相互作用由局部分布函数决定.

Mcnamara 和Zanetti[15]用单粒子分布函数fi取代了格子气细胞自动机中的布尔变量,其LBM演化方程为

其中i表示微观速度的索引,对于D2Q9 模型(见图1),i=0,1,2,···,8;fi(x,t) 是在x位置、t时刻,具有ei速度的粒子的分布函数,δt为时间步长,δx为格子长度;Ω(fi)是碰撞因子,表示碰撞对fi的影响.Chen 等[16]和Qian 等[17]提出用单弛豫时间来替代碰撞因子项.继而,LBM 方程可以写为

图1 D2Q9 模型Fig.1.D2Q9 model.

其中,ρ是宏观密度,u是格点上流体的流速.本文采用LBM 中的D2Q9 模型[18]进行计算,其平衡分布函数具体定义如下式:

式中,c为基准速度,通常情况下取1.ωi为各方向的权重系数,当i=0 时,ωi=4/9 ;当i=1,2,3,4 时,ωi=1/9;当i=5,6,7,8 时,ωi=1/36.ei的形式为

边界条件的选择对于LBM 的计算效率、计算精度和稳定性都起到非常重要的作用.因为在迭代过程中,经过格子的碰撞和流动迁徙过程之后,流场内格点的分布函数已经更新,而边界上仍存在未更新的格点分布函数.因而需要用确定的边界条件,并采用对应的边界处理格式对未更新的格点分布函数进行确定,才能继续进行下一步的迭代.

压力边界条件用来处理存在压差条件下,流体的流动以及边界滑移时流体的流动.如图2 所示,在左入口处,经过流动迁徙后指向流场内部的3 个分布函数f1,f5和f8为未知量.相同的,在右出口处f3,f6和f7也为未知量.需要拟定好出入口密度值,并假设出入口处速度的垂直分量uy=0,并进一步求出出入口速度水平分量ux和未知的分布函数.根据(3)式、(4)式可以得出以下关系[13]:

图2 压力边界条件示意图Fig.2.Schematic diagram of pressure boundary condition.

用LBM 对实际物理体系进行模拟时,必须对模型中的物理量进行量纲转换,即实际物理量与离散格子量之间的转换.在本文的计算中,用L表示长度量纲,T表示演化的时间量纲,G表示质量量纲.在各量纲符号后加上“′”代表格子上的量纲.将颗粒的直径D,流体密度ρ和黏滞系数v作为基准物理量,则:.根据这些基准物理量,可继续求得其他物理量.

3 淋巴管模型

在真实的淋巴系统中,淋巴管段间彼此连接.某段淋巴管段的收缩情况及其内部NO 和Ca2+的分布情况都受到其他连接的淋巴管段的影响,特别是其上游连接管段.建立的多段淋巴管模型包含有3 种淋巴管段,总段数可以通过修改MPI 计算参数进行调整.每种淋巴管段中都含有两个腔内瓣膜和淋巴管壁.如图3(a)所示为初始淋巴管,其直径为0.01 cm,虚线部分表示具有可供组织液自由流入的多孔入口.在多孔入口处设定恒定的反弹比,只允许初始淋巴管中15%的淋巴液渗透回组织.图3(b)为中间淋巴管.图3(c)为出口淋巴管,其尾端表示0.0296 cm 长的静态淋巴管壁,用来降低出口处淋巴管壁的影响.淋巴管段的外部为组织,假设组织处不会流动,只有当淋巴管膨胀,淋巴管壁的位置越过组织时,组织才会转化成流体.当组织转化为流体时,新产生的流体将由其相邻的流体计算得出.通过LBM 中的D2Q9 模型来模拟淋巴管内外部的流体.淋巴管壁、淋巴瓣膜通过曲线边界条件实现,流体作用力通过压力张量积分法计算.组织处、初始淋巴管的入口和出口淋巴管的静态部分通过反向弹回边界条件实现.最外层的黑色加粗壁由平衡分布函数实现.其速度为0,密度为1.为了保护出口淋巴管中尾端的瓣膜,在距离出口0.0016 cm处引入了1 个部分反向弹回边界,流体可以自由流出,若外部流体有从出口处流入的趋势,该边界会阻止85%的流体流入.若尾端的瓣膜完全关闭,部分反弹变为反向弹回.

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

模型中淋巴管壁和瓣膜均被离散化,只能沿着竖直方向移动,二者均受弯曲力作用,可表示为

其中,ym和yn是格子相邻两点y轴上的位置.KB在淋巴管中是个常数,而在瓣膜中用带有上标v(表示瓣膜valve)的表示瓣膜的抗弯强度,此处抗弯强度具体表示为

其中,i表示当前段数,n表示瓣膜的总段数.是在固定点位置处的最大值,是最小值的近似值.系数A调节瓣膜的柔软度,瓣膜在淋巴液回流时要有足够的强度抵挡流体的冲击,因此瓣膜的柔性系数由淋巴管的中间线附近的自由点处向淋巴管壁固定点处逐渐增大.

弹性力FE来自淋巴管外组织处:

式中,R表示淋巴管和瓣膜的半径,R0在淋巴管中是个常数.但在瓣膜中,为了保持瓣膜的抛物线形状,假设瓣膜在初始位置的抛物线形状表示为

式中,yl0表示瓣膜的初始位置,x0表示瓣膜的固定点水平方向位置,正负号分别表示上下瓣膜,B用于调整初始位置是偏向于闭合或者张开.淋巴管的最小半径为Rl.为避免瓣膜张开程度过大,用一个最大B值确定的位置为其极限位置.

因淋巴管和瓣膜都具有黏弹性,黏滞阻力可表示为

负号表示Fr总是作用于壁面沿速度v的反方向.

当两个瓣膜位置接近时,因为缺少流体格点而不能使用LBM 来计算流体力.在这种情况下,使用润滑理论来计算流体力[19],又因为瓣膜只在y方向上运动,只考虑垂直方向上的润滑力,可表示为

式中,l和h分别表示瓣膜片段的长度和瓣膜片段间的距离,va是该段瓣膜的速度.

除了上述弯曲力、弹性力和黏滞阻力外,作用于淋巴管上的力还包括淋巴肌力和流体力.淋巴管由淋巴肌细胞组成,当大量钙离子扩散到淋巴肌细胞细胞质时,会触发淋巴肌细胞收缩.淋巴肌力与Ca2+和NO 的浓度相关,可表示为

式中,CCa和CNO分别表示钙离子和一氧化氮的浓度,KM是决定淋巴肌力强度的常量.Ca2+的反应扩散方程表示为

式中,DCa是Ca2+的扩散系数,是Ca2+的衰减速率,基于NO 通过肌球蛋白轻链磷酸酶作用于平滑肌细胞可以减少收缩力的产生,Ca2+的总衰减速率需乘上 (1+KCa,NOCNO)[5].是Ca2+的生成速率,11 次项是淋巴管拉伸产生Ca2+.δ ↑是非对称Kroneckerδ函数,当CCa从低于阈值Cth增加到Cth时,该函数设定为1,其他情况设定为0;当CCa达到阈值Cth时,它触发由电压门控和钙诱导的钙通道介导动作电位[20,21];表示当CCa＞Cth时CCa的急剧增大.R和RCa分别表示淋巴管的局部半径和非线性项参考半径,Rl是淋巴管的极限最小半径.一氧化氮的产生和扩散可以表示为

式中,CNO表示NO 的浓度,DNO表示NO 的扩散系数,分别表示NO 的衰减速率和生成速率.

4 计算结果与分析

本次计算代码使用C++编写,计算代码在搭载了8 个GPU 的服务器上执行.GPU 型号为NVID IA Quadro GP100,每个GPU 含有3584 个CUDA核心,处理双精度浮点数的能力为5.2 TFLOP/s,搭配16 GB HBM2 显存,理论带宽高达717 GB/s.编译环境为CUDA10.0.每个GPU 用来计算一段淋巴管,淋巴管数分别为3,4,5 段.MPI 用来交换淋巴管段边界间的数据.通过量纲转换将实际淋巴管物理参数与格子物理量进行转换,其详细参数参考文献[22].

图4 表示多段管中的NO 浓度分布和压力分布情况,右侧颜色标尺自下而上的颜色变化表示管内NO 浓度和压力强度的增大.图4(a)、图4(c)和图4(e)分别表示淋巴管的段数为3,4,5 段时NO 的浓度分布情况,瓣膜是产生NO 的主要结构,瓣膜附近的NO 浓度要明显高于其他部分.随着瓣膜数量的增加,NO 在整个管腔内的分布更为复杂,Ca2+的分布变得越来越随机,进而导致Ca2+的浓度达不到淋巴肌细胞收缩的阈值,淋巴管整体的收缩情况也变得混乱.此外,瓣膜附近NO 的浓度还受淋巴管当前收缩阶段的影响.当淋巴管收缩时,瓣膜附近处NO 的浓度将沿着整个淋巴管,自上游向下游逐渐升高.反之,若淋巴管舒张时,瓣膜附近处NO 的浓度将自上游向下游逐渐降低.图4(b)、图4(d)和图4(f)分别表示淋巴管的段数为3,4,5 段时压力的分布情况,淋巴管外红色点表示组织.如图4(d)所示,四段管中的前两段淋巴管管内颜色为红色,表示这部分压力非常大,应处于收缩的初始阶段;后两段淋巴管中压力呈现出逐渐递减的趋势,证明了淋巴管的自发收缩从上游向下游逐渐过渡的过程.四段管整体处在收缩周期的初始阶段,此时淋巴液流动速度较小,作用在瓣膜和淋巴管壁上的剪切力也较小,因而管内部NO 的浓度整体偏低,故图4(c)相较于图4(a)和图4(e)没有出现NO浓度较高的红色区域.图4(b)和图4(f)中,各管内压力强度偏低,且管间压力分布也较为均匀.

图4 淋巴管的段数不同时,(a),(c),(e) NO 浓度和(b),(d),(f)压力分布(颜色标尺顶部红色部分最大压强为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 平均浓度随时间的变化情况.当淋巴管开始收缩时,最左侧瓣膜关闭,淋巴液的流动使作用在瓣膜处的剪切力变大,NO 的平均浓度相应升高.当淋巴管达到舒张期峰值时,淋巴液的流速变小,产生很少的NO,NO的降解或扩散使其平均浓度有所下降.继而进入下一次收缩-舒张循环,NO 的平均浓度呈现出周期性的增减变化,并且随着循环周期的往复,NO 的最低基准浓度持续升高.特别的是,3 段管中平均浓度的峰值经过多次收缩循环后由最低变成了最高,5 段管中平均浓度的峰值变成了最低.

图5 NO 平均浓度随时间变化Fig.5.Average NO concentration changes with time.

图6 为3 段管中3 种淋巴管段的流量随时间变化图.如图6 所示,在初始淋巴管和中间淋巴管中都出现了明显的负流量峰.但在初始淋巴管中,首先出现负流量峰,之后出现了正流量峰.这是因为从组织处流入的组织液首先经过初始淋巴管,使得初始淋巴管最先膨胀.淋巴管壁的拉伸进一步增加了淋巴肌细胞细胞质中钙离子的浓度.因此,多段淋巴管的收缩是从上游的淋巴管开始,并逐渐过渡到下游的淋巴管进行收缩.初始淋巴管最先开始收缩,因为在多孔入口处允许15%的淋巴液渗透回组织,这使管道中的向右的流体阻力大于向左的阻力,导致总体流量向左大于向右,因而出现负流量峰.当初始淋巴管左侧瓣膜完全关闭时,淋巴液只能向右流动,正流量峰继而出现.

图6 三段管中流量随时间变化Fig.6.Flux in three-segment lymphatic vessel changes with time.

对于中间的集合淋巴管,因为钙离子的扩散方向为由左向右,故中间淋巴管的收缩方向与其保持一致,正流量峰出现.钙离子继续扩散到下游,使得右侧下游淋巴管开始收缩,负流量峰出现.最右侧的出口淋巴管,因为淋巴液可以轻松流出,所以出口淋巴管表现出峰值较高的正流量峰.

5 结论

通过LBM,建立了1 个多段淋巴管模型.基于NO 对淋巴管收缩的抑制作用,以及其产生依赖于剪切力的事实,NO 在淋巴管内的分布对淋巴管的整体收缩起到重要作用.在活体的淋巴管中,瓣膜处的直径最小,因而在该区域的剪切力相应的更大.此外,瓣膜的瓣叶与淋巴管壁间淋巴液的流动意味着NO 可以从瓣叶的两侧释放并集中在瓣膜附近.随着淋巴管内瓣膜数量的增多,特别是当淋巴液回流时,管内整体的NO 浓度分布变得十分复杂.但NO 的平均浓度始终随着淋巴管的收缩-舒张循环呈现出周期性的增减变化,并且平均浓度随着循环周期的往复而升高.从流量随时间变化图中可以得出,多段淋巴管的收缩是从上游的初始淋巴管开始,并逐渐过渡到下游的集合、出口淋巴管进行收缩.当淋巴管收缩时,瓣膜附近处NO 的浓度将沿着整个淋巴管,自上游向下游逐渐升高;当组织液流入,淋巴管舒张时,瓣膜附近处NO 的浓度将沿着整个淋巴管,自上游向下游逐渐降低.下一步,将在淋巴管中引入可形变细胞并施加多倍重力,通过细胞可以更直接地观察管内流体流动情况,进而更好地研究剪切力与NO 浓度的变化情况.