为了研究直接心室辅助的生物力学影响以及探究最优的加载模式，本文基于有限元方法建立了心衰患者的左心室模型，并提出了一种维持压迫力峰值的加载模式，从血流动力学和生物力学两个方面与传统的正弦加载模式进行了对比。结果表明，两种模式都能显著提升血流动力学参数，射血分数分别从基线29.33%增加到37.32%与37.77%，峰值压力、每搏量和每搏功等参数都有所增加；且两种模式的应力集中、过度纤维应变等现象均有所改善。然而，当考虑到辅助装置工作周期的相位误差时，本文所提出的辅助模式受到的影响更小，故本文研究或可为直接心室辅助装置的设计和优化提供理论支持。
引用本文： 李宸, 强贤杰, 张盛, 王天博, 刘晓翰, 张悦, 黄刚, 张小刚, 徐俊波, 靳忠民. 基于有限元方法的直接心室辅助加载模式研究. 生物医学工程学杂志, 2024, 41(4): 782-789. doi: 10.7507/1001-5515.202312070 复制
0 引言
心力衰竭（hearth failure，HF）（以下简称：心衰），是心脏结构或功能异常导致心室充盈或射血能力受损的一组复杂临床综合征。2012—2015年进行的中国高血压调查研究显示，我国35岁以上成人心衰患病率约为1.3%，涉及1 370万人[1-2]。目前，基于循证医学证据的主要心衰治疗措施包括“新四联”药物治疗和心脏再同步化治疗等器械治疗[3]。通过药物及现有的器械治疗措施，虽可改善近期血流动力学指标并减轻症状，但有相当比例终末期心衰患者对上述措施低反应[2]，因而需进一步接受心室辅助装置植入[4]。现有的主流心室辅助装置，通过机械泵连接左心室（left ventricle，LV）和主动脉血流来实现心脏泵功能和血液输送及回流[5]。由于其直接接触血液，必须长期服用抗凝血药物，不可避免出现败血症、血栓、心衰加重等问题[6]，且手术创伤大，不利于康复。
直接心室辅助装置（direct cardiac compression，DCC）可包裹住心脏并与心脏的跳动进行同频柔性挤压来辅助泵血，无需在心室及主动脉开孔，大大降低了手术难度及术后创伤[7-8]。虽然DCC已有40年的研究历史，但由于心脏按压的生物力学效应相关研究很少，其临床应用一直受到限制[9]，如何以正确的方式将压迫施加给心脏是一个值得讨论的问题。
目前，已有的仿真模型大多以理想的正弦曲线力加载在心室壁上，且认为压迫峰值与收缩末期重合[10-11]。然而在现实中，由于DCC结构的复杂性，要保持与心脏自然射血周期同频工作是一件很困难的事情，施加压迫力开始时刻与心脏开始收缩的时刻有不可避免的误差。这种误差可能会影响装置的辅助效果，增加患者的负荷。为此，参考外科医生施加压迫的手法[12]，本文提出了在收缩初期开始施加压迫，射血开始时达到压迫力峰值，并在整个射血期间维持压迫力峰值的压迫方式。由于该方式维持了一段时间的峰值压迫力，即使在DCC启动时刻存在滞后的情况下，压迫力峰值与心肌收缩力峰值也会重合，从而减小DCC辅助效果受到上述误差的影响。
为了验证该方式的合理性，本研究通过有限元方法构建人体LV的仿真模型，对比分析正弦加载和在射血阶段维持峰值力的两种加载模式的辅助效果，同时通过改变压迫开始时间来模拟相位误差，进一步对比两种工作模式受误差影响的结果。通过上述研究，本文期望能探究更优的DCC工作模式，以实现在尽量减小心肌的应力和应变的同时最大化提高输出量。
1 材料与方法
本研究使用的心脏核磁共振成像数据来自心脏图谱项目（Cardiac Atlas Project）中的桑尼布鲁克心脏数据集（Sunnybrook cardiac data，SCD），该数据集对公众提供开放访问和数据下载[13]，其研究对象为一名42岁的心衰男性。
本文所描述的三维人体LV模型由以下部分组成：患者个性化LV网格模型、肌纤维结构、被动心肌本构、主动张力模型、体循环参数。求解器采用有限元分析软件Abaqus 2018（SIMULIA Inc., 美国），本构模型通过有限元分析软件Abaqus 2018（SIMULIA Inc., 美国）自带的用户子程序VUMAT实现。
1.1 左心室几何模型及纤维结构的建立
LV几何模型选取舒张初期为初始构型，此时LV的体积和压力都最小且最接近空载构型状态。心衰LV的几何模型基于核磁共振图像由Guan等[14]的公开算法生成，如图1所示，共72 956个四面体单元，16 268个节点。心肌纤维方向定义为从心外膜到心内膜呈线性变化。根据Streeter等[15]先前的实验测量，LV的心内膜肌纤维螺旋角（αendo）、心外膜肌纤维螺旋角（αepi）规定为：αendo = 60°、αepi = − 60°，结果如图1所示。
图1
LV有限元模型
Figure1.
LV finite element model
已有研究通过实验证明，心肌纤维不会沿着平均纤维方向完全排列，而是离散分布的[16-17]。本文遵循纤维离散模型[18]，参考Guan等[14]的参数，心肌纤维在单位半球区域内的平均纤维方向呈正态分布。
1.2 心肌材料本构与参数
本文采用了Holzapfel等[19]提出的基于不变量的纤维增强本构模型，该模型考虑了心肌的形态和结构，将LV心肌的被动特性建模为非均匀的、厚壁的、不可压缩的、正交各向异性的非线性弹性材料，具体参数如表1所示。心肌的主动应力如式(1)所示：
表1
被动本构模型参数
Table1.
Passive constitutive model parameters
 |
