精子运动行为机理的CFD研究进展

Research Progress of CFD on the Mechanism of Sperm Motility

摘要: 精子的使命是寻找卵子将遗传信息传给下一代，其运动规律一直以来都是生殖科学的重点研究对象。由于实验方法的局限性，计算流体力学逐渐在精子运动特性模拟研究方面发挥重要作用。目前已有研究采用有限元、边界元等流体建模方法结合阻力理论、浸入边界等流固耦合方法研究精子的鞭毛运动、导航特性、相互作用等机理。本文对精子运动机理研究现状以及相关CFD模拟研究进行了总结论述，介绍了研究中各类方法及相关应用，并在此基础上对未来精子CFD模拟研究的发展进行了展望。

Abstract: The mission of sperm is to find the egg and pass on the genetic information to the next generation. Its movement rules have always been the focus of reproductive science. Due to the limitations of experimental methods, computational fluid dynamics (CFD) has gradually played an important role in the simulation of sperm motility. At present, fluid modeling methods such as finite element method, boundary element method and so on, combined with fluid structure coupling method which includes resistive force theory and immersion boundary etc., have been used to study the mechanism of sperm flagellum movement, navigation characteristics and interaction. The status of sperm motility mechanism research and relevant CFD simulation research are summarized and discussed in this paper, which will introduce various methods and related applications. On this ba-sis, the future development of CFD simulation research on spermatozoon is prospected.

1. 引言

目前，人们发现某些生物、化学、物理因素能够影响精子的运动方向，如环境流体流动方向 [1] [2] [3]、温度 [4]、某些化学物质浓度 [5] [6]、壁面作用 [7] [8] 以及精子间竞争合作 [9] [10] 等，这些因素为精子成功游向卵子提供了线索。精子运动机理在生殖医学、微机器人等领域有重要应用，如设计微流控芯片进行体外精子筛选 [11] [12] [13]、设计人工精子 [14] [15] 并应用于靶向药递送 [16]。由于精子身长在微米尺度，微观结构力学作用复杂，鞭毛摆动频率高，实验方法只能对精子宏观运动规律进行研究，其微观运动机理一般需要理论推演的方法来进行探索。近年来随着计算机技术发展，计算流体力学(computational fluid dynamics, CFD)方法在航空航天、水力机械、生物游动方面的研究对发展起到了十分重要的推动作用。CFD方法被广泛应用于精子的推进效率 [17] [18]、近壁游动 [19] - [26]、群体游动 [27]、趋化性 [28] [29]、趋流性 [30] [31]、非牛顿流体中游动 [32] [33] [34] [35] 等领域的研究。本文重点就精子数值建模方法与精子运动机理CFD研究两方面进行了总结和评述，本文的工作可为精子运动的CFD研究提供方法及应用层面的参考。

2. 精子研究的数值建模方法

精子游动的雷诺数Re为0.01左右 [36]，粘性力远大于惯性力，N-S方程中对流项可被忽略使方程线性化为Stokes方程，这种微尺度低雷诺数游动通常以求解Stokes方程来求解流体流动。本节针对当前精子CFD研究领域使用的流场求解方法及精子建模方法进行简要概述，并讨论其各自优缺点及所能解决的相应问题。

2.1. 精子研究的流场求解方法

2.1.1. 有限元法

有限元法是一种求解偏微分方程近似解的纯数值方法。其将问题区域分解为相互连接的小子域，然后将子域上的简单方程联系起来使用近似解去推导求解域满足的条件，最终得到问题的解。有限元法可用于研究精子鞭毛摆动的稳定性 [37]、头部形状对游动影响 [38] 等问题，相应的还有COMSOL有限元软件研究精子与圆形细胞在通道中分离状况的例子 [39]。有限元法应用领域较广，由于其有比较成熟的建模软件可以使用，所以在精子研究时建模方便简单易操作。但其适用于精度要求不高、精子数量少、计算量不高的场合，对于无限域内建模及多精子模拟不是很理想，且精度限制了精子游动中的力学机理的研究分析，所以目前有限元法在精子游动模拟中不是特别广泛。近年来有限元法多与实际工程问题结合，探究宏观运动规律、参数特性等，如有学者对声波驱动的精子机器人建模研究，从而参考模拟结果设计材料、声场等参数 [40]。

2.1.2. 边界元法

