1.基于应力失效模型的细胞刺入速度规划方法，其特征在于，包括：
步骤1，建立细胞刺入仿真环境；
步骤2，建立刺入过程中细胞机械损伤的应力失效模型；
步骤3，检测细胞刺入过程中的应力失效比例；
步骤4，基于应力失效模型评估不同速度的细胞刺入过程，对刺入速度进行规划。
2.根据权利要求1所述的方法，其特征在于，所述步骤1具体为：
基于物质点法，将细胞建模为由粘弹性材料点组成的外半径为rpe的球壳包裹着由粘塑性和粘性流体材料点混合组成的半径为rcy的球体，球壳和内球分别表示透明带与细胞质；
吸持针和注射针设置为刚体，形状为空心圆柱体结构；
将吸持针、注射针接触的细胞材料点的相对速度在指向吸持针、注射针表面的法向方向上被投影为0；
定义作为细胞质材料点受到的摩擦力为：
Ffriction＝α·vrelative·m
其中，m为材料点质量，α为摩擦力的比例系数，vrelative为吸持针、注射针与针管周围的细胞质粒子的相对速度。
3.根据权利要求2所述的方法，其特征在于，在吸持针、注射针的针管内和管口附近区域设置吹吸力场：
针管口处的吹吸力场呈球面圆锥状，球面圆锥的顶点位于针管口内的轴向中心线上，球面圆锥的锥面与针管口内壁相交，方向为其中d为球面圆锥内任一点到球面圆锥顶点的向量；
在针管内部区域，吹吸力场方向设置为平行于管壁的方向向量；
定义吹吸力场Fpressure为：
Fpressure＝β·npressure·m
其中，m为材料点质量，β为设定的吹吸力大小，npressure为针管内外的吹吸力场的方向。
4.根据权利要求1所述的方法，其特征在于，步骤2具体为：
定义应力失效准则为，在一次细胞刺入操作过程中，一旦细胞材料点p在任意时刻t承受的应力P(Fp)超过阈值Pfail，该粒子失效且不再恢复：
if:P(Fp)＞Pfail,failp＝1
其中P(Fp)为粒子的第一Piola-Kirchhoff应力张量，使用共旋弹性模型计算：

其中，Fp为粒子的形变梯度，Ψ(Fp)为粒子的弹性势能密度，σpi是Fp的奇异值，μ和λ为粒子的拉梅常数，Jp是粒子体积比，Jp＝det(Fp)，det表示计算矩阵的行列式，为避免数值不稳定，Jp使用下式迭代更新：其中，为粒子的局部速度梯度，tr表示计算矩阵的迹；Rp为计算旋转矩阵，由Rp＝Up*VpT计算得到，其中Up和Vp是由每个粒子的形变梯度Fp进行SVD奇异值分解后的Up,sigp,Vp中得到。
5.根据权利要求1所述的方法，其特征在于，应力失效比例计算公式为：

其中，sum(pfail)为应力失效材料点的总数，sum(p)为所有材料点的总数，使用两者的比值反映细胞中损伤区域的相对大小。
6.根据权利要求1所述的方法，其特征在于，所述步骤4具体包括：
步骤4.1，计算单次的细胞应力失效比例；
初始化细胞、吸持针和注射针；设置细胞质半径和透明带半径，两根微针各自的内外径；将细胞置于吸持针与注射针之间；
吸持针设置吸力场，吸引并固定细胞位置，使得吸持针中心、细胞中心以及注射针中心位于同一条水平线上；在吸持稳定后，注射针以一个固定大小的速度沿水平线朝向细胞中心的方向运动，停留至稳定后，以相同的速度退回初始点，完成细胞刺入仿真；
每一个时间步计算当前每个材料点所受应力，判断材料点是否失效，统计所有材料点的应力失效比例；
步骤4.2，计算不同速度刺入下的应力失效比例；
仿真环境初始设置相同，在注射针接触细胞后的刺入深度相同的情况下，改变注射针刺入速度，输出刺入完成后的应力失效比例随着刺入速度的变化情况；
步骤4.3，基于不同速度刺入下的应力失效比例对刺入速度进行规划；
由于细胞损伤和细胞活性呈负相关性，基于应力失效模型使用不同速度刺入下的应力失效比例评估细胞活性；选择细胞活性最高时相对应的速度进行刺入，从而最小化细胞刺入的机械损伤，提高细胞的后续发育率。