其中，Ta为沿纤维方向产生的主动张力，
，F为变形梯度，f0为平均心肌纤维方向，s0为平均片层纤维方向，n0为s0的法向向量，
，
，
，nf、ns、nn为主动张力在各自方向上的比例，参考Guan等[20]的参数，nf = 0.086，ns = 0.268，nn = 0.646，“⊗”代表张量积。
主动张力的时变弹性模型[21]，如式(2)所示：
 |
其中，Tmax是最大收缩力，代表心肌组织能够产生的最大等距张力。Ca0为细胞内钙离子浓度的峰值，ECa50为钙敏感性，Ct为时间相关性函数，如式(3)～式(5)所示：
 |
 |
 |
其中，(Ca0)max是细胞内钙浓度的最大峰值，B是一个常数，l为当前肌节长度，
。l0是主动张力时的肌节长度，t0是张力峰值的时间，tr是松弛的时间，
。m、b为支配松弛时间与肌节长度关系的常数，lr为无应力肌节长度。E是格林—拉格朗日应变张量，其中下角标f、s和n分别表示肌纤维、片层和片层法线方向。
心衰LV的被动心肌参数分别由Mangion等[22]的模型参数重新缩放，根据舒张末期容积（end-diastolic volume，EDV）确定，如表1所示。其中a、b、af、bf、as、bs、afs和bfs是与纤维方向有关的8个材料常数。另一方面，根据收缩末期容积（end-systolic volume，ESV）缩放主动收缩模型参数的Tmax，所得数值如表2所示。
表2
主动收缩模型参数
Table2.
Active shrink model parameters
1.3 循环系统模型
根据血流动力学与电学系统的等效规律[23]，将简化闭环集总参数循环系统与LV模型耦合，包括主动脉腔（aorta，Ao）和左心房（left atrium，LA）腔。如图1所示循环系统模型中，二尖瓣（mitral valve，MV）和主动脉瓣（aortic valve，AV）的功能通过单向流体交换来实现。CVAV、CVMV、CVSys分别代表AV、MV和体循环的阻抗，CVAV = 1.2、CVMV = 1.0、CVSys = 100，单位为MPa·mm²/t。kAo、kLA为两个接地的弹簧，代表LV和Ao的顺应性，kAo = 0.8、kLA = 1.0，单位为k/(N·mm)。通过有限元分析软件Abaqus 2018（SIMULIA Inc.，美国）的流体腔和流体交换模块实现了简化的循环模型。该模型类似于弹性腔模型（windkessel models）[24]。
参考正常生理参数的研究，设置每个腔室的加载压力作为本研究的初始条件，由于心衰患者的舒张末期压力相比健康受试者通常高出10～25 mm Hg[25-26]，设置初始条件为LV腔：14 mm Hg、Ao：100 mm Hg、LA腔：14 mm Hg。整个过程如下：首先将模型LV预加载至正常LV舒张末期压力，以达到舒张末期状态，然后开始等容收缩，心室的压力增加到超过Ao的压力，AV打开，LV开始射血（t ≈ 0.1 s），整个心动周期为0.65 s。
1.4 不同模式的心外膜压迫
Aranda-Michel等[27]通过在理想椭球体LV上施加压迫得出结论：对LV压迫区域进行细分，对于血流动力学几乎没有影响。因此，本文选择覆盖LV中部大部分区域的压迫范围，覆盖面积为11 278.63 mm²，如图2所示。
图2
压迫模式示意图
Figure2.
Compression mode diagram
传统直接辅助压迫的正弦曲线，开始于LV射血初始时刻，在射血结束时达到峰值，在等容舒张末期附近降为零[27]，本文将这种压迫模式记作模式一；本文提出的另一种压迫幅值曲线在心室收缩时压迫力逐渐增加，直到开始射血时达到峰值，并在整个射血期间，受力均保持在峰值不变，之后逐渐降低为零，本文将该曲线记作模式二；两种压迫模式的曲线如图2所示，计算公式如式(6)～式(7)所示：
 |
 |
