基于流固耦合方法的微流控芯片内细胞变形研究

A Study on Cellular Deformation within Microfluidic Chips Utilizing Fluid-Structure Interaction Analysis

摘要: 本文采用流固耦合(Fluid-Structure Interaction, FSI)方法，针对微流控芯片内细胞的变形行为进行了深入研究。为克服传统大批量细胞分析方法在捕捉单个细胞细微差异方面的局限性，我们设计了一种新型微流控芯片，专用于单细胞的捕获。利用COMSOL Multiphysics多物理场仿真软件，我们对该芯片进行了精细建模与仿真，系统探讨了其对不同弹性模量单细胞的捕获效能。通过模拟细胞在芯片中的流动与捕获过程，我们详细分析了细胞表面应力分布、液体流速、捕获时间以及细胞旋转角度等关键参数。研究结果显示，该微流控芯片能够高效且无损地捕获具有不同弹性模量的细胞。这一发现不仅证实了芯片设计在单细胞捕获方面的有效性，还为微流控技术在单细胞分析领域的更广泛应用奠定了坚实的理论基础，并提供了有力的技术支持。

Abstract: In this paper, we investigated the deformation behavior of cells within a microfluidic chip based on the Fluid-Structure Interaction (FSI) method. In order to overcome the limitations of traditional bulk cell analysis methods that are difficult to capture the subtle differences of individual cells, a microfluidic chip was designed for capturing single cells in this study. The microfluidic chip was modeled and simulated using COMSOL Multiphysics multiphysics field simulation software, and the effect of the chip on the capture of single cells with different elastic moduli was discussed in depth. By simulating the cell flow and capture process in the chip, key parameters such as cell surface stress distribution, liquid flow rate, capture time and cell rotation angle were analyzed. The results show that the microfluidic chip can effectively capture cells with different elastic moduli and the cells will not be mechanically damaged during the capture process. This finding not only verifies the effectiveness of the microfluidic chip design in capturing single cells efficiently and nondestructively, but also provides a solid theoretical foundation and technical support for the further application of microfluidic technology in the field of single cell analysis.

1. 引言

随着生物医学技术的飞速发展，单细胞分析已跃居生命科学研究的前沿，它不仅揭露了细胞间的异质性[1] [2]，为我们探索生命的微观奥秘提供了前所未有的精细视野，还极大地推动了疾病的早期诊断[3]、个性化治疗方案的制定以及新药的研发进程[4] [5]。然而，传统的细胞分析方法因其批量处理的局限性，难以精准捕捉单个细胞间的细微差异，从而制约了研究的深入与拓展[6]。因此，研发高效、精准的单细胞捕获与分析技术，已成为推动生命科学研究与临床医学发展的关键所在[7] [8]。

微流控技术，作为微纳尺度下流体操控的精湛艺术，凭借其独特优势在单细胞分析领域脱颖而出。微流控芯片通过精心设计的微通道、微阀和微反应器，实现了对细胞的精准操控、高通量筛选[9]、快速分析[10]及即时检测，极大地提升了实验效率与数据准确性[11]。随着材料科学、微纳制造技术、生物信息学以及计算流体力学的迅猛发展，微流控芯片的设计日趋复杂且功能愈发全面，从基础的细胞捕获、计数，延伸至复杂的细胞培养、分选、裂解、PCR扩增及测序，其应用范围不断拓宽[12]-[14]。

然而，面对个性化医疗需求的日益增长和生物样本的复杂性，微流控芯片设计仍面临核心挑战：如何在保持高通量的同时，进一步提升单细胞的捕获效率、选择性和灵敏度，并减少非特异性吸附和细胞损伤。流固耦合(Fluid-Structure Interaction, FSI)技术的发展，为微流控芯片的设计优化开辟了新路径。FSI是一种描述流体动力学与结构力学之间相互作用的多物理场耦合现象，在生物医学领域应用广泛[15] [16]。在微流控芯片中，FSI分析能够精确模拟流体流动对芯片内部结构的力学作用[17]，以及这些结构变形对流体流动的反馈效应。这对于深入理解流体–细胞–芯片之间的相互作用机制[18]，优化芯片设计，提高捕获效率具有至关重要的意义。例如，通过FSI模拟，我们可以预测在不同流速、粘度、细胞大小和形状条件下，细胞与芯片捕获结构之间的相互作用力[19] [20]。这有助于我们优化捕获结构的设计，使其既能高效捕获目标细胞，又能有效减少非目标细胞的干扰和细胞损伤，从而推动微流控技术在单细胞分析领域的进一步发展。

