基于3.0T MRA 影像数据的椎动脉开窗几何建模与血流动力学分析

王东侠，江洁，刘鹤，尹燕伟，胡明成★

（1.北京市丰台中西医结合医院 放射科，北京 10000；2.牡丹江林业中心医院 感染科，黑龙江 牡丹江 157011；3.牡丹江医学院附属红旗医院 磁共振科，黑龙江 牡丹江 157011）

0 前言

椎动脉开窗畸形颅内开窗畸形发生率中所占比例极低，据报道据相关文献报道为仅为1.3%-5.3%[1]。椎动脉开窗畸形表现为颅内椎动脉局限性血管重复，表现为椎动脉在走行过程中1 支分成两支，平行走行一段距离再回合[2]。颅内椎动脉开窗畸形由于其弯曲增多，导致血流动力学指标改变，严重损害动脉内膜，易再损伤局部形成血栓，从而导致动脉狭窄或脑栓塞[3]。随着医学影像技术的快速发展，对椎动脉开窗畸形的影像诊断[4]，3.0T MRA 相比1.5T MRI 更具有优势，椎动脉开窗畸形经 CEMRA 或 3D TOF MRA 处理可获得满意的图像。近几年来，应用计算流体力学（Computational fluid dynamics，CFD）方法研究血液系统的血液流变学越来越广泛，凭借CFD 技术的数值模拟血流动力学仿真，可求解得出各种血流动力学指标，对血液流变的形成规律、手术评价等医学研究、治疗提供借鉴和帮助。目前，应用MRA 影像数据构建椎动脉开窗模型与CFD 血流动力学仿真研究不同椎动脉开窗类型的血流动力学的文献极少。本文旨在根据不同几何形态的椎动脉开窗畸形患者的影像资料，构建三维几何模型和CFD仿真技术计算，以探讨椎动脉开窗几何形态对血流动力学指标的影响。

1 资料与方法

1.1 一般资料

检查方法与数据获取：采集核磁科椎动脉开窗畸形患者两例和健康者一例MRA 影像数据。健康患者排除其他脑血管疾病，椎动脉开窗患者经两名副高及以上的影像科医生诊断为右侧椎动脉开窗。裂隙型开窗患者，男，49 岁，主诉头晕伴颈部不适数月、近期明显加重，测血压为130/100mmHg(1mmHg=0.133kPa),心 率85 次/min；凸透镜型开窗患者，男，45 岁，主诉颈部不适，测血压为140/110mmHg，心率90 次/ 分；检查前去除患者随身金属异物，嘱患者扫描过程中平静呼吸、不得随意运动，保持仰卧位头先进，身体长轴与扫描床长轴一致，扫描线圈中心为下颌下缘，常规进行横断位及冠矢状位扫描，范围为听眶线至颅顶。使用飞利浦 Achieva3.0T 磁共振成像系统，采用梯度回波序列及小角度激励，以三维时间飞跃法及磁化转移对比技术重建血管，扫描图像数据统一保存为DICOM 格式。患者知情同意并签署知情同意书。

