用于量化人冠状动脉体内机械特性的血管内超声基于图像的有限元建模方法

摘要

体内 电影血管内超声图像显示对应于不同压力负荷条件的冠状动脉横截面运动。基于有限元模型，采用迭代方案，通过匹配计算模型和医学图像中的冠状动脉运动来确定 体内 冠状动脉的患者特异性机械特性。

摘要

量化冠状动脉壁的机械特性可以为冠状动脉疾病的诊断、管理和治疗提供有意义的信息。由于需要持续监测的患者无法获得患者特异性冠状动脉样本，因此无法直接实验测试血管材料特性。当前的冠状动脉模型通常使用来自现有文献的材料参数，这会导致严重的机械应力/应变计算错误。在这里，我们将介绍一种基于有限元模型的更新方法 （FEMBUA），以根据医学图像量化冠状动脉的患者特异性 体内 材料特性。从冠状动脉疾病患者那里获得了冠状动脉的 体内 电影血管内超声 （IVUS） 和虚拟组织学 （VH）-IVUS 图像。分割显示一个心动周期内血管运动的 Cine IVUS 图像，并选择两个具有最大和最小管腔周长的 IVUS 框架分别代表收缩压和舒张压条件下的冠状动脉几何形状。还对 VH-IVUS 图像进行了分割以获得血管轮廓，并在 VH-IVUS 轮廓上增加了 0.05 cm 的层厚以重建冠状动脉几何形状。使用各向异性 Mooney-Rivlin 材料模型创建了一个计算有限元模型，该模型用于描述血管的机械特性和冠状动脉管腔表面规定的脉动血压条件，以使其收缩和膨胀。然后，采用迭代更新方法通过将计算有限元模型的最小和最大管腔周长与电影 IVUS 图像的最小和最大管腔周长进行匹配来确定各向异性 Mooney-Rivlin 模型的材料参数。这种基于图像的有限元模型更新方法可以成功扩展到确定各种血管床中动脉壁的材料特性，并具有心血管疾病风险评估的潜力。

引言

冠状动脉疾病 （CAD） 是导致死亡和发病的主要原因之一，2019 年全球死亡人数超过 914 万 ^1,2。冠状动脉疾病（如动脉粥样硬化和狭窄）的发展通常伴随着机械力的变化和血管壁材料特性的变化³。冠状动脉的材料特性不仅是决定其对生理负荷的机械反应的基石，也是模拟血管机械行为、预测动脉粥样硬化病变发展和评估各种医疗器械治疗效果的关键要素 ^4,5.因此，对冠状动脉材料特性的深刻理解和准确量化对于早期疾病诊断、精准医学和预后评估具有重要价值⁶。

离体冠状动脉组织的力学实验，如平面双轴测试、压痕测试、膨胀延伸和单轴延伸测试，是量化离体冠状动脉血管壁机械性能的常用方法 ^7,8,9。通过这些方法，从患者或实验动物那里获得冠状动脉样本。进行力学试验，确定不同应力条件下容器壁的应变响应，然后通过拟合实验数据确定材料参数¹⁰.先前的研究表明，冠状动脉特性是高度非线性和各向异性的¹¹。虽然离体实验可以提供准确的材料性能数据，但也存在重大局限性，具体如下：首先，从活体受试者中取出样品后的力学行为会与体内条件下的力学行为不同，这可能会影响测试结果的准确性。其次，由于伦理和实际限制，很难获得大量冠状动脉的正常或病理组织来进行机械测试。

为了克服这些限制，研究人员探索了用于体内、实时和患者特异性冠状动脉材料特性定量的新技术。其中，基于医学图像的基于有限元模型的更新方法 （FEMBUA） 有望解决这些具有挑战性的问题。这种方法利用血管内超声 （IVUS） 和虚拟组织学 （VH）-IVUS 等先进的成像技术来捕获详细的冠状动脉几何形状、组织成分及其运动¹²。通过构建 3D 有限元 （FE） 模型并结合患者特定的生理血压条件，可以通过优化材料参数来匹配图像数据，从而快速准确地量化冠状动脉材料特性，从而恢复心动周期期间的动态血管行为¹³。与离体实验相比，体内 FE 更新方法的优势包括无需组织切除的体内评估、促进大规模评估以及模拟复杂条件下的血管动力学以帮助了解冠状动脉疾病的病理生理学。