COMSOL Multiphysics作为一款业界领先的多物理场仿真软件，不仅全面支持流体动力学、传热、电磁场等多种物理场的精确模拟[21]，还特别提供了功能强大的流固耦合分析模块，为微流控芯片的精细化设计与优化提供了得力助手。借助COMSOL，研究人员能够构建出高精度的微流控芯片模型，并模拟在不同操作条件下芯片内部的流体动力学行为，这涵盖了流体的速度分布、压力变化、涡旋形成以及质量传递等多个方面。更为重要的是，COMSOL还充分考虑了芯片结构在流体作用下的变形、应力分布以及细胞的响应情况[22]-[24]，实现了真正的流固耦合分析。这种综合性的分析方法能够深入剖析芯片设计与单细胞捕获效率之间的内在联系，为芯片设计的优化提供科学依据。同时，它还能指导实验验证过程，推动设计的迭代优化，助力微流控技术在单细胞分析领域取得更多突破。

2. 材料与方法

2.1. 有限元分析模型

2.1.1. 控制方程

在本模型中，将液体视为不可压缩流体，其遵循标准的不可压缩流体动量方程，有：

$\rho \frac{\partial u}{\partial t}+\rho \left(u\cdot \nabla \text{}\right)u=-\nabla p+\nabla \cdot \left[\mu \left(\nabla u+{\left(\nabla u\text{}\right)}^{\text{T}}\right)\right]$ (1)

$\rho \nabla \cdot u=0$ (2)

式中：ρ (kg/m³)为流体密度，u (m/s)为流体流速，p (Pa)为压力，μ (N⋅s/m²)为流体动力粘度。等式右边第二项物理意义为粘性应力张量的散度，表示了粘性应力在流体中的分布和变化。

在仿真中，细胞会收到来自于流体的作用力而移动及变形，该部分控制方程为：

$\rho \frac{{\partial }^2{u}_{\text{solid}}}{\partial t^2}=\nabla \cdot {\left(F S\right)}^{\text{T}}+F_V,F=I+\nabla u_{\text{solid}}$ (3)

其中，u_solid (m)表示固体的位移长度，$t\left(s\right)$ 为时间，F_V (N/m³)表示作用在单位体积流体上的力，I为单位矩阵。上式中S (J)为总应变能，通过应变能公式进行计算：

$S=S_{\text{inel}}+\frac{\partial W_s}{\partial ϵ}$ (4)

其中，S_inel (J)是非弹性应变能，W_s (J/m³)为应变能密度，ϵ为损伤参数，用于描述材料的变形和损伤程度。

$F_V=\left[-pI+K\right]\cdot n$ (5)

该式描述了流体在固体边界上的总作用力，其中−pI表示流体压力p通过单位矩阵I施加在固体边界上的力；K是流体的刚度矩阵，反映了流体的弹性特性；n是固体边界的法向量，指向流体一侧。

$U_{tr}=\frac{\partial u_{\text{solid}}}{\partial t}$ (6)

此公式描述了固体在流体作用下的运动状态。其中：U_tr (m/s)为固体移动速度。

2.1.2. 微流控芯片结构设计

该芯片主体为一种H型捕获单元，其由长1560 μm、宽20 μm的捕获长条和弯曲通道构成，捕获长条中心为细胞捕获位置，其左右两侧分别为输出通道和输送通道，捕获位置后方连接的通道宽度为5 μm，每个捕获单元的细胞捕获位置的圆孔直径为12.5 μm。

Figure 1. Diagram of the practical teaching system of automation major

图1. 捕获单元与捕获单元中心捕获位置尺寸图

其中弯曲通道的作用是提供足够的阻力，使得细胞能够首先流向捕获位置。基于等效电路原理，可将图1中的捕获单元转化为等效流阻图，如图2所示。其等效阻值通过流体阻力公式进行计算，在本模型中，可分为沿程阻力损失及局部阻力损失。沿程阻力损失公式为：