其中，p(t)为压迫力幅值；t=0为舒张末期时刻，心室开始收缩；0～0.1 s为等容收缩阶段；0.1～0.2 s为快速射血期，之后心肌收缩力减小；0.2～0.3 s为减慢射血期，之后AV关闭，射血结束，0.65 s为一个心跳周期的总时长。
为了模拟DCC工作周期的误差，参考第一心音时长约为100 ms，将模式一和模式二从最初的开始时间分别向后移动了10 ms、20 ms和50 ms（第一心音时长的50%），所有的模拟均运行4次循环，LV压力与容积变化趋于稳定。
Obiajulu等[26]通过使用半球形DCC装置的简单模型来估计DCC装置致动器的力需求，参考其施加的压力范围，施加峰值为5 kPa的负载（本模型LV最大血压约为16 kPa）。
1.5 模型有效性验证
如表3所示，本文对比了EDV、ESV和射血分数（ejection fractions，EF）的模拟值和核磁共振图像的测量值，可以看出EDV、ESV和EF的模拟值和测量值误差都在1%以内。
表3
模型结果对比
Table3.
Comparison of model results
另一方面，参考Kuijer等[28]根据短轴和长轴核磁共振图像中标记点的位移，描述了整个人类LV三维位移和应变演变的正常范围。将本文LV三个部位的应变结果与之进行对比，包括周向应变和轴向应变。这些应变计算均以舒张末期作为初始构型，选择单个心动周期内收缩末期作为比较时间节点，所有应变均表示为格林—拉格朗日应变，如式（8）所示：
 |
 |
其中，ec、el分别是周向和纵向方向上的单位向量，由前文计算可得，C为右柯西—格林变形张量。因此，格林—拉格朗日应变定义如式(10)所示：
 |