本文介绍了基于有限元模型的更新方法的关键步骤，包括电影 IVUS 和 VH-IVUS 图像的详细分割和处理、计算薄层结构纯模型的重建、迭代方案的执行以搜索冠状动脉组织的最佳材料参数。该协议的目的是使用 FEMBUA 方法作为演示来量化 CAD 样本患者的冠状动脉的材料特性，尤其是分步方法的插图。我们通过讨论这种 体内 方法的意义和其他方面来结束。

选定的参与者是一名 64 岁的女性，既往无冠状动脉疾病临床病史。该患者在出现胸痛症状后被诊断为冠状动脉疾病。进行冠状动脉造影和 IVUS 扫描以确认诊断。在左前降支中部发现一个狭窄 60% 的斑块病变。经评估，采用最佳药物治疗对患者进行治疗。

研究方案

在获得知情同意的情况下，从东南大学中大医院的一名 CAD 患者那里获得去识别化的临床数据，包括 体内 IVUS 图像和血压数据。样本患者是从中间冠状动脉粥样硬化病变临床研究的患者库中选出的，以展示量化患者特异性冠状动脉物质特性的方法¹⁴。该研究是按照东南大学中大医院临床研究伦理委员会批准的方案（批准号：2017ZDSYLL023-p01）进行的。

1. 数据采集和处理

1. Cine IVUS 和 VH-IVUS 图像采集
   1. 在冠状动脉造影的引导下，将 IVUS 导管放置在动脉粥样硬化病变的远端，并将其拉回近端。沿途生成灰度 IVUS 图像，以可视化冠状动脉血管的横截面。
   2. 使用配备高性能 20 MHz、2.9F 铂导管的 IVUS 成像系统，从 CAD 患者那里获取带有动脉粥样硬化斑块的冠状动脉段的 IVUS 图像（图 1）。
   3. 在 IVUS 图像采集期间，将导管在预先选择的病变部位暂停约 2 秒，以获得一系列称为电影 IVUS 的 IVUS 图像。电影 IVUS 图像清楚地显示了给定斑块部位在心动周期中的动态横截面变化。
   4. 根据心电图上 R 峰值瞬间采集的 IVUS 帧生成 VH-IVUS 图像，以使用成像系统以颜色编码形式可视化斑块成分。
      注意：VH-IVUS图像为动脉粥样硬化斑块中的四个关键斑块成分提供了直观的颜色图：红色的富含脂质的坏死核心（脂质），白色的钙化，深绿色的纤维组织和浅绿色的纤维脂肪组织。
   5. 以 DICOM 格式保存 VH-IVUS 和电影 IVUS 图像，以便离线分析。
