低温–机械耦合效应下的生物组织损伤有限元模型

Finite Element Model of Biological Tissue Damage under Hypothermia-Mechanical Coupling Effect

摘要: 冷冻消融术作为目前治疗肿瘤的重要方法之一，具有安全、高效、微创等优势。然而，脑肿瘤冷冻消融术中引发的颅内压波动现象未得到足够的重视。研究颅内温度与压力相互作用的机制，或许有助于解决这一问题。文章考虑温度和围压的耦合作用，构建了肿瘤冻结后机械损伤的有限元模型，评估低温环境下生物组织冻结后力学性能的影响及损伤规律。并采用COMSOL数值模拟胶质瘤低温消融中颅内温度–力学耦合的响应机制。

Abstract: Cryoablation, as one of the important methods for the treatment of tumors, has advantages such as safety, efficiency, and minimal invasiveness. However, the phenomenon of intracranial pressure fluctuation caused by brain tumor cryoablation has not received enough attention. Studying the mechanism of the interaction between intracranial temperature and pressure may help to solve this problem. Considering the coupling effect of temperature and confining pressure, a finite element model of tumor mechanical damage after freezing was constructed to evaluate the effect of freezing on mechanical properties and the damage law of biological tissues in a low-temperature environment. COMSOL was used to simulate the response mechanism of intracranial temperature-mechanical coupling during glioma cryoablation.

1. 引言

颅内压增高(Increased Intracranial Pressure, ICP)是神经外科常见的临床病理综合征[1]，是颅脑损伤、脑肿瘤、脑出血、脑积水和颅内炎症等所共有的征象。由于上述疾病使颅腔内容物体积增加，导致颅内压持续在2.0 kPa以上，即为颅内压增高。颅内压增高会引发脑疝危象，可使病人因呼吸循环衰竭而死亡，因此对颅内压增高的及时诊断和正确处理十分重要。针对脑肿瘤引起颅内压增高的现象，冷冻消融更适合于脑肿瘤的治疗。该手术在MRI引导下将探针精准插入肿瘤体内，利用冷却剂液化氮使探针加入−150℃ [2]的温度冻结部分肿瘤组织，从而达到消除肿瘤的目的。颅缝闭合后的颅腔容积是固定的，并且颅腔含有脑组织、脑脊液和血液三种内容物。在消融过程中，受到温度影响内容物体积改变，因此颅内压会产生变化。这一现象表明了研究颅内温度–力耦合反应机制对临床手术具有重要的参考价值。

在冷冻消融术中，生物组织处于低于−150℃的环境下，会在短时间内形成冰球，其与低温下的水凝胶[3] [4]、冰、岩石等具有相似的力学性质。然而，关于组织冻结后，温度如何影响颅内压以及这两者的耦合作用机制，还需要做进一步的研究。目前，国内外尚未对脑肿瘤消融过程中温度与颅内压的关系进行深入研究。为了更好地评估消融过程中颅内被冻结组织的机械损伤特点，部分研究者对水凝胶的力学特性进行了深入研究。Lei [5]等揭示了仿生黏弹性水凝胶材料的力学生物学响应机制，以及其调控细胞、组织行为的规律。现阶段，已经研发出具有仿生组织特性的水凝胶，可以应用于消融治疗的相关研究。Tao [6]等通过理论模型和有限元分析研究具有多种变形模式的水凝胶膨胀变形力学，包括正负膨胀、各向同性和各向异性、梯度和弯曲变形模式，得到了水化过程中组织变形与带隙变化的对应关系。这些有关仿生水凝胶温度或力学条件下力学性能的研究，为生物组织在温度–力耦合作用下的性能响应提供了理论依据。Wang [7]等对肺癌的冷冻疗法构建了基于实际肺部肿瘤的3D模型，采用热机械耦合模拟热损伤和应力损伤来描述损伤评估。说明了温度–力耦合条件下组织损伤可使用有限元模拟进行评估。

颅内力学性能相关的有限元模拟和试验大多是基于脑爆炸损伤和车祸等脑意外损伤的力学损伤机制，模型材料参数通常来源于尸体解剖数据。目前的研究状况表明，对于在低温条件下消融手术中颅脑内容物的温度和力学响应问题尚未得到应有的重视，而它对于推动临床手术中患者的颅内压和组织变形研究具有重要的价值。评估手术中冷冻脑组织形成冰球后的持续损伤阶段，需要一个包含组织机械行为的数学模型来表达脑组织中传热损伤的力学规律。同时，颅内力学损伤的应力应变关系是分析颅内压动态变化的基本前提[8]。在生物组织冷冻形成冰球的过程中，其应力会随着组织变形而产生变化，引发有关于大变形的弹性力学特性。深入研究温度和力学的耦合关系，包括热传递、力的传导以及组织的机械变形等方面，对于有效治疗肿瘤和促进人体健康的恢复具有至关重要的作用[9]。