对比结果如表4所示，表4中所有数据均表示为均值±标准差，可以看出，除了在心尖部的轴向应变外，其它所有应变都在临床数据的两个标准偏差范围内，以上结果共同证明了模型的有效性。
表4
收缩末期应变对比（x ± s）
Table4.
Comparison of end-systolic strain（x ± s）
2 结果
2.1 两种辅助模式效果对比
评价血流动力学的最佳方法是使用压力—容积曲线，两种压迫模式与心衰LV的压力—容积曲线如图3所示。
图3
压力容积曲线
Figure3.
Pressure-volume loop
整体上可以看出，两种模式压迫下的压力—容积曲线几乎重叠，说明两种模式的每搏输出量（stroke volume，SV）、EF、峰值压力（peak pressure，PP）和每搏功（stroke work，SW）对于未施加压迫的LV都有所提升，且增加量近似相等，两种压迫模式下的LV血流动力学参数具体值如表5所示。
表5
血流动力学参数
Table5.
Hemodynamic parameters
在收缩末期时心肌收缩力达到最大，此时心肌纤维主应力（以σff表示）分布如图4所示。施加压迫之后，位于心尖部的应力集中有明显改善，且整体的应力分布有所减小。
图4
收缩末期LV心肌纤维应力分布
Figure4.
Fiber stress distribution in the LV end-systole
2.2 考虑相位误差情况下两种辅助模式效果的对比
延迟两种模式压迫开始时间，血流动力学参数结果如图5所示，包括EF、PP和SV。从中可以看出， 随着压迫开始延迟时间的增加，模式一的各项血流动力学参数都有所下降，当延迟时间增大到50 ms时，EF下降了2%，PP下降了6 mm Hg，SV下降了4 mL，而模式二的各项血流动力学参数都基本保持稳定。
图5
血流动力学参数
Figure5.
Hemodynamic parameters
心尖扭转角是LV尖部相对于基部旋转的角度，在射血过程中，心内膜纤维间质变形产生势能，随后扭转变形[29]。模式一和模式二在0.00～0.20 s之间收缩过程中的心尖扭转角变化如图6所示。可以看出，当压迫开始时刻与射血开始时刻同步时，模式一在收缩末期的心尖扭转角为18.6°，模式二在收缩末期的心尖扭转角为20.3°，相差不大。随着压迫开始时间的滞后，模式一的心尖扭转角逐渐减小，在50 ms时为13.7°；模式二的心尖扭转角度依旧保持稳定。
图6
心尖扭转角
Figure6.
Heart apex twist angle
2.3 两种辅助模式对左心室生物力学的影响
不同开始时间点的LV收缩末期相对于心跳周期开始时刻的纤维应变分布如图7所示，图7中数据均为真实应变（以LE11表示）。从图7中可以看出，两种模式的压迫均可以减小LV收缩末期的最大纤维应变。未施加外膜辅助压迫时，最大纤维应变为14.6%。当压迫开始时间没有延迟时，模式一和模式二的最大纤维应变分别为10.0%和9.5%，近似相等。随着压迫开始时间的滞后，两种模式的心肌纤维应变均增大，但模式一增加的更多，当滞后时间为50 ms时，模式一的最大纤维应变为12.9%，模式二的最大纤维应变为10.2%。相对于模式一，改变压迫开始时间对模式二的最大心肌纤维应变几乎没有影响。
图7
不同压迫开始时间的纤维应变分布
Figure7.
Fiber strain distribution at different compression start time
3 讨论
本文针对心衰患者构建了个性化LV有限元模型，在此基础上分别施加了正弦曲线及维持一段时间PP的心外膜辅助模式。
对比来看，在血流动力学参数方面，两种辅助模式下的EF、PP和SV都有所上升，其中LV的EF分别从29.33 %提升到37.32 % 和37.77 %，可以维持患者的正常生理活动（EF降低型心衰的基线EF为40%）。
从生物力学的结果来看，两种模式的压迫都可以增加LV收缩能力及改善应力集中。收缩末期时LV处于收缩状态，心肌纤维长度小于静息长度，如图7所示，最大心肌纤维应变减小说明LV的收缩能力增强。从图4可以看出，由于心尖部心肌壁较薄，所以有明显的应力集中现象，长时间的应力集中可能导致局部心肌梗死，施加心外膜辅助可以很大程度改善应力集中。
本文认为有两个方面的原因导致了应力和应变的减小。第一，由于心肌纤维本身的结构特点，心肌纤维在心外膜处与周向成− 60°，在心内膜处为60°，心内壁区域下心肌纤维更接近于轴向方向，并且心肌纤维在沿着纤维方向上的硬度远远高于横向方向。在心外膜施加压迫时，心室被拉长，纤维有向心尖处伸展的趋势，处于收缩状态的心肌沿纤维方向的加载减小，因此心肌纤维的应力和应变都减小。第二，施加心外膜辅助可增强LV的扭转能力。如图6所示，施加模式一和模式二的压迫后，收缩末期心尖扭转角从8.0°增加到了18.6°和20.3°，恢复正常扭转能力，该数据参考Badano等[29]测量的健康人收缩末期心尖扭转角约为16°。当纤维达到收缩峰值时，沿纤维方向的应力会很大，且心内膜收缩力大于心外膜，但由于心外膜半径较大，心外膜处可以产生更大的扭矩以控制旋转方向。在射血期间，心外膜下心肌较大的扭矩经壁耦合至心肌中层和心内膜下心肌，导致LV整体旋转。在心外膜下，这种扭矩有助于主纤维方向的收缩；在中层，扭转动作增强了圆周方向的缩短；在心内膜下，扭转动作可引起肌纤维重新排列，使得心内膜下肌纤维沿纤维的法向方向更紧密地排列以增厚左室壁。扭转减小会增加心内膜应力和应变，从而增加需氧量，降低LV收缩功能的效率[29]，提高心室的扭转角有助于使纤维应力和肌纤维的缩短沿室壁均匀分布。
进一步推迟压迫开始时间以模拟相位误差的影响，根据血流动力学结果显示，压迫开始时间的滞后会削弱模式一辅助效果，EF、SV和PP都有所下降。具体而言，EF在压迫延迟50 ms时从37.32%下降到35.60%，相对于理想情况下的EF提升，效果降低。从图5可以看出，模式二的各项血流动力学结果基本保持稳定不受影响。这是由于模式二在整个射血期间都维持了压迫力峰值，尽管压迫开始时间滞后，在收缩末期时心外膜压迫力仍与心肌纤维收缩力同时处于峰值，因此对误差的抵抗能力较强。从生物力学角度来看，模式一的心尖扭转能力下降，在压迫延迟50 ms时从18.6°下降到13.7°，模式二基本维持在20.0°左右。模式一的心尖扭转能力减弱，心室收缩效率也降低，心肌纤维需要更大的变形来实现收缩动作，如图7所示，模式一的心肌的纤维最大应变从10.0%增大到12.9%，模式二的心肌纤维最大应变依旧保持不变，说明本文提出的加载模式在减小心肌应变方面能力更强，从而在引导心衰心肌的反向良性重塑方面更有优势。
本文的研究也存在一定的局限性：首先，本文仅考虑了LV的模型，而未对其他心腔的运动情况进行模拟。其次，所建立的简化体外循环系统耦合并不能完全模拟人体的血液循环系统。在下一步工作中，将考虑多个心腔的交互，搭建更为复杂的系统，同时根据不同程度的心衰情况建立不同类型的心脏模型，以进一步比较不同的心外膜辅助方式对心脏生物力学的影响，旨在为DCC的设计与优化、临床治疗和术后康复提供理论支持。
4 结论
本文所提出的维持一段时间峰值压迫力的模式和传统正弦加载压迫力模式相比，带来的血流动力学参数提升近似相等，并且心肌的生物力学结果也近似相同。然而，由于DCC工作周期和心脏跳动周期存在的相位误差，正弦加载压迫力模式的血流动力学参数有所下降，心室的收缩能力也有所降低。相比之下，本文所提出的维持峰值压迫力的模式，可以在误差为50 ms的情况下保持血流动力学参数及心肌纤维应力、应变不受影响，因此可对DCC设计和优化提供理论指导。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：李宸、徐俊波和张小刚构思方案，李宸实施方案并处理数据，张盛负责作图，强贤杰、王天博、刘晓翰、张悦、黄刚、徐俊波和靳忠民讨论并修改了论文。
0 引言
心力衰竭（hearth failure，HF）（以下简称：心衰），是心脏结构或功能异常导致心室充盈或射血能力受损的一组复杂临床综合征。2012—2015年进行的中国高血压调查研究显示，我国35岁以上成人心衰患病率约为1.3%，涉及1 370万人[1-2]。目前，基于循证医学证据的主要心衰治疗措施包括“新四联”药物治疗和心脏再同步化治疗等器械治疗[3]。通过药物及现有的器械治疗措施，虽可改善近期血流动力学指标并减轻症状，但有相当比例终末期心衰患者对上述措施低反应[2]，因而需进一步接受心室辅助装置植入[4]。现有的主流心室辅助装置，通过机械泵连接左心室（left ventricle，LV）和主动脉血流来实现心脏泵功能和血液输送及回流[5]。由于其直接接触血液，必须长期服用抗凝血药物，不可避免出现败血症、血栓、心衰加重等问题[6]，且手术创伤大，不利于康复。
直接心室辅助装置（direct cardiac compression，DCC）可包裹住心脏并与心脏的跳动进行同频柔性挤压来辅助泵血，无需在心室及主动脉开孔，大大降低了手术难度及术后创伤[7-8]。虽然DCC已有40年的研究历史，但由于心脏按压的生物力学效应相关研究很少，其临床应用一直受到限制[9]，如何以正确的方式将压迫施加给心脏是一个值得讨论的问题。
目前，已有的仿真模型大多以理想的正弦曲线力加载在心室壁上，且认为压迫峰值与收缩末期重合[10-11]。然而在现实中，由于DCC结构的复杂性，要保持与心脏自然射血周期同频工作是一件很困难的事情，施加压迫力开始时刻与心脏开始收缩的时刻有不可避免的误差。这种误差可能会影响装置的辅助效果，增加患者的负荷。为此，参考外科医生施加压迫的手法[12]，本文提出了在收缩初期开始施加压迫，射血开始时达到压迫力峰值，并在整个射血期间维持压迫力峰值的压迫方式。由于该方式维持了一段时间的峰值压迫力，即使在DCC启动时刻存在滞后的情况下，压迫力峰值与心肌收缩力峰值也会重合，从而减小DCC辅助效果受到上述误差的影响。
为了验证该方式的合理性，本研究通过有限元方法构建人体LV的仿真模型，对比分析正弦加载和在射血阶段维持峰值力的两种加载模式的辅助效果，同时通过改变压迫开始时间来模拟相位误差，进一步对比两种工作模式受误差影响的结果。通过上述研究，本文期望能探究更优的DCC工作模式，以实现在尽量减小心肌的应力和应变的同时最大化提高输出量。
1 材料与方法
本研究使用的心脏核磁共振成像数据来自心脏图谱项目（Cardiac Atlas Project）中的桑尼布鲁克心脏数据集（Sunnybrook cardiac data，SCD），该数据集对公众提供开放访问和数据下载[13]，其研究对象为一名42岁的心衰男性。
本文所描述的三维人体LV模型由以下部分组成：患者个性化LV网格模型、肌纤维结构、被动心肌本构、主动张力模型、体循环参数。求解器采用有限元分析软件Abaqus 2018（SIMULIA Inc., 美国），本构模型通过有限元分析软件Abaqus 2018（SIMULIA Inc., 美国）自带的用户子程序VUMAT实现。
1.1 左心室几何模型及纤维结构的建立
LV几何模型选取舒张初期为初始构型，此时LV的体积和压力都最小且最接近空载构型状态。心衰LV的几何模型基于核磁共振图像由Guan等[14]的公开算法生成，如图1所示，共72 956个四面体单元，16 268个节点。心肌纤维方向定义为从心外膜到心内膜呈线性变化。根据Streeter等[15]先前的实验测量，LV的心内膜肌纤维螺旋角（αendo）、心外膜肌纤维螺旋角（αepi）规定为：αendo = 60°、αepi = − 60°，结果如图1所示。
图1
LV有限元模型
Figure1.
LV finite element model
已有研究通过实验证明，心肌纤维不会沿着平均纤维方向完全排列，而是离散分布的[16-17]。本文遵循纤维离散模型[18]，参考Guan等[14]的参数，心肌纤维在单位半球区域内的平均纤维方向呈正态分布。
1.2 心肌材料本构与参数
本文采用了Holzapfel等[19]提出的基于不变量的纤维增强本构模型，该模型考虑了心肌的形态和结构，将LV心肌的被动特性建模为非均匀的、厚壁的、不可压缩的、正交各向异性的非线性弹性材料，具体参数如表1所示。心肌的主动应力如式(1)所示：
表1
被动本构模型参数
Table1.
Passive constitutive model parameters
|
其中，Ta为沿纤维方向产生的主动张力，，F为变形梯度，f0为平均心肌纤维方向，s0为平均片层纤维方向，n0为s0的法向向量，，，，nf、ns、nn为主动张力在各自方向上的比例，参考Guan等[20]的参数，nf = 0.086，ns = 0.268，nn = 0.646，“⊗”代表张量积。
主动张力的时变弹性模型[21]，如式(2)所示：
|
其中，Tmax是最大收缩力，代表心肌组织能够产生的最大等距张力。Ca0为细胞内钙离子浓度的峰值，ECa50为钙敏感性，Ct为时间相关性函数，如式(3)～式(5)所示：
|
|
|
其中，(Ca0)max是细胞内钙浓度的最大峰值，B是一个常数，l为当前肌节长度，。l0是主动张力时的肌节长度，t0是张力峰值的时间，tr是松弛的时间，。m、b为支配松弛时间与肌节长度关系的常数，lr为无应力肌节长度。E是格林—拉格朗日应变张量，其中下角标f、s和n分别表示肌纤维、片层和片层法线方向。
心衰LV的被动心肌参数分别由Mangion等[22]的模型参数重新缩放，根据舒张末期容积（end-diastolic volume，EDV）确定，如表1所示。其中a、b、af、bf、as、bs、afs和bfs是与纤维方向有关的8个材料常数。另一方面，根据收缩末期容积（end-systolic volume，ESV）缩放主动收缩模型参数的Tmax，所得数值如表2所示。
表2
主动收缩模型参数
Table2.
Active shrink model parameters
1.3 循环系统模型
根据血流动力学与电学系统的等效规律[23]，将简化闭环集总参数循环系统与LV模型耦合，包括主动脉腔（aorta，Ao）和左心房（left atrium，LA）腔。如图1所示循环系统模型中，二尖瓣（mitral valve，MV）和主动脉瓣（aortic valve，AV）的功能通过单向流体交换来实现。CVAV、CVMV、CVSys分别代表AV、MV和体循环的阻抗，CVAV = 1.2、CVMV = 1.0、CVSys = 100，单位为MPa·mm²/t。kAo、kLA为两个接地的弹簧，代表LV和Ao的顺应性，kAo = 0.8、kLA = 1.0，单位为k/(N·mm)。通过有限元分析软件Abaqus 2018（SIMULIA Inc.，美国）的流体腔和流体交换模块实现了简化的循环模型。该模型类似于弹性腔模型（windkessel models）[24]。
参考正常生理参数的研究，设置每个腔室的加载压力作为本研究的初始条件，由于心衰患者的舒张末期压力相比健康受试者通常高出10～25 mm Hg[25-26]，设置初始条件为LV腔：14 mm Hg、Ao：100 mm Hg、LA腔：14 mm Hg。整个过程如下：首先将模型LV预加载至正常LV舒张末期压力，以达到舒张末期状态，然后开始等容收缩，心室的压力增加到超过Ao的压力，AV打开，LV开始射血（t ≈ 0.1 s），整个心动周期为0.65 s。
1.4 不同模式的心外膜压迫
Aranda-Michel等[27]通过在理想椭球体LV上施加压迫得出结论：对LV压迫区域进行细分，对于血流动力学几乎没有影响。因此，本文选择覆盖LV中部大部分区域的压迫范围，覆盖面积为11 278.63 mm²，如图2所示。
图2
压迫模式示意图
Figure2.
Compression mode diagram
传统直接辅助压迫的正弦曲线，开始于LV射血初始时刻，在射血结束时达到峰值，在等容舒张末期附近降为零[27]，本文将这种压迫模式记作模式一；本文提出的另一种压迫幅值曲线在心室收缩时压迫力逐渐增加，直到开始射血时达到峰值，并在整个射血期间，受力均保持在峰值不变，之后逐渐降低为零，本文将该曲线记作模式二；两种压迫模式的曲线如图2所示，计算公式如式(6)～式(7)所示：
|
|
其中，p(t)为压迫力幅值；t=0为舒张末期时刻，心室开始收缩；0～0.1 s为等容收缩阶段；0.1～0.2 s为快速射血期，之后心肌收缩力减小；0.2～0.3 s为减慢射血期，之后AV关闭，射血结束，0.65 s为一个心跳周期的总时长。
为了模拟DCC工作周期的误差，参考第一心音时长约为100 ms，将模式一和模式二从最初的开始时间分别向后移动了10 ms、20 ms和50 ms（第一心音时长的50%），所有的模拟均运行4次循环，LV压力与容积变化趋于稳定。
Obiajulu等[26]通过使用半球形DCC装置的简单模型来估计DCC装置致动器的力需求，参考其施加的压力范围，施加峰值为5 kPa的负载（本模型LV最大血压约为16 kPa）。
1.5 模型有效性验证
如表3所示，本文对比了EDV、ESV和射血分数（ejection fractions，EF）的模拟值和核磁共振图像的测量值，可以看出EDV、ESV和EF的模拟值和测量值误差都在1%以内。
表3
模型结果对比
Table3.
Comparison of model results
另一方面，参考Kuijer等[28]根据短轴和长轴核磁共振图像中标记点的位移，描述了整个人类LV三维位移和应变演变的正常范围。将本文LV三个部位的应变结果与之进行对比，包括周向应变和轴向应变。这些应变计算均以舒张末期作为初始构型，选择单个心动周期内收缩末期作为比较时间节点，所有应变均表示为格林—拉格朗日应变，如式（8）所示：
|
|
其中，ec、el分别是周向和纵向方向上的单位向量，由前文计算可得，C为右柯西—格林变形张量。因此，格林—拉格朗日应变定义如式(10)所示：
|
对比结果如表4所示，表4中所有数据均表示为均值±标准差，可以看出，除了在心尖部的轴向应变外，其它所有应变都在临床数据的两个标准偏差范围内，以上结果共同证明了模型的有效性。
表4
收缩末期应变对比（x ± s）
Table4.
Comparison of end-systolic strain（x ± s）
2 结果
2.1 两种辅助模式效果对比
评价血流动力学的最佳方法是使用压力—容积曲线，两种压迫模式与心衰LV的压力—容积曲线如图3所示。
图3
压力容积曲线
Figure3.
Pressure-volume loop
整体上可以看出，两种模式压迫下的压力—容积曲线几乎重叠，说明两种模式的每搏输出量（stroke volume，SV）、EF、峰值压力（peak pressure，PP）和每搏功（stroke work，SW）对于未施加压迫的LV都有所提升，且增加量近似相等，两种压迫模式下的LV血流动力学参数具体值如表5所示。
表5
血流动力学参数
Table5.
Hemodynamic parameters
在收缩末期时心肌收缩力达到最大，此时心肌纤维主应力（以σff表示）分布如图4所示。施加压迫之后，位于心尖部的应力集中有明显改善，且整体的应力分布有所减小。
图4
收缩末期LV心肌纤维应力分布
Figure4.
Fiber stress distribution in the LV end-systole
2.2 考虑相位误差情况下两种辅助模式效果的对比
延迟两种模式压迫开始时间，血流动力学参数结果如图5所示，包括EF、PP和SV。从中可以看出， 随着压迫开始延迟时间的增加，模式一的各项血流动力学参数都有所下降，当延迟时间增大到50 ms时，EF下降了2%，PP下降了6 mm Hg，SV下降了4 mL，而模式二的各项血流动力学参数都基本保持稳定。
图5
血流动力学参数
Figure5.
Hemodynamic parameters
心尖扭转角是LV尖部相对于基部旋转的角度，在射血过程中，心内膜纤维间质变形产生势能，随后扭转变形[29]。模式一和模式二在0.00～0.20 s之间收缩过程中的心尖扭转角变化如图6所示。可以看出，当压迫开始时刻与射血开始时刻同步时，模式一在收缩末期的心尖扭转角为18.6°，模式二在收缩末期的心尖扭转角为20.3°，相差不大。随着压迫开始时间的滞后，模式一的心尖扭转角逐渐减小，在50 ms时为13.7°；模式二的心尖扭转角度依旧保持稳定。
图6
心尖扭转角
Figure6.
Heart apex twist angle
2.3 两种辅助模式对左心室生物力学的影响
不同开始时间点的LV收缩末期相对于心跳周期开始时刻的纤维应变分布如图7所示，图7中数据均为真实应变（以LE11表示）。从图7中可以看出，两种模式的压迫均可以减小LV收缩末期的最大纤维应变。未施加外膜辅助压迫时，最大纤维应变为14.6%。当压迫开始时间没有延迟时，模式一和模式二的最大纤维应变分别为10.0%和9.5%，近似相等。随着压迫开始时间的滞后，两种模式的心肌纤维应变均增大，但模式一增加的更多，当滞后时间为50 ms时，模式一的最大纤维应变为12.9%，模式二的最大纤维应变为10.2%。相对于模式一，改变压迫开始时间对模式二的最大心肌纤维应变几乎没有影响。
图7
不同压迫开始时间的纤维应变分布
Figure7.
Fiber strain distribution at different compression start time
3 讨论
本文针对心衰患者构建了个性化LV有限元模型，在此基础上分别施加了正弦曲线及维持一段时间PP的心外膜辅助模式。
对比来看，在血流动力学参数方面，两种辅助模式下的EF、PP和SV都有所上升，其中LV的EF分别从29.33 %提升到37.32 % 和37.77 %，可以维持患者的正常生理活动（EF降低型心衰的基线EF为40%）。
从生物力学的结果来看，两种模式的压迫都可以增加LV收缩能力及改善应力集中。收缩末期时LV处于收缩状态，心肌纤维长度小于静息长度，如图7所示，最大心肌纤维应变减小说明LV的收缩能力增强。从图4可以看出，由于心尖部心肌壁较薄，所以有明显的应力集中现象，长时间的应力集中可能导致局部心肌梗死，施加心外膜辅助可以很大程度改善应力集中。
本文认为有两个方面的原因导致了应力和应变的减小。第一，由于心肌纤维本身的结构特点，心肌纤维在心外膜处与周向成− 60°，在心内膜处为60°，心内壁区域下心肌纤维更接近于轴向方向，并且心肌纤维在沿着纤维方向上的硬度远远高于横向方向。在心外膜施加压迫时，心室被拉长，纤维有向心尖处伸展的趋势，处于收缩状态的心肌沿纤维方向的加载减小，因此心肌纤维的应力和应变都减小。第二，施加心外膜辅助可增强LV的扭转能力。如图6所示，施加模式一和模式二的压迫后，收缩末期心尖扭转角从8.0°增加到了18.6°和20.3°，恢复正常扭转能力，该数据参考Badano等[29]测量的健康人收缩末期心尖扭转角约为16°。当纤维达到收缩峰值时，沿纤维方向的应力会很大，且心内膜收缩力大于心外膜，但由于心外膜半径较大，心外膜处可以产生更大的扭矩以控制旋转方向。在射血期间，心外膜下心肌较大的扭矩经壁耦合至心肌中层和心内膜下心肌，导致LV整体旋转。在心外膜下，这种扭矩有助于主纤维方向的收缩；在中层，扭转动作增强了圆周方向的缩短；在心内膜下，扭转动作可引起肌纤维重新排列，使得心内膜下肌纤维沿纤维的法向方向更紧密地排列以增厚左室壁。扭转减小会增加心内膜应力和应变，从而增加需氧量，降低LV收缩功能的效率[29]，提高心室的扭转角有助于使纤维应力和肌纤维的缩短沿室壁均匀分布。
进一步推迟压迫开始时间以模拟相位误差的影响，根据血流动力学结果显示，压迫开始时间的滞后会削弱模式一辅助效果，EF、SV和PP都有所下降。具体而言，EF在压迫延迟50 ms时从37.32%下降到35.60%，相对于理想情况下的EF提升，效果降低。从图5可以看出，模式二的各项血流动力学结果基本保持稳定不受影响。这是由于模式二在整个射血期间都维持了压迫力峰值，尽管压迫开始时间滞后，在收缩末期时心外膜压迫力仍与心肌纤维收缩力同时处于峰值，因此对误差的抵抗能力较强。从生物力学角度来看，模式一的心尖扭转能力下降，在压迫延迟50 ms时从18.6°下降到13.7°，模式二基本维持在20.0°左右。模式一的心尖扭转能力减弱，心室收缩效率也降低，心肌纤维需要更大的变形来实现收缩动作，如图7所示，模式一的心肌的纤维最大应变从10.0%增大到12.9%，模式二的心肌纤维最大应变依旧保持不变，说明本文提出的加载模式在减小心肌应变方面能力更强，从而在引导心衰心肌的反向良性重塑方面更有优势。
本文的研究也存在一定的局限性：首先，本文仅考虑了LV的模型，而未对其他心腔的运动情况进行模拟。其次，所建立的简化体外循环系统耦合并不能完全模拟人体的血液循环系统。在下一步工作中，将考虑多个心腔的交互，搭建更为复杂的系统，同时根据不同程度的心衰情况建立不同类型的心脏模型，以进一步比较不同的心外膜辅助方式对心脏生物力学的影响，旨在为DCC的设计与优化、临床治疗和术后康复提供理论支持。
4 结论
本文所提出的维持一段时间峰值压迫力的模式和传统正弦加载压迫力模式相比，带来的血流动力学参数提升近似相等，并且心肌的生物力学结果也近似相同。然而，由于DCC工作周期和心脏跳动周期存在的相位误差，正弦加载压迫力模式的血流动力学参数有所下降，心室的收缩能力也有所降低。相比之下，本文所提出的维持峰值压迫力的模式，可以在误差为50 ms的情况下保持血流动力学参数及心肌纤维应力、应变不受影响，因此可对DCC设计和优化提供理论指导。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：李宸、徐俊波和张小刚构思方案，李宸实施方案并处理数据，张盛负责作图，强贤杰、王天博、刘晓翰、张悦、黄刚、徐俊波和靳忠民讨论并修改了论文。
