胡晟^1,2, 叶俊彦¹, 赵勇^1,2
1. 东北大学秦皇岛分校 控制工程学院, 河北 秦皇岛 066004;
2. 东北大学秦皇岛分校 河北省微纳精密光学传感与检测技术重点实验室, 河北 秦皇岛 066004
收稿日期：2022-03-18
基金项目：国家自然科学基金资助项目(61903069)；河北省自然科学基金资助项目(F2020501040)；中央高校基本科研业务费专项资金资助项目(N2223034)。
作者简介：胡晟(1984-)，男，云南景洪人，东北大学秦皇岛分校副教授。
通信作者：HU Sheng, E-mail: husheng@neuq.edu.cn。

摘要：由于红细胞具有双凹形结构和较好的超弹性响应，其能通过大尺度的拉伸和收缩变形穿梭于细长的毛细血管，为维持人体生命活动所需的氧气承担了重要输运工作.本文采用COMSOL有限元软件，依托流-固耦合模块进行了红细胞在不同粗细血管、血浆黏度，以及血液流速等关键影响因素的动态模拟.通过仿真计算可知红细胞能够轻松通过3 μm孔径的毛细血管，相比6 μm孔径的毛细血管承受了4.5倍的流体剪切力.同时5.5 mPa·s黏度的血浆诱使红细胞形成降落伞形状，并且中间凹陷结构消失逐渐沿着流速方向水平凸出.本文对T型和Y型微通道的流速进行了研究，两者的仿真结果都指出非对称流体环境使红细胞形成镰刀弯钩形状.另外，较高流速的支路通道在汇聚口较快诱发红细胞弯钩形状消失，逐渐演变为呈直杆形平躺流动.
关键词：血管芯片红细胞流体动力学超弹性材料COMSOL
Study on Morphology of Red Blood Cell in Micro-channel Towards Vessel-on-Chip
HU Sheng^1,2, YE Jun-yan¹, ZHAO Yong^1,2
1. School of Control Engineering, Northeastern University at Qinhuangdao, Qinhuangdao 066004, China;
2. Hebei Key Laboratory of Micro-nano Precision Optical Sensing and Measurement Technology, Northeastern University at Qinhuangdao, Qinhuangdao 066004, China