使用的软件及图形工作站：医学图像处理软件：MIMICS 21.0（比利时，Materialise 公司）；医学正向工程建模软件：3-Matic(比利时，Materialise公司）；网格划分软件：Fluent Meshing；CFD 仿真软件：CFX 2019R3(美国，ANSYS 公司）；图像后处理软件：EnSight 2019R3(美国，ANSYS 公司）；血流动力学指标求解程序采用ccl 格式编写；工作站：戴尔（DELL）T5820 图形工作站（CPU I9-10900X 10核3.7GHz，内存64G，硬盘512G+4T，显卡RTX3080 10G）。

1.2 方法

1.2.1 三维模型构建

使用 MIMICS 软件对DICOM 格式MRI 图像文件进行阈值分割，采用蒙版编辑工具和三维重建算法进行三维重建建立初步模型[5]。将三维重建获取的初步模型导入3-matic 软件中进行包裹、光顺、重构三角面片及切面等处理并以stl 格式导出。

1.2.2 有限元网格划分

将三个stl 格式模型导入Fluent Meshing 软件中进行四面体非结构化网格划分，为保证计算精确设定6 层边界层。

1.2.3 计算设置

由于椎动脉属于大血管，根据雷诺数确定血流模式为层流，不考虑重力影响，计算按牛顿流体计算，血液粘度系数为0.0035Pa·s,密度为1060kg/m3。边界条件为入口速度，入口速度-时间公式根据体积流量转化（式1、式2）[6-7]。

其中d 为血管直径，t 为时间，Vmean为进口平均速度，q 为入口流量。r 代表所在点位置的半径，v为入口速度[8]。在垂直于血流方向的血管截面上的血流速度分布为越接近入口中心处速度越快，越远离中心处速度越慢，即血流是充分发展的（图1），入口体积流量曲线如图2 所示，出口采用压力出口，出口曲线如图5 所示。

图1 充分发展的血流

图2 体积流量随时间变化曲线

2 结果

2.1 血流动力学参数选取

本研究选择时间平均壁面切应力（time average wall shear stress，TAWSS）、平均壁面切应力梯度（time average wall shear stress grade,TAWSSG）、剪切震荡系数（oscillatory shear index，OSI）、粒子滞留时间(relative retention time，RRT)等参数作为观察指标[7-9]。壁面切应力（wall shear stress,WSS）是指血液流动时在血管壁表面上引起的切向的动态摩擦力[9]。对于脉动流，在一个心脏周期内用每个节点上积分WSS 量的值来计算TAWSS：

其中wssi是瞬时剪切应力矢量，T 是周期的持续时间。同时提出了壁面切应力梯度（wall shear stress,WSSG），可以更明显观测WSS 数值变化，WSSG的时间平均值即平均壁面切应力梯度（TAWSSG）：

在一个心动周期内，震荡剪切系数（OSI）可以描述WSS 方向的变化程度，OSI 值介于（0，0.5）之间，OSI 数值越大[10]，表示WSS 方向的变化也就越大，但其值的大小与WSS 没有必然联系，其表达式为：

2.2 OSI 云图分析

由图3 所示，正常与开窗畸形椎动脉的在椎动脉主干血管上OSI 呈现不同的变化趋势。凸透镜型开窗畸形椎动脉在椎动脉末端分叉区域出现高OSI，相较于裂隙型高OSI 区域面积较小，说明裂隙型开窗畸形椎动脉的血液振荡区域小。无论是椎动脉裂隙型开窗畸形还是凸透镜型开窗畸形的椎动脉的高OSI 区域都比正常的椎动脉高OSI 区域面积大得多，由此可见开窗畸形造成WSS 方向变化剧烈。

图3 不同开窗畸形椎动脉与正常椎动脉OSI 分布云图

2.3 TAWSS 云图分析

图4 为不同几何形态的椎动脉开窗TAWSS 分布云图。由图可见，在一个心动周期内，裂隙型开窗畸形椎动脉在椎动脉主干血管壁的低TAWSS 区域面积最大，其次为凸透镜型开窗椎动脉，而正常椎动脉主干血管壁TAWSS 最高。TAWSS 的分布大小具有临床意义，TAWSS 如果低于4dynes/cm2，则提示主干血管容易出现动脉硬化，需要临床干预。

图4 不同开窗畸形椎动脉与正常椎动脉TAWSS 分布云图

2.4 TAWSSG 云图分析

图5 为TAWSSG 分布云图。TAWSSG 云图提示在椎动脉开窗区域，存在TAWSSG 局部区域较高情况，推测可能由于血流经过开窗畸形部位，血管较细，导致局部血流速度较快，对血管壁造成剧烈冲击，形成高TAWSSG 区域。TAWSSG 在正常椎动脉主干血管分叉处呈现出数值较其他区域较大，但TAWSSG 过大，则易造成血管损伤。

图5 不同开窗畸形椎动脉与正常椎动脉TAWSSG 分布云图

3 讨论

在一个心动周期内使用每个节点WSS 量值进行积分即为TAWSS，TAWSS 为一个心动周期内血管壁面的平均WSS 大小。高WSS 能导致内皮细胞功能紊乱致使内皮细胞损伤，主要通过活性保护因子分泌减少及分泌血管损伤因子增多[9-16]；并且血管内皮细胞平行血流方向排列且倾向于变长是在WSS 大于1N/m2的条件下，内皮细胞功能失调且内皮细胞增加渗透性，出现淤滞且停留时间显著延长的情况。在综合情况共同影响下，血管壁周边的内皮细胞吸收边流中的致动脉粥样硬化因子增多形成斑块等病变。为了更加清晰地表示出动脉粥样硬化易发的高风险区域，使用时均参数TAWSS 进行分析，从上述云图可看出，从三种不同几何形态学的椎动脉开窗中，裂隙型具有较大的TAWSS 区域，并且高TAWSS 区域与高速涡流血流区域及高OSI 区域有部分共同区域，这也提示该区域有较大风险形成血管壁的粥样硬化。OSI 是依据一个心动周期内动脉血管壁上不同位置WSS 方向的变化计算出来的，表示的是一个心动周期内WSS 方向振荡的情况，当OSI 值较大时即表明该周期内WSS 振荡情况越严重；同时OSI 也可表示血液流动方向的情况，当OSI 值越大时即血液流动方向可能发生频繁的方向改变甚至有可能回流[12]。当OSI 值为0 时可见在该心动周期内血管内皮细胞所受剪切方向恒定，具有较低的粥样硬化斑块形成的风险[14]；而当OSI 值不为0 时，即该OSI 值区域内的内皮细胞所受的剪切力有振荡[16-18]；数值越大，越易导致内皮细胞功能障碍，从而极大促使斑块的发生。裂隙型和凸透镜型均于开窗分支血管终末段连接处呈现高OSI 区域；重复型开窗的高OSI 主要集中于开窗分支起始处及开窗外侧支拐角小片状区域，该区域内出现粥样硬化性斑块损伤的可能性较大。