$\Delta P_1=\lambda \frac{l}{d}\frac{\rho v^2}{2}，\lambda =\frac{64}{R_e}，R_e=\frac{\rho vd}{\eta }$ (7)

其中，λ为沿程阻力系数，与流体雷诺数Re相关，l (m)为直管长度，d (m)为管路内径，ρ (kg/m³)液体密度，v (m/s)为平均流速，η (Pa·s)为液体运动粘度。

局部损失公式为：

$\Delta P_2=\xi \frac{\rho v^2}{2}$ (8)

其中，ξ为局部阻力系数，与局部障碍物的形状和尺寸有关，在本模型中主要涉及到45˚与90˚“折管”类型的局部阻力系数，其值分别为0.31与1.12。

Figure 2. Equivalent flow resistance diagram

图2. 等效流阻图

2.1.3. 细胞捕获模型建立

在该仿真模型的构建中，本文使用了层流、动网格、固体力学以及流固耦合多物理场。为了提高计算效率，将模型进行了二维简化并对求解域进行了进一步的分割以优化计算。

在捕获通道入口处引入了直径4 μm的圆形细胞模型，并与周围区域一起构成动网格区域[25]；在仿真过程中，细胞会在流体的作用下移动，动网格可确保仿真顺利执行。除了细胞区域外其他所有区域都被设定为层流物理场，并且定义了一个阶跃函数用以提高模型的收敛性[26]；在固体力学物理场中，设置好细胞相关的物理参数[27]；最后添加流–固耦合多物理场，并将两者间的耦合方式设置为全耦合(图3)。

Figure 3. Schematic diagram of domain Settings in the model

图3. 模型内域设置示意图

为了适应复杂的几何形状并保持较高的网格质量[28]，将模型几何划分为三角形网格，特别是捕获位置的网格更加密集，最小网格尺寸仅有0.3 μm，最终模型的总网格数约为56万，平均网格质量达到了0.96 (图4)。

Figure 4. Grid division results

图4. 网格划分结果图

2.2. 材料参数

本文选择甲基纤维素溶液[29]作为流体材料，其运动粘度为15 mPa⋅s，密度为1000 kg/m³；Li等人通过AFM压痕方法测定了人类乳腺上皮细胞的弹性模量位于0.41~1.2 kPa 之间[30]，罗达通过一种改进的挤压变形法测出了前列腺癌细胞PC-3弹性模量参数为0.62~1.89 kPa；LNCaP弹性模量为0.56~1.6 kPa。Esteban等通过改进后的微管吸附方法测得淋巴细胞的泊松比为0.33 ± 0.11 [31]。因此，本文模型中将细胞仿真参数分为四组，前三组细胞模型区分细胞核与细胞质，最后一组不做区分。如表1所示：

Table 1. Cell tissue material parameters

表1. 细胞组织材料参数

组别	细胞核弹性模量[kPa]	细胞核泊松比	细胞质弹性模量[kPa]	细胞质泊松比
I	1	0.3	0.5	0.4
II	1.5	0.3	1	0.4
III	2	0.3	1.5	0.4
O	1.5	0.4	1.5	0.4

3. 结果与分析

为了探究不同组别细胞在被捕获过程中的状态，分别设计了以下四组仿真组仿真实验，实验中芯片入口的液体流速均为200 um/s。

Figure 5. Cell surface stress distribution and flow velocity during capture

图5. 捕获过程中细胞表面应力分布及流速图

图5展示了四组实验中，细胞在微流控芯片不同部位时的表面应力及液体流速分布情况。每组实验占据一行，而每列则对比了不同组细胞在芯片同一位置的结果。从I、II、III组实验结果可以观察到，在细胞未进入捕获位置左侧的狭窄通道前，细胞质区域的表面应力主要分布在130~190 N/m²范围内，而细胞核区域的应力则多集中在250~310 N/m²。当细胞在流体作用下进入芯片左端的狭缝时(如图最后一列所示)，细胞与芯片接触的点出现了显著的应力集中现象，同时细胞核上的应力也相应增大，如第III组实验中，细胞核的应力增大至370~400 N/m²。