2. 图像分割和处理
   1. 使用查看器打开 DICOM 文件，双击 相应的序列名称 以打开图像，然后单击 导出>导出图像 以将每个电影 IVUS 帧或 VH-IVUS 帧保存为 BMP 格式的单个图像。每个 BMP 图像都包含 500 x 500 像素，如原始 DICOM 文件所示。
   2. 逐帧检查电影 IVUS 图像，以查找在一个心动周期期间在预选斑块部位获得的连续帧。该样本斑块部位在一个心动周期内生成了 26 个电影 IVUS 帧。
   3. 检查所有生成的 VH-IVUS 图像，以找到在给定斑块部位获得的 VH-IVUS 图像。此处使用的 VH-IVUS 图像是在一次心动周期中使用电影 IVUS 帧中的一个 IVUS 帧创建的。
   4. 使用 ImageJ 软件分割 VH-IVUS 和电影 IVUS 图像，以获得血管边界和斑块成分边界的轮廓（参见 图 1C）。
      1. 选择 “直线”>“分割线 ”选项卡，并在电影 IVUS 和 VH-IVUS 图像上手动描绘管腔、冠状动脉外边界和斑块成分的轮廓。对于电影 IVUS 图像，仅分割管腔和外血管边界的轮廓，而对于 VH-IVUS 图像，分割管腔、外血管边界和斑块组件边界的轮廓。
      2. 为简单起见，仅保留大的脂质组分以生成有限元模型，而忽略小的、孤立的脂质组分。该噬菌斑样品中仅存在一种脂质。使用 Image > Overlay 将 描绘的轮廓叠加到原始图像上 > Add Selection 选项卡。
      3. 导航到 到 ROI 管理器 菜单以管理等值线、调整属性并将颜色和线宽设置为适当的值，以实现更好的可视化效果。选择 Properties 选项卡，将描边颜色设置为其他颜色，并在 width 处填充线宽。在这里，将管腔、外界和脂质等值线的线颜色分别设置为绿色、蓝色和红色，并将线宽设置为 3。
      4. 使用命令栏中的 Edit > Selection > Fit Spline 平滑轮廓 选择特定轮廓以获得平滑轮廓。此作利用样条曲线拟合技术自动平滑轮廓。
      5. 单击 “XY 坐标”选项卡中>“文件”>另存为 “，将每个轮廓（如流明、外边界和斑块组件）的点坐标保存在单独的 txt 文件中。此文件包含构成等值线的点的 x 和 y 坐标值，以像素为单位。
   5. 记录 DICOM 文件中电影 IVUS 和 VH-IVUS 图像中每个像素的真实物理大小（表示为像素大小）。此处使用的 IVUS 数据中一个像素的实际距离为 0.002 cm。此信息将用于将以像素为单位的点坐标转换为以 cm 为单位的实际距离。
3. 等值数据处理
   1. Cine IVUS 轮廓数据处理
      1. 使用 MATLAB 在一个心动周期中从所有电影 IVUS 图像中读取管腔轮廓的 txt 文件。
      2. 将所有 Lumen 轮廓乘以像素大小，以获得 Lumen 轮廓的实际大小。
      3. 计算所有管腔轮廓的管腔周长，并确定具有最大 （Cmax） 和最小 （Cmin） 管腔周长的 IVUS 框架，分别代表舒张压和收缩期。
   2. VH-IVUS 轮廓数据处理
      1. 使用 MATLAB 从 VH-IVUS 图像中读取管腔、外界和斑块成分轮廓的 txt 文件。
      2. 将所有等值线乘以像素大小，得到所有等值线的实际大小。
      3. 将每个等值线重新划分为 100 个等距点并执行 2D 平滑，以获得新的 VH-IVUS 等值线数据来替换旧的等值线数据。

2. 有限元模型

1. 冠状动脉血管几何重建
   1. 通过为 VH-IVUS 轮廓的所有点（包括管腔、外界和脂质）添加 z 坐标值，并在三维空间中创建一层轮廓，并为所有点设置 z = 0（图 2）。
   2. 通过为 VH-IVUS 等值线的所有点添加 z 坐标值来创建另一层等值线，并为所有点重置 z = 0.05 cm。
      注意：这两层等值线通过向 VH-IVUS 等值线添加固定的 0.05 cm 层厚，为仅薄层结构模型重建了 3D 冠状动脉血管几何形状（图 2）。
2. 有限元网格生成
   1. 通过线性插值每层权重为 1/3 和 2/3 的管腔和外边界轮廓（图 3A）来创建两个辅助轮廓。
   2. 通过将管腔/外边界连接到脂质轮廓上最近的点（例如，图 3B 中的点 A 和 B）或两个带有径向线的辅助轮廓，将血管区域分为 8 个圆周部分和 3 个径向部分（参见图 3B）。
   3. 用直线连接各层之间的所有点，形成具有 3 x 8 体积的 3D 结构（图 3B）。使用六面体单元划分每个体积，以生成有限元网格（图 3C）和不同的材料组（图 3D）。
   4. 通过将网格密度细化 10% 直到溶液的变化< 5% 来执行网格分析。