Youngren等提出了一种新的“边界积分方程”方法既解决了流场中形状不规则的stokes流问题，又解决了有限元等纯数值解法难以处理较大区域(如半无限域)计算量问题 [41]，以stokes流积分方程的解为基础的方法大多可以称为边界元法。近年来许多精子CFD研究应用了边界元法，如研究精子近壁游动规律 [21] [23] [25] [26]、趋流性 [30]、游动效率与流体粘弹性关系 [34] 等。不同于有限元法在整个连续区域内划分单元，边界元法仅需将求解区域边界离散划分为单元，并用微分算子解析解作为边界上积分方程的核函数，所以称为半数值半解析方法。相较于有限元法，离散方便、精度高、便于处理无限域等优点，使得边界元法成为近年来精子CFD研究多喜欢采用的方法，多用于研究单精子的水动力学特征。但其也有许多不足，如需要计算相当数量的奇异积分、应用领域不如有限元广、计算量依然很大难以大规模计算等。

2.1.3. 格子玻尔兹曼方法

有限元、边界元等传统方法是计算流体宏观物理量(质量、动量、能量等)的宏观连续性方法。而格子玻尔兹曼方法(Lattice Boltzmann Method, LBM)作为一种介于微观粒子与宏观参量之间的数值模拟方法，其将微观粒子运动的统计平均结果来描述流体运动，避免了直接求解复杂的N-S方程，流体的宏观参量也可由微观粒子统计平均值得到，目前已有研究采用这种方法模拟纤毛细丝在流场中运动 [42] [43] [44]。相较于有限元、边界元等传统方法，LBM能够避免网格划分复杂性建立简易网格，结合浸入边界法可以更好地处理流固耦合中的细节问题。低雷诺数下微流场的模拟方面LBM有着巨大优势，能够深入揭示水动力学机理，如Liu等详细分析了精子的近壁面吸引力产生的水动力学机理 [22]。LBM目前在流体力学等领域得到了广泛关注，具有算法简单易懂、编程容易、可大规模并行式计算等优点。目前对于较低雷诺数下的精子游动模拟有着天然优势与适用性，但对于高马赫数、高雷诺数的气动力学仍然较难使用。

2.1.4. 多粒子碰撞动力学

以上有限元等方法属于有网格方法，均是基于一套已经预设好的网格来划分流场。而多粒子碰撞动力学(Multi particle collision, MPC)属于一种无网格方法，将流体建模为一组任意分布粒子而不再预设一套网格系统，通过粒子间碰撞迁移模拟流体运动，在处理复杂边界作用结构时有很大优势 [45]。鉴于这种优势，Elgeti等得以研究多精子绑定聚集时复杂的多鞭毛边界相互作用，能够很好地揭示其中的水动力学机理 [46]，其之后又探究了三维情况下精子近壁吸引现象 [19] 及精子在直角通道转向时壁面吸引机理 [47]。MPC方法有些类似于LBM同样包含粒子迁移与碰撞的计算步骤，相较于LBM不需要划分网格，能够更好地处理不规则复杂边界。但是诸如MPC等基于粒子的无网格方法，其需要大量的粒子来保证模拟的求解精度与结构的精细度，这无疑大大增加了计算量与成本。目前其在精子研究方面无法满足更多精子的计算模拟，但未来计算集群算力提高、程序并行化优化更好的情况下，MPC等无网格方法将在计算流体领域有更好前景，其应用于精子研究也成为可能。

2.2. 精子建模中的流固耦合问题

目前在精子模型的构建方面，通常将精子视为一根连续分布离散点连接构成的摆动细丝，指定各点之间相对运动规律来产生鞭毛波动效果，处理这种波动与流体的相互作用的方法，称之为流固耦合。常见的流固耦合方法有阻力理论(Resistive Force Theory, RFT)，斯托克斯子方程(Stokeslet)与浸入边界法(Immersed Boundary Method, IBM)。