本研究旨在探讨脑胶质瘤冷冻消融过程中，颅内内容物冻结部分的温度损伤和力学响应机制。冷冻组织的损伤研究主要内容包括温度损伤机制、力学响应机制[10]。因此，采用有限元法[11] (Finite Element Method, FEM)对脑胶质瘤冷冻治疗手术进行数值模拟，分析脑胶质瘤膨胀对颅脑各部位产生的压力数值和变形情况。结果表明，所建立的三维有限元模型能够较好地反映组织冻结后在持续低温损伤和机械损伤耦合下力学性能的变化规律，为脑肿瘤冷冻消融手术的临床可行性预测提供了一定的推动意义。

2. 脑胶质瘤消融数值模拟

本文开发了一种用于模拟脑胶质瘤冷冻消融的模型，采用有限元法求解生物热方程和机械变形方程的耦合效应，并利用COMSOL Multiphysics实现。研究冷冻消融过程中颅内热传递和力学的耦合，并获得应力与应变之间的关系。

2.1. 数学方程

生物组织的热传递采用Pennes生物热方程进行描述[7]，该方程能够以最简化的方式准确捕捉组织中的温度分布，定义如下：

$\rho c\frac{\partial T}{\partial t}=\nabla \left(\text{k}\nabla \text{T}\right)+{\omega }_b{c}_b\left(T_b-T\right)+Q_m+Q_r$ (1)

其中，$T$ 、$\rho$ 、$c$ 、$\text{k}$ 分别为温度、密度、比热容、热传导率。单位分别为$℃$ 、$\text{kg}/{\text{m}}^3$ 、$\text{J}/\left(\text{kg}\cdot ℃\right)$ 、$\text{W}/\left(\text{m}\cdot ℃\right)$ 。${\omega }_b$ 为血液灌注率，$c_b$ 为血液比热容，$T_b$ 为血液温度，$Q_m$ 为组织代谢产热率($\text{W}/{\text{m}}^3$ )，$Q_r$ 为外部热源项($\text{W}/{\text{m}}^3$ )。

固体力学物理场仿真基于质量守恒和欧拉公式，并根据柯西应力张量求解，动力平衡方程表达为：

$\rho \frac{{\partial }^{2}u}{\partial t^2}=\nabla \cdot S+F_V$ (2)

其中，$\rho$ 表示已变形组织的真实密度，$S$ 为应力，体积力$F_V$ 为单位变形体积上的力。

仿真模型主要包括生物传热和固体力学两个数学模型，对力学和传热两个物理场进行热膨胀耦合，求解得到脑组织冷冻消融过程中的温度分布、损伤程度、应力分布和应变分布。

2.2. 结果与讨论

模拟结果得到了冷冻过程中不同时间点脑内容物的温度分布，消融结束时治疗区域形成了一个5 mm的冰球。提取距离冻结源不同位置的7个点的温度数据，图1显示了这7个位置随时间的温度变化，并根据细胞死亡温度来判断其消融的准确性。同时，利用这些温度来求解脑组织机械行为的力学变形方程，然后利用得到的脑组织变形来重新计算颅内压值。最后获得脑肿瘤冷冻后的应力–应变关系。

通过模拟出的温度变化情况与细胞受损死亡的温度结合，确定出合理的热源、消融参数和边界条件。脑组织冷冻形成冰球后，冰球持续低温，随着脑实质的体积增大，而引起颅内压呈弥漫性增高。腔内各部位及各分腔之间压力均匀升高，不存在明显的压力差，因此脑组织无明显移位。

对比距离冷冻源7个不同位置下的脑组织应力分布(见图2)，我们发现：在脑组织区域，离冷冻源越近，则温度越低，同时应力较小且均匀分布。颅骨和脑组织的边界处，应力最大，且向外逐渐减小。同时，根据组织应变随时间变化的情况(见图3)，在组织的弹性阶段，距离冷冻源不同位置的应变非常接近，且都非常小。这也意味着组织功能区挤压变形很小，保证了消融手术的安全性。在87 s左右，应力和应变都达到了峰值点，此时治疗区域基本完全处于冻结状态。此后的消融中，治疗区域内的温度呈小幅下降，应力和应变的波动幅度也越来越小，逐渐趋于稳定状态。结合温度分布图，在5 mm半径肿瘤消融结束时，距离肿瘤冷冻源中心5 mm处的温度为269 K (−4℃)，距离肿瘤3~4 mm处冰体温度达到217~252 K (−56℃~−21℃)。