图中的流线清晰地反映了芯片内液体的流速分布。在芯片较宽的通道内，中心流速明显低于较窄处的流速。特别是在细胞进入左侧窄缝之前(即细胞捕获位置处)，其周围的流线均匀且流畅；而细胞部分进入窄缝后，流线出现了被阻断的现象。相较于I、II、III组实验，O组细胞在进入左侧狭窄通道前，其表面应力大于前三组细胞质的表面应力，但小于细胞核的表面应力。当细胞部分进入左侧狭窄通道后，细胞表面的最大应力达到390 N/m²，与前三组的结果较为接近。同时，芯片内部的流线也被阻断，并且在细胞上侧出现了环形流线，这表明液体在此处形成了涡旋。这些结果为进一步理解细胞在微流控芯片中的捕获机制提供了重要的依据。

为进一步深入探究细胞表面的应力变化规律，本文测定了整个捕获过程中细胞表面的平均应力，并将结果展示在图6中。图中最引人注目的特点是曲线前后两端应力迅速爬升的两个阶段。在0~0.2秒之间，细胞在流体的推动下开始移动，由于细胞从静止状态开始加速，因此加速度较大，在流体曳力的作用下，细胞表面的平均应力急剧增大。而在1.2秒之后，细胞开始进入左侧的狭缝，受到流体与芯片通道的共同作用，细胞受到挤压而变形，导致细胞表面应力再次增大。

Figure 6. Nucleus-plasma surface stress variation diagram

图6. 细胞核–质表面应力变化图

值得注意的是，在0.2~1秒之间，细胞核表面的应力始终远大于细胞质表面的应力，且前者几乎为后者的两倍。这是因为细胞质的弹性模量相对较小，在相同的外力作用下，细胞质会产生较大的形变来分散和减小应力。通过对比O组与II组细胞核的平均应力分布，我们发现，在弹性模量相同的情况下，泊松比越大，平均应力也越大，这意味着泊松比大的细胞更容易出现应力集中现象。此外，对比O组与III组细胞质的平均应力，我们发现O组的细胞质平均应力略高于III组，尽管两组的细胞质弹性模量和泊松比数值均相同。我们推测，这可能是由于III组细胞核的泊松比较小，在受到相同外力作用时，III组的细胞核变形量小于O组，因此III组的细胞质变形量也应当小于O组。同时，考虑到III组细胞质的弹性模量高于O组，这解释了为何III组的细胞质平均应力会高于O组的现象。这一推测与王玉娇等人在宫颈癌细胞力学特性研究中的发现相一致[32]。

图中的三个彩色条带分别代表了细胞进入“分岔口”、“通道折点”和“细胞捕获位置”的时间区间。我们可以发现，细胞在通过分岔口和通道折点时，细胞核与细胞质表面的平均应力都会产生较为明显的波动。这是因为细胞的运动方向发生了较大改变，导致细胞表面各微小单元之间产生了较大的应力。而在细胞进入捕获位置时，由于运动角度未发生变化，因此应力没有产生明显变化。

为更精确地研究细胞在芯片中的运动状态，本文将细胞的运动分解为x和y两个方向进行单独分析，结果如图7和图8所示。观察这四组细胞的运动状态，我们发现它们具有高度一致性。在进入分岔口之前，细胞在x方向的位移基本保持为0，而在y方向上则匀速增大；这表明细胞在这一阶段主要是沿着y方向(即通道的主方向)移动。进入分岔口后，细胞在x方向的位移开始逐渐增大，同时y方向上的位移也继续匀速增大，但增长斜率相较于之前有所减小。这说明细胞在分岔口处开始发生横向(x方向)的偏移，同时仍然保持着沿y方向的移动，但速度有所放缓。经过折点后，细胞在x方向的位移增速明显变大，而y方向的位移则保持不变。这表明细胞在折点处发生了较大的横向偏移，而沿通道主方向(y方向)的移动则暂时停止或减缓。

Figure 7. Cell displacement diagram (x-axis)

图7. 细胞位移图(x轴)