3. 材料特性定义
   1. 使用改进的各向异性 Mooney-Rivlin 材料模型来描述冠状动脉血管壁的材料属性。假设冠状动脉血管和斑块成分是超弹性、各向异性、几乎不可压缩和均匀的材料，改进的各向异性 Mooney-Rivlin 材料模型的应变能密度函数为：
      (1)
      (2)
      (3)
      其中 I₁ 和 I₂ 是右柯西-格林变形张量 C 的第一和第二不变量，定义为 c = [c_ij] = X^TX，X = [X_ij] = []，（X_i） 是当前位置 （a_j） 是原始位置，I₄ = c_ij（n_c）_i（n_c）_{j，n c} 是容器圆周方向的单位向量。c₁， c₂ 、 D₁ 、 D₂ 、 K₁ 和 K₂ 是患者特异性的材料参数。
   2. 根据离体双轴测试结果分配患者特异性冠状动脉的材料参数的初始值，即 c₁ = −1,312.9 kPa，c₂ = 114.7 kPa，D₁ = 629.7 kPa，D₂ = 2.0，K₁ = 35.9 kPa 和 K₂ = 23.5（图 4A-B）^13,15。
   3. 为板状组件分配材质参数（如果存在）。更具体地说，对于脂质，c₁=0.5 kPa，c₂=0，D₁=0.5 kPa，D2=1.5;对于钙化，我们使用 c1=920 kPa、c2=0、D1=360 kPa 和 D₂=2.0（图 4B）¹⁶。
      注：假设斑块成分（脂质和钙化）是超弹性、各向同性的且几乎不可压缩的，它们的机械性能由各向同性 Mooney-Rivlin 材料模型描述，公式 （2） 中给出了应变能密度函数。
4. 控制方程和边界条件设置
   1. 定义仅薄层结构模型的控制方程，其中包括运动方程、非线性 Cauchy-Green 应变-位移关系和冠状动脉血管材料模型¹¹。
   2. 在管腔表面规定患者特定的血压波形，以模拟真实的生理条件（图 4C）。为了获得患者特定的血压波形，请使用袖带测量的收缩压和舒张压值缩放典型的主动脉压波形（图 4D）。

3. 基于有限元模型的更新方法，用于患者特定的冠状动脉材料属性

注意：确定患者特异性冠状动脉材料特性的迭代过程如图 5 所示。

1. 通过以 95% 的固定收缩率轴向收缩从 VH-IVUS 图像重建的冠状动脉几何形状，并以最初设置为 98% 的圆周收缩（表示为 S）圆周收缩，确定对应于零压力条件的空载几何形状作为计算模型的初始几何形状。
   注意：由于从 VH-IVUS 图像重建的冠状动脉几何形状是在 体内 条件下，在管腔上规定血压，并且从栓系的远端和近端冠状动脉段轴向拉伸，因此 体内冠 状动脉几何形状应在圆周和轴向收缩以获得零压力几何形状。
2. 在以下步骤中，将轴向收缩率固定为 95%，并更新圆周收缩率。
3. 定义材料比率（表示为 k）以分配冠状动脉的患者特定材料特性：即 c₁ = k*（−1,312.9） kPa，c₂ = k*114.7 kPa，D₁ = k*629.7 kPa，K₁ = k*35.9 kPa，固定 D₂ = 2.0 和 K₂ = 23.5。
   注意：由于仅获得两个数据点（对应于舒张压和收缩压的最小和最大管腔周长）来确定未知参数（圆周收缩率 S 和 Mooney-Rivlin 模型的材料参数），我们通过假设冠状动脉的 体内 患者特异性材料特性与初始猜测成正比，材料比率表示为 k，从而减少了未知参数的数量： 也就是说，c₁ = k*（−1,312.9） kPa，c₂ = k*114.7 kPa，D₁ = k*629.7 kPa，K₁ = k*35.9 kPa，而 D₂ = 2.0 和 K₂ = 23.5 是固定的。