Figure 1. Temperature profile of the brain as a function of cryoablation time

图1. 颅脑随冷冻消融时间变化的温度分布图

Figure 2. Maps of stress distributions over time at different locations

图2. 不同位置的应力随时间变化分布图

Figure 3. Plots of strain distributions over time at different locations

图3. 不同位置的应变随时间变化分布图

结合图1~3可以看出，脑组织在温度100 K左右时，应力和应变的波动较为平稳，而90~100 K之间应力和应变的波动幅度明显，此处处在冻结温度的临界点。

3. 分析与讨论

Figure 4. Stress-strain plots at different positions

图4. 不同位置下的应力–应变关系图

通过应力–应变曲线的分析，研究在消融过程中颅内压与组织位移之间的关系。根据图4可知，在消融结束后距离冷冻中心不同位置的7个点对应的应力–应变关系，曲线的形状反映了肿瘤冰体在外力(颅内压)作用下发生的弹性形变过程。从图中可以看出，在距离冷冻源1 mm~7 mm的范围内，应变会随着应力的增大而增大。当距离冷冻源1 mm至4 mm时，距离越大，相同应变对应的应力就越大，而相同应力下产生的应变越小。然而，在距离冷冻源5 mm以外，随着距离越远，相同应变对应的应力越小，而相同应力下产生的应变也越大。这种现象的出现是由于生物组织的温度变化导致了弹性模量的改变，图4中生物组织冻结后不同位置下的曲线反映了不同温度下的应力–应变关系，距离冷冻源越近，温度越低。在距离冷冻源1~4 mm处弹性模量随温度降低而减小，在距离冷冻源5 mm外弹性模量则随着温度的降低而增大。同时，细胞内水分从液态到固态的转化程度也会对生物组织的应力–应变关系产生影响。在距离冷冻源1~4 mm处，组织基本已经冻结，并且达到了组织坏死温度，对生物组织的破坏力较强。而在距离冷冻源5 mm以外，为了保护正常组织不被损伤，温度较为安全，仅有少部分细胞被冻结。通过以上分析可知，应力–应变的分析适用于生物组织冻结后的损伤评估。

观察冷冻后的生物组织在持续低温下的体积应力–体积应变(见图5(a))和等效偏应力–等效偏应变(见图5(b))随时间的变化趋势，可以更直观地了解生物组织应力与应变之间的关系。体积应变是指组织在低温产生的载荷作用下组织的体积变化值与原体积的比值，可以看出在100 s时体积的比值达到了峰值，这也对应了生物组织在温度为100 s内迅速降温冻结从而产生的极速膨胀，并在之后趋于小幅度降温。等效偏应变则是组织在受力作用下产生的变形程度，同样在100 s内处于变形程度的峰值。这两种力学属性从体积变化和变形程度两个方向描述了生物组织的弹性性质。

Figure 5. (a) Trends in volumetric stress and volumetric strain over time; (b) Trend of equivalent deviational stress and strain over time

图5. (a) 体积应力和体积应变随时间的变化趋势；(b) 等效偏应力和等效偏应变随时间的变化趋势

数值模拟结果可以为本构关系的构建提供力学参数与理论参考。并利用仿真得到的应力–应变关系，与本构模型的理论结果进行对比验证，目的是建立一个用于生物组织冷冻后的机械损伤评价体系。

4. 力学特性

利用“固体传热”与“固体力学”接口耦合实现温度与力学响应的计算，这是研究损伤机理的基础。关于脑组织冻结后力学性质的描述包括以下内容：1) 黏弹性[12]：脑组织的黏弹性表现为应力–应变曲线关系、滞后关系、蠕变特征、应力松弛等。2) 各向异性[6] [7] [13]：由于大脑组织是由神经纤维有序排列组成的，具有方向性。脑组织在不同方向反应出来的不同力学特性被称为各向异性，并且排列方向不同会对大脑组织的力学性能造成影响。脑组织的力学特性非常复杂，本文作为脑肿瘤消融模拟的初步研究，将模型简化为各向同性材料。模拟结果为之后更进一步的研究做了理论铺垫。