图9展示了细胞在通道内运动的平均速度图像。从图中可以更明显地看出，细胞在进入分叉口之前的速度大约在200 um/s左右。在经过分岔口时，细胞经历了明显的减速，随后速度又提升至约155 um/s。有趣的是，细胞在经过折点时，速度发生了一阵较小的波动。我们推测这可能是由于流体在流过折点时，方向在短时间内发生了较大的变化，从而影响了细胞的速度。在1秒之后，细胞进入捕获位置。捕获位置是一个直径为19 um的圆形区域，远大于通道宽度。液体在进入捕获位置后速度会迅速放缓，因此细胞的移动速度也同样降低，最低约为100 um/s。随后，细胞移动速度又开始增大，这是因为捕获位置左端的通道宽度仅为5 um，液体在进入这一较窄通道后流速会高于在8 um通道内的流速，因此细胞的移动速度会出现增大的现象，并且速度将会高于200 um/s。

Figure 8. Cell displacement diagram (y-axis)

图8. 细胞位移图(y轴)

Figure 9. Variation of cell migration velocity

图9. 细胞移速变化图

李文姣等人在其研究中发现，细胞在靠近壁面移动时会出现“跳跃”现象[33]，细胞跳跃高度在0.5~0.9 um之间。然而，在本文的模型中，细胞同样在近壁面运动，却并未出现该跳跃现象。仔细对比后发现，该作者文中细胞尺寸与流体腔尺寸之比为0.12，而本文中的比例约为0.94。我们推测，由于空间过小，导致无法观察到跳跃现象。在本文的模型中，细胞与通道壁面的间隙较小，可能限制了细胞的跳跃行为。

在仿真实验中，我们观察到细胞在流体作用下不仅会在x和y方向上发生位移，同时还会发生旋转。如图10所示，展示了细胞在被捕获过程中的角度变化结果。O组、I~III组的细胞角度变化范围分别为：−1.8˚~−0.8˚、−1.8˚~−0.6˚、−1.5˚~−0.5˚、−2˚~−0.5˚。其中，第III组的细胞旋转角度变化最大，达到了1.5˚。这是因为第III组细胞的细胞质与细胞核的弹性模量都是最大的，这意味着细胞更难以变形，因此在相同的作用力下会产生更为剧烈的旋转。

Figure 10. Variation of cell rotation angle

图10. 细胞旋转角度变化图

四组细胞在0.2秒至0.4秒期间均呈现出旋转角度的下降趋势。这是因为细胞在进入分岔口时，其左侧的流体速度小于右侧，从而造成了顺时针旋转。当细胞完全通过分岔口之后，其两侧液体流速变得相同，此时在弹性力的作用下，细胞角度开始回转。在这个阶段中，O组细胞的回转角度最大，甚至超过了0˚。同时，我们注意到I、II、III组细胞的回转角度大小与预期不符。为了更深入地探究这一现象，本文计算了四组细胞的加权平均模量，分别为：1.5 kPa、0.63 kPa、1.13 kPa、1.63 kPa。我们发现，第I组的加权弹性模量最小，但其回转角度却大于第II组；而第III组的加权弹性模量最大，其回转角度却小于O组。这些结果似乎表明，细胞的旋转角度不仅与其弹性模量相关，还可能受到其他因素的影响。

因此，我们可以推断，细胞的旋转行为是一个复杂的物理过程，不仅受到细胞自身力学特性的影响，还可能受到流体动力学特性、细胞与芯片壁面的相互作用以及细胞形状和大小等多种因素的共同影响。未来的研究需要更全面地考虑这些因素，以更准确地预测和控制细胞在微流控芯片中的旋转行为。

4. 讨论

本研究通过精心设计的四组仿真实验，对细胞在微流控芯片中被捕获过程中的应力分布、运动状态及旋转角度进行了全面而深入的探究，得出了一系列富有意义的结论。这些结论不仅深化了我们对细胞在微流控环境中力学行为的理解，也为微流控芯片在细胞生物学研究中的应用提供了重要的理论支持和实践指导。

首先，实验结果清晰地揭示了细胞在通过微流控芯片不同部位时应力分布的变化规律。在未进入捕获位置左侧的狭窄通道前，细胞核区域的应力普遍高于细胞质区域，这反映了细胞核相对坚硬的物理特性。当细胞进入狭窄通道时，应力集中现象显著增强，特别是细胞与芯片接触的点以及细胞核区域，这进一步验证了微流控芯片在细胞捕获和操控中的有效性。通过对细胞表面平均应力的进一步分析，我们发现细胞核表面的应力始终远高于细胞质表面，这主要归因于细胞质和细胞核在弹性模量上的差异，细胞质具有较大的弹性，能够有效分散应力，而细胞核则相对坚硬，更易出现应力集中。这一发现不仅为我们提供了关于细胞在微流控环境中应力分布变化的宝贵信息，也为我们理解细胞在微流控芯片中的力学行为提供了新的视角。