技术领域
[0001]本发明属于细胞级别的显微操作技术领域，具体涉及一种基于应力失效模型的细胞刺入速度规划方法。
背景技术
[0002]细胞刺入技术是生物医学工程中的一项关键环节，被广泛应用于细胞核移植、卵胞浆内精子注射、胚胎微注射等复杂操作。该技术常涉及使用两根微针，一根较粗的吸持针通过负压吸持细胞用于固定，另一根较细的注射针刺入细胞并进行后续的吸取或注射。然而，这一刺入过程不可避免地导致细胞内部被挤压受力，造成的机械损伤直接影响操作后细胞的发育，甚至可能导致细胞死亡。
[0003]细胞刺入操作的成功率和刺入后细胞的发育率直接影响包含该操作的复杂生物细胞实验的成功率。尽管机器人化微操作技术的发展提高了细胞刺入操作的成功率，但被刺入后细胞的发育率较手动方式并没有得到显著提高，因此，细胞刺入技术仍需在减小细胞损伤、提高细胞后续发育率方面进一步优化。
[0004]生物研究者的发现表明，细胞的刺入速度影响细胞的发育率。当细胞受到外部力或应力作用时，其内部的结构和生物分子可能会受到影响，导致细胞受到损伤。然而，当前仍缺少关键指标衡量细胞刺入速度与细胞机械损伤之间的关系。因此，迫切需要一种将细胞刺入速度与细胞机械损伤相关联的模型，并基于模型规划细胞刺入速度以提高细胞后续发育率。
发明内容
[0005]本发明的目的是针对现有技术的不足，提供一种基于应力失效模型的细胞刺入速度规划方法，以揭示细胞刺入速度和细胞机械损伤之间的关联，提高细胞后续发育率。
[0006]为了实现上述目的，本发明通过以下技术方案实现：
[0007]基于应力失效模型的细胞刺入速度规划方法，包括以下步骤：
[0008]步骤1，建立细胞刺入仿真环境；
[0009]步骤2，建立刺入过程中细胞机械损伤的应力失效模型；
[0010]步骤3，检测细胞刺入过程中的应力失效比例；
[0011]步骤4，基于应力失效模型评估不同速度的细胞刺入过程，对刺入速度进行规划。
[0012]进一步，所述步骤1具体为：
[0013]基于物质点法，将细胞建模为由粘弹性材料点组成的外半径为rpe的球壳包裹着由粘塑性和粘性流体材料点混合组成的半径为rcy的球体，球壳和内球分别表示透明带与细胞质；
[0014]吸持针和注射针设置为刚体，形状为空心圆柱体结构；
[0015]将吸持针、注射针接触的细胞材料点的相对速度在指向吸持针、注射针表面的法向方向上被投影为0；
[0016]定义作为细胞质材料点受到的摩擦力为：
[0017]Ffriction＝α·vrelative·m
[0018]其中，m为材料点质量，α为摩擦力的比例系数，vrelative为吸持针、注射针与针管周围的细胞质粒子的相对速度；
[0019]进一步，在吸持针、注射针的针管内和管口附近区域设置吹吸力场：
[0020]针管口处的吹吸力场呈球面圆锥状，球面圆锥的顶点位于针管口内的轴向中心线上，球面圆锥的锥面与针管口内壁相交，方向为其中d为球面圆锥内任一点到球面圆锥顶点的向量；
[0021]在针管内部区域，吹吸力场方向设置为平行于管壁的方向向量；
[0022]定义吹吸力场Fpressure为：
[0023]Fpressure＝β·npressure·m
[0024]其中，m为材料点质量，β为设定的吹吸力大小，npressure为针管内外的吹吸力场的方向。
[0025]进一步，步骤2具体为：
[0026]定义应力失效准则为，在一次细胞刺入操作过程中，一旦细胞材料点p在任意时刻t承受的应力P(Fp)超过阈值Pfail，该粒子失效且不再恢复：
[0027]if:P(Fp)＞Pfail,failp＝1
[0028]其中P(Fp)为粒子的第一Piola-Kirchhoff应力张量，使用共旋弹性模型计算：
[0029]
[0030]
[0031]其中，Fp为粒子的形变梯度，Ψ(Fp)为粒子的弹性势能密度，σpi是Fp的奇异值，μ和λ为粒子的拉梅常数，Jp是粒子体积比，Jp＝det(Fp)，det表示计算矩阵的行列式，为避免数值不稳定，Jp使用下式迭代更新：其中，为粒子的局部速度梯度，tr表示计算矩阵的迹；Rp为计算旋转矩阵，由Rp＝Up*VpT计算得到，其中Up和Vp是由每个粒子的形变梯度Fp进行SVD奇异值分解后的Up,sigp,Vp中得到。
[0032]进一步，应力失效比例计算公式为：
[0033]
[0034]其中，sum(pfail)为应力失效材料点的总数，sum(p)为所有材料点的总数，使用两者的比值反映细胞中损伤区域的相对大小。
[0035]进一步，所述步骤4具体包括：
[0036]4.1计算单次的细胞应力失效比例；
[0037]初始化细胞、吸持针和注射针；设置细胞质半径和透明带半径，两根微针各自的内外径；将细胞置于吸持针与注射针之间；
[0038]吸持针设置吸力场，吸引并固定细胞位置，使得吸持针中心、细胞中心以及注射针中心位于同一条水平线上；在吸持稳定后，注射针以一个固定大小的速度沿水平线朝向细胞中心的方向运动，停留至稳定后，以相同的速度退回初始点，完成细胞刺入仿真；
[0039]每一个时间步计算当前每个材料点所受应力，判断材料点是否失效，统计所有材料点的应力失效比例；
[0040]4.2计算不同速度刺入下的应力失效比例；
[0041]仿真环境初始设置相同，在注射针接触细胞后的刺入深度相同的情况下，改变注射针刺入速度，输出刺入完成后的应力失效比例随着刺入速度的变化情况；
[0042]4.3基于不同速度刺入下的应力失效比例对刺入速度进行规划；
[0043]由于细胞损伤和细胞活性呈负相关性，基于应力失效模型使用不同速度刺入下的应力失效比例评估细胞活性；选择细胞活性最高时相对应的速度进行刺入，从而最小化细胞刺入的机械损伤，提高细胞的后续发育率。
[0044]相对于现有技术，本发明具有的有益效果为：
[0045]本发明实现了一种基于应力失效模型的细胞刺入速度规划方法。该方法较传统方法，能够对大范围的动态刺入速度下的细胞机械损伤结果进行评估，从而规划出最优的刺入速度，提升细胞后续发育率，通过和现实实验中结果进行对比验证，证明本细胞刺入速度规划方法的有效性。
附图说明
[0046]为了更清楚地说明本发明的实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是示例性的，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图引申获得其它的实施附图。
[0047]图1为本发明实施例的方法流程图；
[0048]图2为本发明实施例的细胞刺入环境示意图；
[0049]图3为本发明实施例的吹吸力场示意图；
[0050]图4为本发明实施例的应力失效材料点示意图；
[0051]图5为本发明实施例的单次应力失效比例曲线图；
[0052]图6为本发明实施例的不同速度刺入下应力失效比例曲线图；
[0053]图7为本发明实施例的不同速度刺入下细胞卵裂率曲线图。
具体实施方式
[0054]下面将参照附图更详细地描述本发明的示例性实施例。虽然附图中显示了本发明的示例性实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本发明。
[0055]本实施例为利用该方法对细胞刺入过程中的机械损伤进行评估和刺入速度规划。
[0056]本实施例在Python中使用太极编程语言进行细胞刺入仿真与数据处理。本实施例所使用的细胞是从当地屠宰场所取的家猪卵母细胞。卵母细胞的获取方法如下：
[0057]猪卵巢从屠宰场取出后，在两个小时之内用装有35°到37°的生理盐水的保温瓶运送到实验室。然后立刻用37°的含有100IU/L的盘尼西林和50mg/L的链霉素的无菌生理盐水清洗两次。卵母细胞从卵巢上的2到6毫米直径的小囊中抽取，将抽出的细胞用TL-Hepes-PVA冲洗三次后在39°、二氧化碳浓度5％的培养箱内体外成熟培养(IVM)42小时。经过IVM后，将细胞用0.1％的透明质酸酶进行脱卵。最后用M199将细胞清洗三次，得到的这些细胞为本例中所用卵母细胞。
[0058]下面结合附图对本发明做进一步说明，由图1所示，本发明提供的实现细胞刺入速度规划的具体实施情况包括以下步骤：
[0059]1.建立了细胞刺入仿真环境。
[0060]1.1建立复合性质的三维细胞模型。
[0061]由于卵细胞结构的复杂性，细胞的弹性性质主要由外部的透明带体现，而内部细胞质则更类似于流体。因此，本发明基于物质点法，将细胞建模为由粘弹性材料点组成的外半径为rpe的球壳，包裹着由粘塑性和粘性流体材料点混合组成的半径为rcy的球体，从而分别表示透明带与细胞质的性质。为平衡计算量和精度，本实施例将仿真中边长为1的立方体空间划分为96*96*96的欧拉背景网格，比例尺为1:500μm，按实际细胞比例设置细胞质半径为50μm，透明带球壳外半径为75μm，使用65536个材料点随机初始化组成细胞。
[0062]为保证透明带与细胞质材料点的随机均匀分布，采用以下算法初始化三维球壳区域内材料点的位置(对于球体，可视为内半径为0的球壳)：
[0063]设球壳内采样点与球体中心的距离为r，大小在球壳内半径r0和外半径r1之间，即r0＜r＜r1。为保证三维采样的均匀性，则采样点概率的累积分布函数需满足在[0,1]之间的均匀分布函数。因此，倒推回随机生成采样点的距离R可用下列式子计算，以保证均匀采样：
[0064]
[0065]其中，u为在[0,1]之间均匀采样的值。
[0066]则随机生成的一个采样点相对球体中心的三维坐标可表示为：
[0067]
[0068]其中，θ为在[0,2π]之间均匀采样的值，v为在[-1,1]之间均匀采样的值。
[0069]根据上述算法构建细胞模型，对于透明带的球壳区域，外半径为透明带半径rpe，内半径为细胞质半径rcy；对于细胞质的球体区域，外半径为细胞质半径rcy，内半径为0。分别多次均匀采样材料点，最终生成透明带与细胞质，建立三维细胞模型。
[0070]1.2设置吸持针与注射针。
[0071]两根微针设置为刚体，形状为空心圆柱体结构，本实施例按照实际比例关系，将吸持针的内径设为50μm，外径设为160μm，注射针的外径设为20μm。考虑到微针质量相对于细胞可视为无穷大，细胞与微针发生接触时由于接触力的作用被挤压，因此将微针接触的细胞材料点的相对速度在指向微针表面的法向方向上被投影为0。考虑到在相对速度较小时，粘稠流体的阻力的大小与相对速度成正比，定义细胞质材料点受到的摩擦力为：
[0072]Ffriction＝α·vrelative·m
[0073]其中，m为材料点质量，α为摩擦力的比例系数，本实施例设置为0.05，vrelative为微针与针管周围细胞质粒子的相对速度。
[0074]为了体现微针管内气压对操作的影响，在针管内和管口附近区域设置吹吸力场。如图3所示，针管口处的吹吸力场呈球面圆锥状，球面圆锥的顶点位于针管口内的轴向中心线上，球面圆锥顶点到针管口平面的垂直距离为a，本实施例中a取100μm，球面圆锥顶点到球面的半径为b，本实施例中b取150μm，球形圆锥的锥面与针管口内壁相交，方向为其中d为球面圆锥内任一点到球面圆锥顶点的向量；在针管内部区域，吹吸力场方向设置为平行于管壁的方向向量。定义吹吸力场Fpressure为：
[0075]Fpressure＝β·npressure·m
[0076]其中，m为材料点质量，β为设定的吹吸力大小，本实施例设置为2500以保证吸持成功，npressure为针管内外的吹吸力场的方向。
[0077]建模完成的细胞刺入仿真环境如图2(b)所示，图2(a)为实际细胞刺入环境图像。
[0078]2.建立刺入过程中细胞损伤的应力失效模型。
[0079]仿真环境中可同时计算细胞材料点在每一个时间步上所受的应力，考虑到细胞内部受力与细胞机械损伤息息相关，基于应力数据建立刺入过程中的应力失效模型。定义应力失效准则为，在一次细胞刺入操作过程中，一旦细胞材料点p在任意时刻t承受的应力P(Fp)超过阈值Pfail，该粒子失效且不再恢复，本实施例设置Pfail为2000，即：
[0080]if:P(Fp)＞Pfail,failp＝1
[0081]其中P(Fp)为粒子的第一Piola-Kirchhoff应力张量，可使用共旋(corotated)弹性模型计算：
[0082]
[0083]
[0084]其中，Fp为粒子的形变梯度，Ψ(Fp)为粒子的弹性势能密度，σpi是Fp的奇异值。μ和λ为粒子的拉梅(lame)常数，。
[0085]Jp是粒子体积比，Jp＝det(Fp)，det表示计算矩阵的行列式，为避免数值不稳定，Jp使用下列式子迭代更新：其中，为粒子的局部速度梯度，tr表示计算矩阵的迹。
[0086]Rp为计算旋转矩阵，由Rp＝Up*VpT计算得到，其中Up和Vp是由每个粒子的形变梯度Fp进行SVD奇异值分解后的Up,sigp,Vp中得到。
[0087]从而，可以判断在操作过程中的特定材料点p是否受到过大应力而被损伤。为进一步宏观评估操作中细胞机械损伤量，定义应力失效比例为：
[0088]
[0089]其中，sum(pfail)为应力失效材料点的总数，sum(p)为所有材料点的总数，使用两者的比值反映细胞中损伤区域的相对大小。
[0090]图4为建立刺入过程中细胞机械损伤的应力失效模型，图中绿点表示失效的细胞材料点。
[0091]3.对细胞刺入过程的应力失效比例进行检测。
[0092]单次刺入的细胞应力失效比例曲线图如图5所示。
[0093]放置使吸持针中心、细胞中心以及注射针中心位置位于同一条水平线上，等待200个时间步至吸持稳定，开启细胞刺入操作。注射针以一个固定大小的速度沿水平线朝向细胞中心的方向运动150μm，并停留50个时间步，之后以相同的速度退回初始点，完成细胞刺入仿真。本实施例使用1μm/时间步的速度刺入细胞并退回，图5记录了细胞应力失效比例随操作的变化情况，在细胞与注射针接触的过程中，细胞应力失效比例快速上升，反应了动态刺入过程对细胞造成的机械损伤。
[0094]4.基于不同速度刺入下的应力失效比例对刺入速度进行规划，实现了基于应力失效模型的细胞刺入速度规划方法。
[0095]图6为不同速度刺入下应力失效比例曲线图。
[0096]仿真环境初始设置相同，在注射针接触细胞后的刺入深度相同的情况下，改变注射针刺入速度，输出刺入完成后的应力失效比例随着刺入速度的变化情况。本实施例测试了13组不同的刺入速度对应的应力失效比例结果，每组重复20次后取平均值。图6记录了不同速度刺入的应力失效比例曲线图，随着刺入速度的增加，细胞的应力失效比例先增加再减少再单调增加。考虑到过小的刺入速度会导致更长的操作时间，实验中应尽量减少操作时间以降低细胞处于不适宜的操作外环境(温度、湿度、气体浓度等)带来的不利影响，因此过高或过低的刺入速度都不利于细胞后续发育，最优刺入速度为图中导数为0时的极小值点对应的刺入速度。
[0097]图7为实验所得的不同速度刺入下细胞卵裂率曲线图。
[0098]为了验证本发明方法的有效性，使用实际的细胞刺入实验的得到的细胞发育结果进行比较验证。本实施例共使用192枚卵母细胞，对9组不同的速度的进行刺入实验，每组约包含20个细胞。刺入实验使用20μm的标准注射针完成，刺入时的细胞姿态随机，刺入完成后进行孤雌激活，继续培养两天后统计细胞卵裂率，最终结果如图7所示。将图7(a)不同速度刺入实验得到的卵裂率变化与图7(b)红框中本发明方法预测的细胞活性变化结果相对比，对应关系良好，细胞卵裂率随着刺入速度的增加呈现了先增加后减小的趋势。证明本发明方法有效且对实际的实验操作有着一定的指导意义。
[0099]以上通过实施例对本发明进行了详细说明，但所述内容仅为本发明的示例性实施例，不能被认为用于限定本发明的实施范围。本发明的保护范围由权利要求书限定。凡利用本发明所述的技术方案，或本领域的技术人员在本发明技术方案的启发下，在本发明的实质和保护范围内，设计出类似的技术方案而达到上述技术效果的，或者对申请范围所作的均等变化与改进等，均应仍归属于本发明的专利涵盖保护范围之内。应当注意，为了清楚的进行表述，本发明的说明中省略了部分与本发明的保护范围无直接明显的关联但本领域技术人员已知的处理表述。