大脑中的颅骨、脑脊液、大脑皮质和脑组织有着不同的材料属性和力学参数[13]。颅骨属于刚性较大、变形较小的材料，这里选用刚性材料表示。脑组织主要由水组成，质地柔软，体模量较大，其变形特性取决于剪切模量，呈现非线性和时间变化属性，且与血液间有相互作用，其材料特性采用黏弹性本构来表示软组织(白质、灰质)。脑脊液在大脑内部循环流动，存在于脑组织与颅骨之间的间隙，使颅骨与脑组织结构之间紧密连接，其材料特性采用黏弹性本构来表征。本文颅脑有限元模型的材料参数参考了国内外相关研究人员的实验研究数据[14]-[17]，见表1、表2。

张[18]等根据人体脑组织振动试验确定了剪切模量函数及其参数，并且通过应力–应变曲线的斜率获得了弹性模量。本研究中脑组织的密度和体积模量参数取自其指定的最精确的值，表1列出了灰质和白质的黏弹性力学参数。

Table 1. Mechanical parameters of viscoelastic materials of white matter and gray matter

表1. 白质与灰质的黏弹性材料力学参数

	短期剪切模量(KPa)	长期剪切模量(KPa)	衰变常数(s^−^1)
灰质	41.0	7.8	40
白质	34.0	6.4	40

由于白质和灰质的各项力学参数没有明显的差异，所以在计算过程中不对其进行区分。取脑组织的短期剪切模量为36 KPa，长期弹性模量为7.2 KPa。对于小变形响应，在线性应变状态下可采用广义Maxwell模型[19]描述松弛效应，剪切模量表达方程式用Prony级数表示，如下所示[20]：

$G\left(t\right)=G_{\infty }+{\displaystyle \sum }_{i=1}^n{G}_i\text{exp}\left(-t/{\tau }_i\right)$ (3)

其中，$G_{\infty }为$ 长期剪切模量，$G_i$ 为剪切模量，${\tau }_i$ 为松弛时间。

Table 2. Finite element material parameters and properties of the head

表2. 头部有限元材料参数及属性

组织	弹性模量 $E$ (MPa)	泊松比 $\nu$	短期剪切模量 $G_0$ (MPa)	长期剪切模量 $G_{\infty }$ (MPa)	衰变常数 $\beta$ (s^−^1)	材料属性
颅骨	6500	0.45	—	—	—	弹性
脑脊液	0.2	0.45	1.0 × 10^−^3	0.9 × 10^−^3	80	黏弹性
大脑皮层	16.7	0.45	—	—	—	黏弹性
脑组织	22	0.45	7.0 × 10^−^3	1.4 × 10^−^3	80	黏弹性
脑肿瘤	22	0.45	7.0 × 10^−^3	1.4 × 10^−^3	80	黏弹性
血液	0.2	0.45	1.0 × 10^−^3	0.9 × 10^−^3	80	黏弹性

通过确定大脑各部分的力学特性参数，实现对颅脑低温动态损伤的有限元模型构建。

由于简化模型忽略了脑组织的各向异性和非线性特性，可能造成温度以及应力–应变变化有些许差异，这将在一定程度上影响模拟的准确度。未来研究应着力去补足材料力学属性，使模型更贴近真实脑部形态和特性。

5. 结论

本文首先进行了肿瘤冷冻消融的数值模拟，通过模拟生物组织在结冰后持续低温下的力学演变特性，对现行的冷冻消融技术进行了力学上的评估体系建立。构建了肿瘤冻结后持续在低温与围压耦合作用下的有限元模型，从而深入研究组织破坏过程中的损伤演化特性。

在肿瘤冷冻消融术中，在持续低体温和颅内压的环境下，治疗区冷冻组织的机械损伤特性是提高临床手术精准度需要关注的研究课题。此外，温度和力并不是简单的单一叠加作用下产生的损伤，而是一种复杂的耦合效应。

此外，本文根据脑组织结冰后的持续低温损伤模拟，研究了脑组织冰体在瞬态压缩下的力学性能，揭示了肿瘤冰体压缩应力–应变曲线具有较强的应变率相关性和敏感性。并建立了一个用于评估冷冻消融手术中生物组织冻结后产生力学损伤机制的体系。

该模型的开发有助于推动实现医疗工作者在脑肿瘤冷冻消融手术中的精准预测，将肿瘤结冰后的颅内压与变形程度形成关联，从而开发出更有效的治疗机制，建立更加科学的治疗手段。后续研究应致力于改进颅脑消融手术，提供更准确的理论依据。