浸入边界法(Immersed Boundary, IB)最早于1972年提出并用于模拟血液在可收缩心脏瓣膜中的运动 [48]。这种方法将浸没在流体中的复杂边界看作以拉格朗日网格表示的无质量边界，将其上每个点的力模化为体积力，以正则化的插值函数实现流体与边界的相互作用。目前已有研究应用该方法模拟精子游动 [22] [24] [49]，由于IB能够很好处理复杂边界，所以能很好地研究精子游动的力学机理、多精子问题，是较新、可行性较高的方法。并且其具有计算效率高的优势，在精子三维模拟、可压缩流体及非牛顿流体模拟方面，未来将有着十足的发展空间，而提高收敛性与稳定性也是IB需要解决的问题。阻力理论(RFT)通过计算固体微元和本地流体微元同向运动速度之差，然后加权得到流体对本微元的粘性阻力，是由细长体理论发展而来的代数近似方法 [38]。使用RFT可以研究剪切流中 [50] 与壁面附近运动规律 [26]、精子转向机制关系等 [51]，但其将精子当作鞭毛与头部的线性叠加，所以RFT多适用于研究鞭毛自身特性机理及头尾相互耦合作用可以忽略的精子游动 [52] [53] [54]。也正是因为这种相互作用的忽略，使得RFT并不能完美、准确地探究精子力学机理，如其将粘性力作用简化为当地速度处理，无法具体分析许多现象的力学本质。此外，当前RFT仅限于模拟二维情况下鞭毛波动，无法对三维情况下螺旋形摆动进行建模分析 [38]，所以泛化性与发展前景较差。对于正则化斯托克斯子方法(RSM)，公式(1)表示stokes流场在精子表面附近的边界积分方程，在处理细长鞭毛时可以忽略后一项的双层积分，只保留前一项单层积分。RSM将中心节点力平滑到表面导出压力速度表达式，再引入截止函数对当地流体速度正则化来解决 $J\left(x,y\right)$ 导致的积分奇异性问题 [55]。

$V\left(x\right)=-\frac{1}{8\text{π}\mu }{\displaystyle \int J\left(x,y\right)\cdot q\left(y\right)\text{d}s}\left(y\right)+\frac{1}{8\text{π}}{\displaystyle \int v\left(y\right)}\cdot T\left(x,y\right)\cdot n\left(y\right)\text{d}s\left(y\right)$ (1)

此方法大多与边界元法结合，且其可以很好地适用于高粘度微流动环境下的柔性边界，故在精子研究方面也有许多应用，尤其在无界流场中的精子游动，如研究精子游动速度与效率 [36]、建立趋化因子耦合模型研究趋化性 [28] [29]、精子趋流性 [31] 等。相较于RFT方法，RSM能处理头尾相互耦合作用，且可以用于三维鞭毛螺旋摆动的建模，未来也有望与更多的流场建模方法结合来应用于精子研究。

3. CFD在精子运动机理中应用

3.1. 研究精子鞭毛运动

正常人类精子体长约60 μm，其中鞭毛长约55 μm，能在流体中产生有节律的波动推动自身迁移，相关研究表明Ca²⁺等化学物质会影响鞭毛摆动与推进模式 [56]。由于鞭毛是精子推力的主要来源，且其力学模型容易建立、不必考虑复杂的头与鞭毛耦合关系、计算简单易实施，所以鞭毛的力学致动机理、数学模型、运动特性等是当下研究的重点。力学致动机理方面，Olson等将Ca²⁺及刚度耦合到鞭毛模型中并探究其对精子运动轨迹的影响 [29]。而刚度变化会导致鞭毛摆动波形的变化，从而使得精子游动活性产生巨大差异 [52]，鞭毛波形中的二次谐波的变化又会影响精子的转向 [51]。所以，如何准确复现鞭毛摆动波形并构建合适、稳定的鞭毛模型依然是有待解决的问题。尽管已有研究对不同数学模型下鞭毛摆动稳定性进行探究 [37]，但摆动中的力学机理是否能够准确反映鞭毛实际运动效果、以及相关趋化因子的耦合作用如何影响其摆动波形也需要深入研究。且当前研究将精子简化为鞭毛忽略了头部的影响，不能准确揭示运动机理，如Rorai等研究认为忽略头部影响导致模拟结果不准确 [38]，所以将精子简化为鞭毛研究其运动特性虽然有许多优点(如能够广泛的应用于鞭毛类生物的群集效应研究)，但能否准确揭示其运动依然有待商榷。此外不同于二维空间波动模型，三维空间螺旋状摆动的鞭毛模型需要进一步构建。

3.2. 研究精子导航机理