其次，泊松比对细胞应力分布的显著影响是本研究的一个重要发现。泊松比越大，细胞在相同外力作用下的应力集中现象越明显。这一发现不仅为深入探究细胞力学特性提供了新的线索，也提示我们在设计微流控芯片时，需要充分考虑细胞的力学特性，以确保芯片能够实现对细胞的精确操控。

此外，本研究还证实了捕获过程中细胞表面应力未超过400 Pa，远低于H-22型肝癌细胞可承受的最大微流剪应力阈值697 Pa [34]。这一发现表明，本文所设计的芯片不会对细胞造成机械损伤，为微流控芯片在细胞生物学研究中的应用提供了安全保障。这一结论不仅增强了我们对微流控芯片安全性的信心，也为我们进一步探索和优化微流控芯片的设计提供了重要的参考依据。

在捕获过程中，细胞速度的变化反映了流体在微流控芯片中的复杂流动特性。细胞在进入分岔口时经历减速又回升，以及经过折点时速度出现波动。这些发现不仅有助于我们优化微流控芯片的设计，提高细胞捕获效率，也为我们理解流体在微尺度下的流动行为提供了新的视角。这一部分的讨论加深了我们对流体在微流控芯片中流动特性的理解，为优化芯片设计和提高捕获效率提供了重要的科学依据。

最后，实验还揭示了细胞在流体作用下不仅发生位移，还会产生旋转。不同组别的细胞在旋转角度上存在较大差异，这主要与细胞的弹性模量有关。然而，我们也发现细胞的旋转角度不仅与其弹性模量相关，还可能受到流体的流动特性、细胞与芯片壁面的相互作用等其他因素的影响。这一发现为深入理解细胞在微流控芯片中的动态行为提供了新的视角，有助于我们更全面地掌握细胞在微流控环境中的运动规律。这一部分的讨论拓宽了我们对细胞在微流控芯片中动态行为的认识，为未来的研究提供了新的研究方向。

综上所述，本研究通过仿真实验对细胞在微流控芯片中的应力分布、运动状态及旋转角度进行了全面而深入的探究。这些发现不仅为我们提供了关于细胞在微流控环境中力学行为的重要信息，也为微流控芯片在细胞生物学研究中的应用提供了有力的理论支持和实践指导。未来，我们将继续深入探索细胞在微流控环境中的力学行为，为微流控芯片的设计和优化提供更多的科学依据。

5. 结论

本文基于细胞的表面应力分布、移动速度以及旋转角度这两个核心维度，运用多物理场仿真实验方法，对所设计的微流控芯片在单细胞捕获方面的可行性和有效性进行了深入验证。主要研究结论概述如下：

1) 当细胞通过微流控芯片中的狭窄通道时，细胞核区域展现出显著的应力集中现象。具体而言，细胞核表面的应力水平(约300 N/m²)普遍高于细胞质表面(约150 N/m²)，这一差异主要归因于细胞质与细胞核在弹性模量上的不同。此外，泊松比也对细胞的应力分布产生重要影响，泊松比增大时，细胞更易出现应力集中现象。

2) 细胞在微流控芯片中的运动状态深受流体流动特性的影响。特别是在分岔口和通道折点等关键位置，细胞的运动方向和速度发生显著变化，进而导致细胞表面应力产生明显波动。这些发现对于优化微流控芯片的设计和提高细胞捕获效率具有至关重要的意义。

3) 细胞的旋转角度与其弹性模量存在紧密关联，弹性模量较大的细胞在相同作用力下会呈现出更为剧烈的旋转。然而，细胞的旋转角度还受到其他多种因素的共同影响。在本文所设计的微流控芯片中，细胞旋转角度的最大值仅为2˚，这表明在整个捕获过程中，细胞保持了相当稳定的状态。

NOTES

^*通讯作者。