Abstract: Red blood cell (RBC), undertaking the significant work of oxygen transport to maintain human life, could shuttle spindly capillaries by their own stretching and contractile properties due to their biconcave discoid shape and hyperelastic responses. In this paper, the coupled fluid-solid module in finite element software COMSOL was used to study the dynamic simulation of RBC when the three factors, including capillary width, plasma viscosity, and blood flow rate were considered. The results implied that the RBC could easily pass through capillary with hole size of 3 μm, and the aqueous shear stress was 4.5 times greater than that in the capillary with 6 μm hole size. Meanwhile, the plasma viscosity of 5.5 mPa·s induced parachute-like formation of RBC. Furthermore, biconcave structure in the center gradually disappeared and turned into bulge along the flow direction. With regard to the study on fluidic velocities in T-shaped and Y-shaped micro-channel, both results indicated that the asymmetric structure caused the falcate shape of RBC. In addition, the higher fluidic velocity in the other branch channel, the easier the convergence leaded to the disappearance of falcate shape and the rod posture lying flat flowing in the solution.
Key words: vessel-on-chipred blood cell (RBC)hydrodynamicshyperelastic materialCOMSOL
红细胞(red blood cell, RBC)是哺乳动物胚胎发育过程中最早产生的造血细胞.胚胎期的红系造血对于哺乳动物的妊娠过程及个体出生后的生命活动具有重要意义.同时，红细胞是血液中数量最多的细胞，约占人体细胞总数的1/4，其贯穿于人体全身的组织器官来行使重要生理功能——输送氧气.所以，红细胞独特的双凹形结构，使其能够最大限度地从周围环境摄取氧气和二氧化碳.最为重要的是红细胞薄膜强大的柔韧性和变形能力，使其受流体挤压可轻松通过人体毛细血管(直径2~3μm)或狭窄的脾窦(< 2μm).因此，对红细胞的形变学习和结构分析能有助于了解红细胞膜内骨架中2个重要蛋白(跨膜蛋白和骨架蛋白)复合体的机械形变模式.血液传染性疾病(如：疟疾、流行性乙型脑炎或登革热)，都以红细胞为宿主，方便进行人体全身器官的交叉传染和病毒扩散.对红细胞弹性形变的在线检测，以及红细胞受外界刺激发生的机械响应都将加快生物医学对血液代谢和抑制细菌药物的研发^[1-2].
当前，凭借微流控芯片制备出的微纳尺寸沟道可以有效模拟红细胞在复杂血管中的迁移、黏附、变形和机械疲劳等生理或病理表征^[3].根据微流控技术衍生出的血管芯片为镰刀状细胞贫血病、血液微血管阻塞、疟疾感染红细胞黏附功能，以及遗传性椭圆形红细胞增多症提供丰富的实验模拟环境.在红细胞形变研究方面，Ni等^[4-5]研究了不同酒精浓度下，红细胞在血液流动过程表现的特殊应力应变.他们从理论上揭示了饮酒会破坏红细胞的固有运行规律，低剪切率和高黏度血浆诱导红细胞自发旋转，从而降低红细胞的定向迁移能力，提高了红细胞在血管中发生堵塞的几率.Ye等^[6]和Cetin等^[7]分别采用格子玻尔兹曼方法和有限体积法进行流-固耦合模型的建模，使其应用于红细胞在毛细血管中的形变模拟.Cetin等对流体计算主要采用边心非结构化有限体积法处理任意拉格朗日-欧拉(arbitrary Lagrangian-Eulerian, ALE)方程，处理生物细胞则使用经典伽辽金有限元方法.针对毛细血管直径、红细胞膜厚度，以及红细胞之间的间距都做了较好的理论预测和成因分析.Taraconat等^[8]则对库特细胞计数器内部流动的红细胞形态进行了研究报道.微流通道内红细胞的运动姿态、路径和变形都影响传感信号的识别.所以他们的理论研究不仅能有效提供输出电信号的细胞电学特征表达，还可从测量信号的电场特征值判定细胞形变特征.
国内对红细胞的研究热度也呈现逐年上涨的发展态势，Jiao等^[9]对红细胞在确定性侧向位移微流控芯片分离方面进行了深入研究.他们报道了两种模式(位移和锯齿模式)下，红细胞不同运动轨迹.另外对柱形障碍物间距和错位距离进行了仿真优化.Liu等^[10]使用浸入边界-格子玻尔兹曼方法进行了红细胞在收缩微腔中形态模拟研究.他们结合CPU和GPU计算模式可以有效减少红细胞的模拟时间.同时，张朦等^[11]结合显微镜和微流控芯片技术，开展了红细胞多角度形态学分析的研究工作.他们的实验结果也证实了图像处理算法可以准确识别红细胞数量，后期以红细胞形态为参考可辨识出正常人和糖尿病患者之间的红细胞弹性差异.另外，Wang等^[12]也在红细胞结构形变的仿真中引入神经网络方法，对这种非线性细胞弹性薄膜的体剪切模量进行机器学习.结果表明训练后的四层神经网络对病变细胞具有较好的预测能力.然而目前鲜有****采用COMSOL有限元软件对红细胞在微通道内部的运动轨迹和变形特征作出更为详细的研究.COMSOL有限元分析软件内置的流-固耦合(fluid-structure interaction, FSI)模块可以实现红细胞受流体剪切力作用的动态模拟，该方法不但可以减少繁琐的编程调试工作，还解决了流体场和固体应力场相互耦合产生的多物理场微分方程求解问题.
1 应力介绍微流控芯片构建的血管微环境，内部血液满足不可压缩性流体的假设，即流体密度ρ与时间和空间变化无关.低雷诺数条件下(Re<1300)，腔体内部的流体运动满足层流，相应动量与质量守恒可用Navier-Stokes方程进行描述：

(1)