精子导航机制包括趋化性、趋温性、趋流性、近壁作用几个方面，其中趋化性指精子可以感受某些趋化物分子并与之结合，使其朝向或远离趋化物源运动。有研究发现Ca²⁺梯度浓度可使获能精子转向 [5]、孕酮可使精子游动加速 [6]。当前CFD对趋化性的研究方法为将Ca²⁺的化学影响耦合到鞭毛摆动的力学模型中来探究精子运动规律 [28] [29]。趋温性是精子能够感受周围环境温度并朝温度更高地方游动的特性。Bahat等发现了这种特性，之后测出了人类精子趋温范围与敏感度 [4]。趋流性则是精子感知周围流体的流动状态来调节自身运动的一种特性，精子具有逆流而上的负趋流性 [2]，且已有CFD研究进行了验证 [31]。趋流性是精子游动的主要导航机制 [1]，相关实验 [3] 及CFD研究 [30] 证实这种特性是一种被动物理行为。精子的导航机理探究是目前一个重要的研究方向，目前关于趋化性、趋流性的相关研究较少，且多集中于验证精子游动实验结果的测定规律，并未很好地揭示其力学机理及相关现象成因。趋温性机理发现较晚，且由于温度会改变流体密度、粘度、可压缩性等，流体建模过程中方程求解相对困难，对于采用CFD研究趋温性规律、机理是一个不小的挑战。而未来在相应模型理论进一步成熟的基础上，单个精子及精子群体在温度梯度中的运动规律、相互作用机理也是一个新的研究发展方向。精子在有界区域游动过程中倾向于靠近壁面游动，进而使得精子能够在生殖道内依靠壁面引导游动方向，也可归为一种导航机制 [8]。近壁现象最早在1963年Rothschild等发现，并认为可能与水动力学有关 [57]。目前相关CFD研究主要集中于近壁水动力学机理 [19] [22]、游动稳定性 [20] [21] [26]、近壁游动效率 [24] [25] 及阶跃型边界附近游动规律 [23] 等。时下对于精子近壁游动规律CFD研究是精子机理探索的热点，但仍有进一步探索空间；且不同于二维情况下波浪式摆动，三维情况下鞭毛实际为螺旋摆动，虽然Elgeti使用三维模型对精子壁面附近游动规律进行了探究 [19]，但未来依然需要进一步研究三维情况下壁面附近游动规律。

3.3. 研究精子在非牛顿流体中运动

非牛顿流体也称为粘弹性流体，其粘度与剪切速率不成正比。一般CFD研究将流动环境简化为牛顿流体，简化产生的水动力学差异目前缺乏深入研究与认知。如Ishimoto等模拟研究认为精子近壁游动时的水动力学特征与牛顿流体中没有明显差别 [32]，但有研究发现非牛顿流体会促进精子间吸引与聚集 [33] [35]，这种聚集是否因为不同性质的流体中水动力差异导致的依然有待研究。此外精子游动速度及效率也会受到非牛顿流体的粘弹性系数的影响 [34]，但依然未能揭示其中的力学机理。随着更确切的鞭毛力学模型的完善，导航机制与运动特性的探明，精子运动将进一步结合非牛顿流体特性来探究更真实环境下的游动规律，从而为人工精子在人体环境中游动进行靶向药递送提供更好的理论支撑。

3.4. 研究多精子相互作用与群集

精子在寻找卵子的过程中面临着竞争，为了更快前进需要合作形成群体，研究发现鹿鼠精子 [9] 与鸟精子 [10] 中的合作会加速游动。当下CFD研究发现水动力能够使鞭毛摆动同步最终聚集成群体，如图1 [46]。

同样Schoeller等发现鞭毛厚薄、频率会使得鞭毛摆动同步从而聚集成群 [27]，不同的是他们将精子模型放到了“伪”三维空间，精子鞭毛摆动及整体运动处于一个固定的空间平面内。目前多精子作用研究主要关心诸如效率、速度、功率等运动参量 [17] [18] 或聚集成因与绑定现象 [27] [58]，揭示精子间相互作用力学机理的研究较少。而未来随着计算机算力的提高、非牛顿流体数理模型的完善，非牛顿流体中多精子相互作用的水动力机理、粘弹性与精子聚集间的力学作用关系，有望得以揭示，也将成为较为前沿的研究方向。

Figure 1. Synchronization of sperm attraction and binding and the phenomenon of sperm aggregation and binding

图1. 精子互相吸引绑定同步与多精子聚集绑定现象

4. 结语

对于实验无法量化的特性，精子CFD研究能起到较好效果并揭示运动机理，对精子筛选、人工精子设计等提供理论依据，起到指导作用。虽然在精子运动机理方面已经有CFD模拟进行研究，但依然存在现有精子模型难以准确复现实验结果、鞭毛摆动模型研究相对滞后、导航机制中水动力学机理有待深入探究、精子间相互作用及粘弹性流体中运动特性相关研究较少的问题。未来精子CFD研究需要构建合适模型解决这些问题，并将其综合起来总结研究、论证、对比，相应研究结果也将辅助设计微流控芯片筛选精子、制造人工精子进行靶向药递送等。总之，精子运动机理的CFD研究将为我们科学深入地了解生命诞生过程、鞭毛类微生物游动群集机理等提供巨大帮助。

基金项目

国家自然科学基金资助(编号：81771935)。

NOTES

^*通讯作者。