第八章 颈动脉斑块应力分析：患者特异性建模问题
作者：Hao Gao（格拉斯哥大学数学与统计学院）、Quan Long（布鲁内尔大学生物工程研究所）
出处：L. Saba 等主编，《Multi-Modality Atherosclerosis Imaging and Diagnosis》，Springer，2014 年，第八章
第一节 章节概述
本章聚焦于颈动脉粥样硬化斑块应力分析中的患者特异性建模问题，旨在探讨如何将有限元方法与医学影像技术相结合，以实现对斑块易损性的准确评估。作者 Gao 和 Long 指出，获取准确的斑块内部应力场在生物力学界仍是一项重大挑战，尤其是涉及患者特异性几何形状、边界条件和材料模型等方面的工作。本章分为两大部分：第一部分基于高分辨率体内颈动脉 MRI，对一例有症状患者和一例无症状患者进行斑块应力分析；第二部分则利用各向异性材料本构模型进行斑块应力分析，该模型源自体内高分辨率颈动脉 MRI 数据。研究表明，有症状患者的斑块壁面最大应力显著高于无症状患者（227.7 kPa 对比 134.9 kPa），且有症状患者的最小纤维帽厚度更薄、脂质核心体积更大。本章的工作为斑块易损性评估提供了重要的生物力学依据，也为临床决策支持系统的开发奠定了基础。
第二节 关键问题与研究动机
颈动脉粥样硬化斑块破裂是导致缺血性脑卒中的主要机制之一。目前，斑块应力分析已从二维结构模型发展到三维流体固耦合模型，从理想化模型发展到患者特异性模型，但在实际应用中仍面临诸多挑战。首先，在几何建模方面，虽然高分辨率多光谱 MRI 技术能够可视化和重建斑块成分，但图像空间分辨率有限，纤维帽厚度可能低于 MRI 分辨率下限（本研究为 0.39 mm），导致纤维帽重建存在不确定性。其次，在材料模型方面，动脉壁具有非均匀、非线性、各向异性和粘弹性特征，缺乏准确的患者特异性材料参数数据是现有文献中最大的不确定性来源。第三，如何将患者特异性几何形状与患者特异性材料模型相结合，构建完整的患者特异性斑块力学模型，仍是一个尚未解决的关键问题。本研究的目标是基于体内 MRI 数据，建立一套从几何重建到材料参数识别再到应力分析的完整流程，为斑块易损性评估提供更为准确的生物力学预测工具。
第三节 主要公式与推导
3.1 Mooney-Rivlin 本构模型
动脉壁假设为非线性、各向同性、不可压缩材料，采用 ANSYS 11.0 中的三维非线性 Mooney-Rivlin 模型描述，其应变能密度函数 \(W\) 为：
其中 \(I_1\) 和 \(I_2\) 分别为第一和第二应变不变量，\(J\) 为材料变形体积与未变形体积之比，\(d\) 为材料不可压缩参数，\(C_{10}\)、\(C_{01}\)、\(C_{20}\)、\(C_{11}\)、\(C_{02}\) 为材料常数。本研究中采用的参数值为：\(C_{10} = 50.445\) kPa，\(C_{01} = 30.491\) kPa，\(C_{20} = 40\) kPa，\(C_{11} = 120\) kPa，\(C_{02} = 10\) kPa，\(d = 1.44 \times 10^{-7}\)。脂质核心则假定为软材料，弹性模量为 2 kPa，泊松比为 0.49。
3.2 动脉壁运动学
考虑厚壁圆筒假设下的 CCA 截面运动学，引入三种构型：\(\Omega_0\) 为零应力和切除后构型 \((R, \Theta, Z)\)，\(\Omega_1\) 为完整但无载荷构型 \((\lambda, \beta, \zeta)\)，\(\Omega_2\) 为体内有载荷构型 \((r, \theta, z)\)。通过 Humphrey 提出的关系式建立各构型之间的映射：
其中 \(t\) 为时间，\(\Theta_0\) 为动脉壁切除时的开角，\(\xi\) 为考虑从 \(\Omega_0\) 到 \(\Omega_1\) 轴向残余应力的轴向拉伸比，\(\phi\) 为体内载荷引起的轴向拉伸比（假设在心动周期中恒定）。
变形梯度张量定义为 \(F = \partial x / \partial X\)：
左、右 Cauchy-Green 张量为：
局部体积比为：
假设动脉壁不可压缩，则 \(J = 1\)。由上述关系可推导出径向方向的伸长比：
3.3 Gasser 应变能函数
采用 Gasser 等人提出的应变能函数描述具有双族胶原蛋白纤维增强的动脉壁组织：
其中 \(a\)、\(D\)、\(k_1\)、\(k_2\) 为材料参数，\(\varepsilon\) 表示纤维方向分散度（本研究取 \(\varepsilon = 0\)，表示纤维完全对齐）。\(N = 2\) 表示考虑两族纤维，\(I_1 = \text{trace}(C)\) 为第一应变不变量，\(I_{4(f,f)} = a_f C a_f^T\) 为伪不变量。
对于非线性弹性不可压缩材料，第二 Piola-Kirchhoff 应力为：
Cauchy 应力张量可由 \(S\) 计算得到：
3.4 动脉动力学平衡方程
在无体积力条件下，\(\Omega_2\) 构型下理想化 CCA 截面的运动方程为：
其中 \(\rho\) 为 CCA 截面的密度，\(a_r\) 为径向加速度。由于惯性项贡献很小，采用准静态假设。腔内压力可由下式计算：
其中 \(P_i(t)\) 为计算得到的腔内压力，\(P_a(t)\) 为外膜层对中膜层的压力类贡献（周向应力）。
第四节 关键算法与建模方法
4.1 MRI 图像采集与斑块成分分割
研究招募了一例有症状患者和一例无症状患者，均来自神经血管专科门诊。有症状患者近期发生过视网膜或皮层短暂性缺血发作，并在事件发生后 6 个月内接受扫描；无症状患者在成像前无任何症状。采用 1.5 T MRI 系统（GE Diagnostic Imaging）和四通道相控阵颈线圈（PACC）进行多对比度成像。扫描范围从颈总动脉分叉下方 12 mm 至狭窄远端 12 mm，覆盖颈动脉分叉的大部分区域。成像序列包括：中间 T2 加权脂肪抑制序列（ImT2W FatSat）、T2 加权序列（T2W）、短 T1 反转恢复序列（STIR）。视野为 \(10 \times 10\) cm\(^2\)，矩阵大小 \(256 \times 256\)，层厚 3 mm，像素尺寸为 \(0.39 \times 0.39 \times 3\) mm。
使用 MATLAB 编写的内部程序进行斑块分割，识别脂质核心、动脉壁和管腔。斑块区域基于 T2W、ImT2W FatSat 和 STIR 图像进行重建，正常动脉部分基于 TOF 图像重建。
4.2 流体固耦合有限元模拟
结构模型采用约 90,000 个十节点三维四面体单元进行网格划分，颈内动脉（ICA）和颈外动脉（ECA）出口平面的计算节点在所有方向上固定，并施加 11% 的轴向预拉伸。流体域采用约 1,000,000 个三维四面体单元进行网格划分，血液假定为不可压缩牛顿流体，粘度为 \(4 \times 10^{-3}\) Pa·s，密度为 \(1,067\) kg/m\(^3\)，流动假定为层流。瞬态模拟采用随时间变化的 CCA 入口压力和 ICA、ECA 出口质量流率作为边界条件。
4.3 材料参数拟合
基于相位对比 MRI 数据，通过非线性最小二乘法拟合材料参数。将轴向残余应变相关拉伸比 \(\xi\) 取为 1（文献值），体内载荷引起的轴向拉伸比 \(\phi\) 取为 1.1，开角 \(\Theta_0\) 取为 \(130^\circ\)，静水压力 \(P_a(t)\) 取为 2 kPa。中膜-内膜界面半径 \(r_m\) 初始值设为 0.6 mm，并利用动脉壁面积恒定（不可压缩性）约束计算心动周期中 \(r_m\) 的变化。最终拟合得到的材料参数为：\([a, k_1, k_2, \phi, R_m] = [32.1 \text{ kPa}, 9.6 \text{ kPa}, 3.8, 45.7^\circ, 4.8 \text{ mm}]\)。
第五节 主要结论
本章研究表明，患者特异性斑块应力分析对于评估斑块易损性具有重要价值。通过对有症状和无症状患者的对比分析，发现有症状患者的斑块壁面最大主应力显著高于无症状患者（227.7 kPa 对比 134.9 kPa），与 Tang 等人的研究结果（破裂斑块 \(247.3 \pm 121.4\) kPa 对比未破裂斑块 \(109.1 \pm 21.0\) kPa）高度一致。在形态学特征方面，有症状患者的最小纤维帽厚度（FCT）明显更薄（0.087 mm 对比 0.177 mm），脂质核心体积明显更大（247 mm\(^3\) 对比 70 mm\(^3\)）。FCT 与最大主应力之间存在显著的负相关性（有症状患者 \(r = -0.434, p < 0.05\)；无症状患者 \(r = -0.692, p < 0.05\)），表明薄纤维帽区域对应高应力区域，这与斑块破裂的生物力学机制相符。
基于相位对比 MRI 的患者特异性材料模型拟合方法成功应用于斑块应力分析，拟合压力与测量压力吻合良好，最大差异约为 0.288 mmHg。该方法考虑了动脉壁的各向异性超弹性特性和胶原蛋白纤维增强效应，能够更真实地反映体内动脉壁的力学行为。
第六节 挑战与开放问题
尽管本章在患者特异性斑块应力分析方面取得了重要进展，但仍存在若干挑战和未解决的问题。首先，MRI 图像分辨率仍然是限制因素。目前最高分辨率的体内 MRI 仍无法清晰显示厚度小于 0.39 mm 的纤维帽，这与临床观察到的临界破裂厚度（200 μm）存在显著差距。未来需要发展更高分辨率的 MRI 技术或采用其他成像方式（如光学相干断层扫描）来弥补这一不足。
其次，材料模型的区域变异性尚未解决。即使对于同一受试者，不同位置的动脉壁材料特性也可能存在差异，这增加了建立准确患者特异性模型的难度。第三，目前的分析主要关注静态或准静态应力分析，而斑块破裂是一个动态过程，涉及疲劳加载和时间依赖性材料行为。第四，本研究仅包含两例患者的小样本研究，需要积累更多不同斑块原型的数据来验证和完善该流程。
第五，临床转化仍面临诸多障碍，包括如何将生物力学分析结果与现有临床标准（如狭窄程度）相结合，以及如何实现自动化、标准化的分析流程。
第七节 个人评述与思考
本章在患者特异性斑块应力分析领域做出了系统性贡献，将医学影像、计算力学和生物力学有机结合，为斑块易损性评估提供了新的思路。作者采用了完整的端到端流程，从 MRI 图像采集、斑块分割、三维几何重建、流体固耦合有限元分析，到患者特异性材料参数识别，展示了该领域的最新研究水平。
笔者认为，本章的亮点在于第二部分提出的基于体内 MRI 的材料参数识别方法。该方法巧妙地利用相位对比 MRI 测得的管腔面积变化数据，结合非侵入性血压信息，通过优化算法反演材料参数，避免了传统方法中需要体外实验或组织样本的局限性。这一思路对于其他生物力学问题的患者特异性建模具有重要的借鉴意义。
然而，本章也存在一些值得进一步探讨的问题。例如，在材料模型方面，虽然 Gasser 模型考虑了胶原蛋白纤维的各向异性，但动脉壁的实际结构更为复杂，还涉及弹性纤维、蛋白多糖基质等多种成分，以及血管平滑肌的主动收缩能力。在 FSI 耦合方面，血液被视为牛顿流体，而实际血液具有非牛顿特性，尤其是在狭窄区域的高剪切率环境下。此外，残余应力的测量和建模也值得更深入的研究。
总体而言，本章为颈动脉斑块应力分析的患者特异性建模提供了重要的方法论框架。随着医学影像技术的进步和计算能力的提升，基于生物力学模拟的斑块易损性评估有望在临床决策中发挥越来越重要的作用，为精准医疗时代的心脑血管疾病防治提供有力工具。
公式汇总
| 编号 | 名称 | 形式 | 物理意义 | 类型 |
|---|---|---|---|---|
| (8.1) | Mooney-Rivlin 应变能密度函数 | \(W = C_{10}(I_1 - 3) + C_{01}(I_2 - 3) + \cdots + \frac{1}{d}(J - 1)^2\) | 描述动脉壁非线性、各向同性、不可压缩超弹性行为 | (T) |
| (8.4) | 变形梯度张量 | \(F = \frac{\partial r}{\partial R} e_r \otimes e_r + \frac{\lambda r}{\Theta_0 R} e_\theta \otimes e_\theta + \phi e_z \otimes e_z\) | 描述材料从参考构型到当前构型的变形 | (T) |
| (8.5) | 左、右 Cauchy-Green 张量 | \(B = FF^T, \quad C = F^TF\) | 描述变形的应变度量 | (T) |
| (8.6) | 局部体积比 | \(J = \det(F) = \lambda_r \lambda_\theta \lambda_z\) | 描述材料体积变化 | (T) |
| (8.9) | 径向伸长比 | \(\lambda_r = \frac{\Theta_0}{\phi \xi} \cdot \frac{\sqrt{r_m^2 - r^2}}{R_m}\) | 描述径向方向的拉伸程度 | (T) |
| (8.10) | Gasser 应变能函数 | \(W = a(I_1 - 3) + \frac{k_1}{2k_2}\sum_{f=1}^{N}(e^{k_2 E_f^2}-1) + \frac{1}{D}(\frac{J^2-1}{2}-\ln J)\) | 描述含双族胶原蛋白纤维增强的动脉壁各向异性超弹性 | (T) |
| (8.14) | 第二 Piola-Kirchhoff 应力 | \(S = -pC^{-1} + 2\frac{\partial W}{\partial C}\) | 描述非线性弹性不可压缩材料的应力 | (T) |
| (8.15) | Cauchy 应力张量 | \(\sigma = FSF^T J^{-1}\) | 描述真实物理应力状态 | (T) |
| (8.20) | 动脉动力学平衡方程 | \(\frac{\partial \sigma_{rr}}{\partial r} + \frac{\sigma_{rr} - \sigma_{\theta\theta}}{r} = \rho a_r\) | 描述径向力平衡 | (T) |
| (8.21) | 腔内压力计算公式 | \(P_i(t) = P_a(t) + \int_{r_i(t)}^{r_m(t)} \frac{\sigma_{\theta\theta}(r, t) - \sigma_{rr}(r, t)}{r} dr\) | 通过应力差积分计算腔内压力 | (T) |
注：(T) = 理论推导，(E) = 经验公式