式中：u为液体流速；p为压强；η为黏度.红细胞薄膜由磷脂双分子、血影蛋白、脂锚定蛋白和整合蛋白等组成，包含着内部细胞基质，满足超弹性结构方程，如下所示：

(2)

式中：ρ_mem为薄膜密度；v为细胞膜的位移矢量；σ为柯西应力张量；▽为矢量微分算子.刚度矩阵C=F^TF，其中F=▽_xv+I描述参考系x方向的变形梯度张量.因此柯西应力张量σ有如下关系:

(3)

式中：J=det(F)是变形梯度张量的行列式；P=?W/?(▽_xv)=FS为第一类Piola-Kirchhoff应力张量，其中S表达如下^{[4, 9]}：

(4)

(5)

式中，S和ε分别代表第二类Piola-Kirchhoff应力张量和Green-Lagrange应变张量.红细胞薄膜满足不可压缩超弹性材料特征，因此应变能密度W服从neo-Hookean方程：

(6)

式中: G₀为初始固体剪切模量，约等于μ₀/h₀(μ₀为初始剪切模量，h₀为膜厚)；κ为初始本体模量；I₁=λ₁²+λ₂²+λ₃²(λ_i为主拉伸系数，i=1, 2, 3)，此外设定的细胞薄膜属于不可压缩弹性材料，因此满足λ₁λ₂λ₃=1.式(6)等号右侧第二项为二级neo-Hookean参量，描述非轴向张力作用λ₁的最大拉伸比.J_el为与温度有关的弹性体积比率，I₁和J_el表达式分别如下：

(7)

(8)

式中，tr表示求取刚度矩阵的迹.固体受力发生的形变量主要与应变能密度W相关.
细胞膜表面的流速u_mem=?v/?t，流体牵引力应该满足边界连续条件，

(9)

式中，n_p和n_f分别表示在参考系和空间系中，细胞界面指向流体的单位法线向量.因此，流体和细胞边界应力张量σ_f可表示为

(10)

正常红细胞呈双凹曲面圆盘形状，具有较好的轴对称特性，绕z轴旋转，柱坐标系(ρ=)下曲线方程如式(11)所示:

(11)

正常红细胞一般具有三种空间结构^[13].为了研究方便，本文仅考虑其中一种常规结构，相应的参数R₀，c₀，c₁和c₂分别为3.91 μm，0.207 161，2.002 558和-1.122 762.根据式(11)，可画出红细胞的空间几何形状，如图 1所示.
图 1(Fig. 1)

图 1 正常红细胞的几何示意图Fig.1 Schematic for a normal erythrocyte geometry

2 仿真结果与形态分析根据上述的流-固方程介绍，本文采用COMSOL有限元软件内置FSI模块进行红细胞运动形变2D模拟仿真.针对流体部分，微通道内部血液的血浆密度和黏度分别为1 060 kg/m³和1.64 mPa·s；细胞内部细胞质的密度和黏度分别为1 050 kg/m³和6 mPa·s.同时考虑超弹体方程，如式(3)~(6)所示.细胞膜涉及弹性参数取值参考文献[4-5, 14]，相应的μ₀和κ分别为6.8 μN/m，2.18×10⁹ N/m²；细胞膜密度与细胞质密度近似相等，为1 050 kg/m³，膜厚h₀为100 nm.基于此，在不同微通道尺寸下，进行黏度和流速等影响因素的红细胞运动建模与仿真分析.
2.1 束腰微通道形态分析首先对红细胞在小尺寸束腰微通道的运动进行模拟计算，本文设定微通道入口的初始流速为0.24 mm/s.两种微通道结构如图 2所示，右侧色度条颜色代表流速大小，主通道高度10 μm，总长50 μm，两种不同束腰的最窄位置高度分别为3 μm和6 μm.白色虚线框代表红细胞在t=0 s时刻的初始位置.当流体驱使细胞自左向右运动，红细胞受流体作用逐渐发生弯曲变形.根据图 2的模拟结果可知，由于流体受到挤压，速度色度条显示出3 μm束腰位置的流速高于6 μm位置的流速.此时红细胞薄膜受到流体剪切力作用发生弯曲，且弯曲方向与流速方向一致.红细胞的双凹形状受流体作用非常容易在此区间发生扭曲.相比球形或椭球形细胞，双凹形结构有利于红细胞在遇到较窄的血液通道时提前实现上下边缘弯曲，方便快速且顺利地通过狭缝.
图 2(Fig. 2)