4. 在以下迭代过程中，将 k 值集更新为初始值 k 1 以及圆周收缩率 S。
5. 运行软件对计算模型进行求解，得到数值结果。
   1. 使用 MATLAB 将用于创建薄层仅结构模型的所有命令写入批处理文件（补充文件 1）中。
   2. 使用高级用户界面 （AUI） 加载此批处理文件以生成模型（图 6A）。通过单击 Data File/Solution 求解薄层仅结构模型，并将其保存为 .dat 文件（图 6C）。模拟三个心动周期，并采用最后一个周期的解来呈现数值结果。
   3. 通过导航 List > Value List > Zone ，然后在 Coordinate 下的 Variables to List 中选择 X-POSITION、Y-POSITION 和 Z-POSITION ，将节点坐标的结果导出到 txt 文件。单击 Apply and Export 以导出坐标结果。
   4. 将与舒张压和收缩压条件相对应的管腔轮廓数据保存到.txt文件中，以便计算管腔周长。
6. 将舒张压条件下由 FE 模型（仅薄层结构模型）计算的管腔周长与体内电影 IVUS 数据 （Cmin） 进行比较，并检查相对误差是否为 <1%。如果满足条件，则转到下一步，或者使用 secant 方法更新材料比率 k，然后转到步骤 3.3 以再次运行^17,18。
   注意：在第一次迭代中，使用 Newton 方法而不是割线方法更新材料比率。
7. 将 FE 模型在收缩压条件下计算的管腔周长与 体内 电影 IVUS 数据 （Cmax） 进行比较，并检查相对误差是否为 <1%。如果是，则停止迭代过程，或者更新圆周收缩率 S 并重定向到步骤 3.4 以再次运行。
   注意：在第一次迭代中，使用牛顿法而不是割线法来更新圆周收缩率。
8. 记录最佳 S 和 k 值，并计算 Mooney-Rivlin 材料模型的相应材料参数。
9. 绘制冠状动脉血管的圆周和轴向应力-拉伸比曲线（图 7），可以导出如下：
   (4)
   其中 σ 代表柯西应力，λ 代表拉伸比，i = c，a 代表圆周和轴向。
   1. 要在特定方向上绘制材质曲线，请将另一个方向的拉伸比率固定为 1。计算圆周和轴向的有效杨氏模量（分别标记为 YMc 和 YMa）作为拉伸比区间 [1.0， 1.1] 处材料曲线的刻度函数的斜率，以反映冠状动脉¹³ 的一般材料刚度：
      (5)
10. 随时提取斑块应力/应变分布，并记录收缩期和舒张期的节点分布和最大应力值（图 8）。

结果

我们详细描述了 FEMBUA 方法，该方法能够在实时 IVUS 成像后对冠状动脉斑块进行快速斑块材料和应力分析，并且可以确定斑块的体内材料特性和生物力学结果。表 1 提供了该冠状动脉的 Mooney-Rivlin 材料模型的体内材料参数。有限元模型的仿真结果（包括冠状动脉中的应力/应变分布）如图 8 所示。结果的详细分析如下。

讨论

协议中的关键步骤
基于有限元模型的更新方法中最关键的步骤在于迭代过程。在该方法中，有限元模型应从 体内 电影 IVUS 图像中准确恢复血管横截面上的冠状动脉血管运动。为此，本研究采用了最小化有限元模型和 体内 图像之间的管腔周长差异，以找到合适的材料特性。协议中还有其他关键步骤，包括图像分割和处理、有限元模型构建和?...

材料

Name	Company	Catalog Number	Comments
Bee DICOM Viewer	SinoUnion Healthcare Inc.	Version 3.5.1	A DICOM image reader software
ADINA	Adina R & D	Version 9.0	Finite element solver
ImageJ	National Institutes of Health		Segmented IVUS contours
MATLAB	MathWorks	Version R2018a	Commercial programming platform
Volcano s5 imaging system	Volcano Company		Intravascular ultrasound imaging system