图 2 不同束腰高度的红细胞流动形变仿真Fig.2 Simulation on RBC deformation under different waist height (a)—束腰高度3 μm; (b)—束腰高度6 μm.

剪切应变率反映了流体与固体在相互作用过程中抵抗剪切变形的能力.微观条件下，它的产生与流-固接触面的黏附力分子有关.如图 3所示，求解红细胞的剪切应变率可知3 μm束腰处具有较大剪切应变率.细胞发生的形变远高于6 μm束腰.同时细胞远离束腰位置，两者的剪切应变率近似，细胞受到大小相同的剪切应力开始恢复到最初进入束腰时的“降落伞”形状.根据第四强度理论得到的von Mises应力，其计算公式如下：
图 3(Fig. 3)

图 3 红细胞流动在微通道所受的剪切应变率Fig.3 Shear rate of RBC flowing inside different micro-channel

(12)

式中, σ_i(i=1, 2, 3)为柯西应力张量的一、二、三主应力分量，反映了细胞膜弹性形变抗屈服破坏的能力.对红细胞边界积分可得细胞行进受到的von Mises应力大小，如图 4所示.3 μm束腰宽度使红细胞在流体和通道双重挤压发生的形变量远高于6 μm束腰宽度.与初始进入通道前的双凹形状细胞相比，3 μm束腰流出的细胞需要一定的时间恢复到前期弯曲状态.
图 4(Fig. 4)

图 4 红细胞流动在微通道的von Mises应力Fig.4 von Mises stress of RBC flowing inside different micro-channel

2.2 黏度形态分析人体器官具有较大差异，检测各自的血液指标也存在很多不一致的参数.其中黏度是影响红细胞正常流动的关键参数之一，血液流变学指出过高的血液黏度不但会对血管造成不可逆转的机械损伤，还会造成动脉粥样硬化和红细胞压积增高.本文因此设定两种不同的黏度值1.64 mPa·s和5.5 mPa·s.血液的初始流速0.24 mm/s和密度1 060 kg/m³保持不变，分别对红细胞的运动形变进行模拟，结果如图 5a~5b所示.红细胞的扭曲方向和形态大体类似，最终形变结果的对比如图 5c所示，可发现黏度较高的血液驱动，红细胞沿流体方向移动的中心位置表现出较大的拉伸形变.另外受黏度力作用的细胞边缘位置有一定的形变收缩.
图 5(Fig. 5)

图 5 不同黏度条件下的红细胞形态变化Fig.5 Dynamic deformation of RBC under different fluidic viscosity (a)—η=1.64 mPa·s; (b)—η=5.5 mPa·s; (c)—不同黏度的红细胞最终形变对比.

红细胞在不同黏度溶液中von Mises应力分析如图 6所示，高黏度环境下的红细胞承受较大的屈服效应，并且随运动距离增大而趋于稳定.此时在较高的液体剪切力作用下，细胞膜的前后方向表现出较大的膨胀形变.从黏度5.5 mPa·s溶液拟可知，细胞处在高黏度血液中运动，它的机械形变一直处于高紧缩的曲胀姿态.如此长时间的运动也容易导致机械疲劳而降低红细胞的生命时间.
图 6(Fig. 6)

图 6 红细胞在不同黏度溶液中von Mises应力Fig.6 von Mises stress of RBC suspending medium with different viscosity

2.3 流速形态分析本节研究不同液体流速对红细胞形态的影响，针对血管网络具有拓扑形状的多通道进行几何建模，初始流速的T型微通道基本结构如图 7所示.通道的宽度与上文一致，仍为10 μm，2个输入端口的流速分别为u₁和u₂.改变流速比(u₂₁=u₂/u₁)可得对应的红细胞运动轨迹和形变如图 7a~7b所示.
图 7(Fig. 7)

图 7 T型微通道内不同流速条件下的红细胞形态Fig.7 Dynamic deformation of RBC under different fluidic velocity in T-shaped micro-channel (a)—初始流速u1=0.24 mm/s，u2=0.24 mm/s; (b)—初始流速u1=0.24 mm/s，u2=0.48 mm/s; (c)—不同流速比的红细胞运动变形对比.

在顶部支路通道中，细胞形状与图 5相似，双凹形结构仍然受流体作用发生弯曲向下流动.但在与水平通道流体交汇处，受流体汇聚挤压造成的推托效应，细胞靠近上表面血管管壁.然而非对称的流体挤压作用使得细胞在之后的形变表现出非对称弯曲的特点.由于远离流体交汇处，细胞靠近血管壁的上端受到的剪切力较小，最快呈现舒展姿态.反之处于转角下端位置的红细胞仍保持“镰刀”勾形姿势.图 7a和图 7b的对比发现细胞长轴逐渐与流体方向平行，呈现出平躺姿态向下游动.将运动前后的细胞形状放于图 7c对比，可知红细胞中间凹陷部分消失，转而在下游通路交汇处发生勾形弯曲.同时图 7b的流速比较大，细胞更倾向平躺流动，而且细胞弯曲量也相对较低.这与下通道位置流体挤压细胞有关.
由此可见，血管通道交汇容易造成细胞发生较大的弯曲变形，同时流速比失调也能引起细胞形变位置的改变.从时间上对细胞的von Mises应力进行分析，如图 8所示.与上述图 7c的形变结论相似，较高的下游流速导致细胞的von Mises应力增大.但是离开混合通道交界处后，在流速比较高的环境中，红细胞承受von Mises应力逐渐减少，其形变恢复较好.由此也证明了红细胞在高速流动的血液环境仍能较好地调整自身的形态特征.另外，对Y型微通道也进行了分析，更改流速比u₂₁=1或2，与图 7a~7b的仿真结果进行对比.细胞仍然发生非对称的弯钩形结构向下游流动，如图 9所示.与T型通道流动相同，在高流速u₂与低流速u₁的液体作用下，红细胞更容易随着液体驱动而呈现伸展形态.
图 8(Fig. 8)

图 8 红细胞在T型微通道内，不同通道流速的von Mises应力曲线Fig.8 von Mises stress curves of RBC calculated on different velocities in T-shaped channel

图 9(Fig. 9)

图 9 Y型通道内不同流速条件下的红细胞形态Fig.9 Shapes of RBC under different fluidic velocity in Y-shaped channel (a)—初始流速u1=0.24 mm/s，u2=0.24 mm/s; (b)—初始流速u1=0.24 mm/s，u2=0.48 mm/s; (c)—不同流速比的红细胞运动变形对比.

对Y型通道的红细胞形变计算可得图 10的von Mises应力曲线.与上述结论相同，细胞在汇聚位置依然承受较强的流体应力而后趋于稳定，最大应力强度低于300 μN/m.本文仿真的细胞形态变化也与文献[5, 7-8]基本相似，该流-固耦合模型为病变红细胞的形状与姿态识别研究提供了一定的理论价值，具有开发速度快、无需编程和建模简单的优点.
图 10(Fig. 10)

图 10 红细胞在Y型微通道内，不同通道流速的von Mises应力曲线Fig.10 von Mises stress curves of RBC calculated on different velocities in Y-shaped channel

3 结论1) 层流驱使红细胞呈降落伞式姿态向下游移动，因红细胞具有较好的弹性形变特征，通过狭窄通道后的红细胞能基本复原为进入前的弯曲形态.较大黏度的流体加速了细胞中心腹部凸起，也引起两端凸形边缘坍缩.
2) 针对非对称流体结构，如T型和Y型微通道，细胞的流体应力分布不均，红细胞表现为镰刀弯曲形态.勾形方向指向另外一侧支路通道，当该支路液体速度大于红细胞通道速度，红细胞勾形较快消失，并伸展为直杆形态.
参考文献

[1]	Czerwinska J, Wolf S M, Mohammadi H, et al. Red blood cell aging during storage, studied using optical tweezers experiment[J]. Cellular and Molecular Bioengineering, 2015, 8(2): 258-266. DOI:10.1007/s12195-015-0380-0
[2]	Mignon T, Mendez S A. Theoretical investigation of the frisbee motion of red blood cells in shear flow[J]. Mathematical Modelling of Natural Phenomena, 2021, 16: 1-28. DOI:10.1051/mmnp/2020054
[3]	戚晓菁, 李学进. 微流控芯片技术在血细胞变形和流动性分析研究中的应用进展[J]. 实验流体力学, 2020, 34(2): 1-10. (Qi Xiao-jing, Li Xue-jin. Research progress on mechanical and flow properties of blood cells in microcirculation using microfluidic devices[J]. Journal of Experiments in Fluid Mechanics, 2020, 34(2): 1-10.)
[4]	Ni A, Cheema T A, Kwak M K, et al. Two-dimensional numerical simulation of the red blood cell floating in a plasma-alcohol solution through stenosis in a microvessel[J]. The Korea-Australia Rheology Journal, 2014, 26(3): 293-301. DOI:10.1007/s13367-014-0033-1
[5]	Ni A, Cheema T A, Park C W. Numerical investigation on the structural characteristics of multiple RBCs in a stenotic microcapillary under plasma-alcohol solution[J]. The Korea-Australia Rheology Journal, 2015, 27(2): 163-171. DOI:10.1007/s13367-015-0016-x
[6]	Ye H L, Shen Z Q, Li Y. Shape-dependent transport of microparticles in blood flow: from margination to adhesion[J]. Journal of Engineering Mechanics, 2019, 145(4): 1-18.
[7]	Cetin A, Sahin M. A monolithic fluid-structure interaction framework applied to red blood cells[J]. International Journal for Numerical Methods in Biomedical Engineering, 2019, 35(2): e3171.1-e3171.24.
[8]	Taraconat P, Gineys J P, Isebe D, et al. Numerical simulation of deformable particles in a coulter counter[J]. International Journal for Numerical Methods in Biomedical Engineering, 2019, 35(11): 1-22.
[9]	Jiao Y Y, He Y, Jiao F. Two-dimensional simulation of motion of red blood cells with deterministic lateral displacement devices[J]. Micromachines, 2019, 10(6): 1-15.
[10]	Liu Z X, Liu H C, Huang D M, et al. The immersed Boundary-lattice Boltzmann method parallel model for fluid-structure interaction on heterogeneous platforms[J]. Mathematical Problems in Engineering, 2020, 2020: 1-13.
[11]	张朦, 孟晓辰, 祝连庆. 基于动态显微成像的红细胞多角度形态学分析方法[J]. 光学技术, 2021, 47(4): 417-421. (Zhang Meng, Meng Xiao-chen, Zhu Lian-qing. Multi-angle morphological analysis of red blood cells based on dynamic microscopic imaging[J]. Optical Technique, 2021, 47(4): 417-421.)
[12]	Wang Y, Sang J, Ao R, et al. Numerical simulation of deformed red blood cell by utilizing neural network approach and finite element analysis[J]. Computer Methods in Biomechanics and Biomedical Engineering, 2020, 23(15): 1190-1200.
[13]	Evans E, Fung Y C. Improved measurements of the erythrocyte geometry[J]. Microvascular Research, 1972, 4(4): 335-347.
[14]	王带领, 谭建平, 程立志, 等. 基于液固耦合的光镊拉伸红细胞变形特性[J]. 中国医学物理学杂志, 2017, 34(4): 404-409. (Wang Dai-ling, Tan Jian-ping, Cheng Li-zhi, et al. Deformation of red blood cells stretching by optical tweezers based on fluid-structure interaction[J]. Chinese Journal of Medical Physics, 2017, 34(4): 404-409.)