A Magnetic Resonance Imaging-based Computational Protocol for Analysis of Plaque Morphology and Hemodynamics in Patients with Carotid Artery Stenosis

Summary

Assessment of internal carotid artery (ICA) stenosis is based on the estimation of percentage stenosis, which does not account for physiologically relevant risk factors for stroke such as plaque composition and hemodynamics. This protocol leverages quantitative magnetic resonance imaging and computational fluid dynamics to characterize ICA plaque composition and hemodynamics.

Abstract

Current assessment and management of internal carotid artery (ICA) stenosis is based on the estimation of percentage stenosis via duplex ultrasound (DUS) or computed tomography angiography (CTA), which does not account for physiologically relevant risk factors for stroke, such as plaque vulnerability and hemodynamics. Knowledge of the composition of the carotid plaque and hemodynamic loads on the plaque can be used to provide a much more complete assessment of the embolic potential of the plaque rather than using percentage stenosis alone. Through pairing magnetic resonance imaging (MRI) and patient-specific computational fluid dynamics (CFD), differences in both hemodynamics across an ICA stenosis and plaque composition can be identified. Quantitative multi-contrast atherosclerosis characterization (qMatch) MRI allows for detailed analysis of plaque composition. CFD models can be created using phase contrast (PC) MRI, which can be used to obtain flow waveforms and CTA and/or time-of-flight (TOF)-MRI anatomy. After creating a 3D geometric model of the carotid bifurcation, PC-MRI derived waveforms are prescribed to the common carotid artery inflow and external carotid artery outflow. A three-element Windkessel model, which is iteratively tuned to match the patient's blood pressure, is then prescribed to the ICA. Finally, solutions to the incompressible Navier-Stokes equations are obtained to provide high-resolution velocity and pressure and thus capture hemodynamics across the carotid bifurcation and ICA stenosis. This article provides a detailed protocol that allows for non-invasive and patient-specific characterization of plaque composition and hemodynamic loads of patients with ICA stenosis.

Introduction

Internal carotid artery (ICA) stenosis is a leading cause of stroke, long-term disability, and death^{1,2,3,4,5,6,7}. Current assessment and management of ICA stenosis is based on the estimation of percentage stenosis via duplex ultrasound (DUS) velocities or cross-sectional anatomy [computed tomography angiogram (CTA) and/or magnetic resonance imaging (MRI)]. However, percentage stenosis does not account for physiologically relevant risk factors for stroke such as plaque vulnerability and hemodynamic loads across the plaque^{8,9,10,11,12,13,14}. Although reduced stroke risk following carotid endarterectomy (CEA) has been demonstrated in symptomatic patients with greater than 50% stenosis, the benefit of CEA in asymptomatic patients is debated^3,4. In fact, many surgeons reserve operative intervention for those with stenotic lesions >80% and/or in cases with high-risk (vulnerable) plaque morphology¹⁵. Improved methods of determining which ICA stenoses are at risk of plaque embolism and thus would benefit from CEA are warranted.

Quantitative multi-contrast atherosclerosis characterization (qMatch) is an MRI technique that utilizes low-rank modeling to enable high-resolution 3D imaging that provides co-registered multi-contrast dark-blood and bright blood images, and relaxometry images for comprehensive and quantitative assessment of carotid arterial plaques^16,17. qMatch has improved 3D isotropic resolution, large anatomic coverage, and quantitative assessment of carotid artery plaque burden compared to conventional MRI. Patient-specific computational fluid dynamics (CFD) can be used to characterize the hemodynamic loads on the plaque, thus providing unique information regarding the hemodynamic and biomechanical risk of cerebrovascular embolic events^{18,19,20,21,22,23}. Knowledge of the composition of the carotid plaque and hemodynamic loads on the plaque could be used to provide a more comprehensive assessment of embolic potential than percentage stenosis alone. In this work, we present a protocol that uses both qMatch MRI and MRI-informed CFD to identify differences in plaque composition and hemodynamics across an ICA stenosis.

Protocol

The study was approved by the Institutional Review Board of the University of Michigan and informed consent from each study subject was obtained. This protocol uses CRIMSON, a validated, open-source computational hemodynamics framework that performs key computational modeling tasks such as mesh generation, boundary condition specification, and finite element analysis^24,25. To download CRIMSON and/or review modeling tutorials, visit the website (https://crimson.software). The CRIMSON GUI requires a Windows operating system. The CRIMSON flow solver is available for both Windows and Linux.

1. Patient Recruitment and Patient-specific data acquisition

1. Recruit adult patients who have a diagnosis of severe ICA stenosis demonstrated on DUS and/or CTA (as defined by the North American Symptomatic Carotid Endarterectomy Trial (NASCET) criteria)²⁶. Include patients who do not have a known contraindication to MRI (i.e., metallic implants) or MRI intolerance (i.e., claustrophobia, inability to lie flat/remain still). Exclude patients if they are pregnant or have a contraindication to MRI. Obtain informed consent, which should include a discussion and understanding of the study procedure, risks, benefits, assurances of confidentiality, duration of study, and right to withdraw from study.
2. Obtain retrospective and/or prospective patient data to inform CFD models. Use CTA, MRI, and/or angiography images for patient anatomy.
   NOTE: Boundary conditions will be discussed in more detail later. However, in general, data to inform boundary conditions often include non-invasive or invasive pressure, DUS velocities, and/or phase contrast (PC)-MRI-derived flow.
3. Prior to MRI, conduct a detailed prescreening MRI safety form for each enrolled patient to identify any contraindications to MRI. Review MRI safety forms with 2+ study team members. Instruct enrolled subjects to remove all metallic items and provide them with a gown.
4. Position the subject supine on a 3T MRI system, provide hearing protection and blanket for patient comfort, and position a head and neck coil.
5. After performing initial localizing sequences to establish proper orientation over the carotid bifurcation, perform the following three sequences:
   1. Obtain a 3D time-of-flight MRI of the head and neck for anatomic characterization of the vasculature from the common carotid artery (CCA) at C5 to the distal ICA the foramen magnum.
   2. Obtain a 2D cardiac-gated PC-MRI at the level of the CCA at C5 and above the carotid bifurcation at the proximal external carotid artery (ECA) and mid ICA distal to the lesion to measure volumetric blood flow waveforms. Patient-specific velocity encoding (Venc) is based on the peak systolic velocity (PSV) at each vessel (CCA, ECA, and mid ICA distal to the lesion) measured via DUS. In general, aim for a Venc ~20% higher than that of the PSV at the vessel of interest.
   3. Use qMatch MRI sequence localized over the carotid bifurcation to obtain detailed information on plaque composition and plaque vulnerability.

2. Obtaining flow waveforms from PC-MRI

1. After obtaining the 2D cardiac-gated PC-MRI at the above locations, obtain volumetric flow waveforms using the built-in software on MRI scanner.
   1. On the MRI scanner, identify and use respective flow quantification software to obtain PC-MRI-derived flow waveforms.
   2. Select each vessel of interest (i.e., CCA, ECA, and ICA) and place a contour around the specified vessel to provide an automated flow waveform. Manually edit the contours to ensure the accurate area of the vessel.
   3. Export the flow waveforms from the respective software.
      NOTE: Flow quantification software may differ between different MRI manufacturers.
2. Utilize a Fourier Transformation to interpolate and create a flow waveform that is smooth, continuous, and has a larger number of data points, thus allowing for a more refined flow profile for CFD simulations.
   NOTE: imposing flow into CRIMSON²⁴ (which will be discussed later) it is important that the waveform function is continuous: both the function itself and its derivatives exist and are continuous for all values of time. The Fourier interpolation generates a continuous waveform based on any arbitrary combination of measured flow data points (PC-MRI) and desired time points (for CFD analysis).
3. To ensure conservation of mass between inflow and outflow faces, compare the average flow of the CCA, ECA, and ICA after Fourier Transformation.
   1. In cases where conservation of mass (i.e., CCA flow = ECA flow + ICA flow) is not within 10% do not move forward and proceed with troubleshooting.
   2. First, verify that an accurate was PSV used for Venc and check if the ECA PC-MRI-derived flow waveform was measured after a large branch (or branches).
   3. In cases where the PC-MRI flow waveform was obtained after large ECA branches,increase the flow to the ECA and re-check conservation of mass.

3. Computational fluid dynamics modeling: geometry

1. Import de-identified DICOM image data for patient-specific anatomy (CTA, MRI, angiography) into CRIMSON using the import button in the data manager.
2. Use the Geometry Modeling window to select Vessel Path Editing and create a vessel tree consisting of the anatomic range of interest (CCA, ECA, and ICA).
3. Use the Vessel Path Editing window to place centerline points along the length of each vessel in the anatomy of interest (CCA, ECA, and ICA).
   1. The CCA centerline is typically started at the level of C5, corresponding to the location where the flow waveform from PC-MRI was obtained.
   2. The ICA centerline is typically ended 1–2 cm distal to the stenosis, corresponding to the location where the flow waveform from PC-MRI was obtained.
   3. The ECA the centerline is typically ended proximal to the first-order branches off the ECA, corresponding to the location where the flow waveform from PC-MRI was obtained.
4. Using the Vessel Re-slice window, the length of centerline points along each vessel is visualized. This window will appear after at least two points along the vessel centerline have been added and contains a cross-sectional view along (perpendicular to) the centerline.
   NOTE: Vessel centerlines may also be imported in CRIMSON (they must be in the VTK file format).
5. Use the Vessel Re-slice window to specify the boundaries of the vessel wall by adding vessel contours (using a circle, ellipse, or manual contour). The Vessel Re-slice window provides a view of the vessel along the centerline so that accurate contours can be defined. Contours are added manually by the user across varying points of the vessel centerline in the Vessel Re-slice window.
   NOTE: On the left-hand side of the vessel re-slice window the original image is displayed. On the right-hand side of the vessel re-slice window the gradient of the image is displayed. The gradient image view can be helpful when defining contours, since it may show the boundary of the lumen more clearly.
   1. Place contours frequently enough along the centerline to fully capture the curvature and changing geometry of the vessel while not too close to overfit or produce artifacts.
6. After contours have been placed across the vessels of interest, use the Loft button in the Vessel Contour Modeling window to create a combined 3D solid model of each geometry via a process known as lofting.
7. Select the Vessel Blending window to generate a single solid geometry vessel. The most common algorithm for blending is the fillet. The typical fillet size is between 0.3 to 1 mm.

4. Computational fluid dynamics modeling: meshing

1. Select the Meshing and Solver Setup window and use the meshing button to visualize meshing options and select specific mesh parameters.
   NOTE: A mesh consists of multiple tetrahedral elements and is required to run a simulation as the Navier-Stokes equations for velocity and pressure are solved at each point (node) across the mesh. A basic mesh can be defined using global and/or local features. Specifically, the mesh can be defined by element size (i.e., a smaller element size leads to a smaller or more refined mesh), curvature refinement (which adds a more mesh elements to areas with higher curvature), or other local mesh refinement features. Specific meshing strategies may differ based on different geometries of interest. In the setting of the present geometry of interest (i.e., the CCA, proximal ICA, and proximal ECA), utilize both global and local mesh features.
2. Use the global options window to set the global element size to be an absolute value ranging between 0.5 mm and 0.75 mm.
3. Use the global options window to specify the boundary layer type as geometric growth. Set the total number of layers to be 3, the first layer thickness to be 0.2 mm, and the total layer thickness to be 1.0 mm, thus allowing for a finer mesh along the outside of the face and a less fine mesh along the middle of the face.
4. Lastly, use a curvature refinement to add more mesh elements at areas with curvature (i.e., at the stenosis).
   NOTE: Local mesh refinement options can also be used to create a finer mesh at specific vessels, bifurcation areas, or inlet/outlet faces.
5. Review the mesh elements by clicking on the Mesh Information button after right-clicking on the mesh.
   NOTE: A final mesh should contain elements with appropriate aspect ratios (ratio of the largest side to the smallest side of a given tetrahedral element, smaller is better), a distribution of elements that capture flow features in critical areas (i.e. the stenosis, vessel outlets, boundary layers), and avoid excessive distortion or sharp changes in cell size.
   Final meshes of the present geometry of interest should contain 400,000-700,000 elements.
   Figure 1A depicts the critical steps pertaining to patient geometry and meshing.

5. Computational fluid dynamics modeling: boundary conditions

1. To specify boundary conditions, select the Meshing and Solver Setup window and then select the Solver Setup icon. In the Solver Setup window, add a boundary condition set (referred to as a "BC Set") and then select a specific boundary condition using the BC icon.
   NOTE: Boundary conditions are used to represent the pressure and blood flow beyond the boundaries of the segmented model. Decision of which boundary conditions to use and where they are prescribed is arguably the most important and critical aspect of any CFD model and should be made deliberately and be supported by physiologically relevant meaning. Boundary conditions should be selected and tuned to match patient-specific values, and in cases where patient-specific values are not available literature data can be used to inform the computational model.
2. Observe the boundary conditions that are currently available in CRIMSON:
   1. Inlet: Pressure, Prescribed Velocity (Flow Waveform), custom lumped parameter circuit (any arbitrary combination of resistors, capacitors, inductors, pressure nodes, and custom circuit elements defined via a Python script).
   2. Wall: No slip (refers to a rigid or non-deformable wall), Deformable.
   3. Outlet: Pressure, RCR, Prescribed Velocity (Flow Waveform), custom lumped-parameter circuit.
3. Click the BC icon to select a specific boundary condition. First, select No Slip to implement rigid, non-deformable walls and apply this to all walls using the Apply to all walls button.
4. Next, click the BC icon and select prescribed velocity to import the previously defined inflow waveform (i.e., the PC-MRI derived CCA flow after Fourier Transformation). In the boundary condition window, map the parabolic velocity profile to the inlet of the CCA.
   NOTE: In CRIMSON, the convention is for inlet flows to be negative and outlet flows to be positive.
5. Similarly, import the pulsatile ECA outflow waveform (prescribed velocity) reconstructed from PC-MRI and map the parabolic velocity profile to the outlet of the ECA.
6. Select the BC icon | RCR to populate a three-element Windkessel model (RCR), which consists of a proximal resistance (R_p), a distal resistance (R_d), and a capacitor (C). Map the RCR to the outlet of the ICA. Calculate approximate patient-specific RCR values using the PC-MRI flow data and the patients' blood pressure.
   1. The total arterial resistance is R_T= P_mean/Q_T, where the mean blood pressure P_mean = 1/3 P_systolic + 2/3 P_diastolic, and Q_T is total cardiac flow entering the model (in this case CCA flow).
   2. The total arterial compliance is C_T = (Q_T,max-Q_T,min)/(P_systolic-P_diastolic)*Δt, where Q_T,max and Q_T,min are maximum and minimum values of CCA inflow, and Δt is the time lapse between these values.
   3. Initial estimates for the Windkessel model parameters are informed by patient-specific imaging and are obtained by distributing a fraction of R_T and C_T on the ICA outlet.
      NOTE: Figure 1B depicts the boundary conditions used in the present modeling scheme. The present study utilizes the aforementioned boundary condition set; however, other boundary condition sets could be utilized.  

6. Computational fluid dynamics modeling: simulation

1. In the Meshing and Solver Setup window, select the Solver Setup icon | Sovler Parameters to specify the solver parameters within CRIMSON.
   1. Run simulations using a time step size of 0.1 ms for four cardiac cycles.
      NOTE: The residual required for a solution to be considered converged for each time step is 1 x 10^-4. Because high-grade ICA stenoses have regions of complex and re-circulatory flow model blood as an incompressible Non-Newtonian fluid using the Carreau-Yasuda model. This can be done by adding a viscosity constant model to the solver input file (see 6.3.1). Set the density of blood to be 1,060 kg·m⁻³.
      A stabilized finite-element formulation for the incompressible Navier-Stokes equations solves for blood flow velocity and pressures in the models.
2. To start a simulation, prepare simulation files using Solver Setup in CRIMSON. Specifically, generate files containing the flow data (bct.dat), inlet flow at each time step (bctFlowWaveform.dat), information on the mesh and boundary conditions (geombc.dat), information for the face on which each boundary condition is applied (faceinfo.dat), the first time step number of the simulation (numstart.dat), 3-element Windkessel data (rcrt.dat), files containing information on pressure and velocity at every point in the mesh (restart files), and the instructions for the flowsolver (solver.inp).
   1. Add the Carreau-Yasuda model into the solver.inp and add to the simulation files to allow for blood to be modeled as a Non-Newtonian fluid.
3. To run simulations, choose from the following:
   1. For the simplest way to run the CRIMSON Navier-Stokes flowsolver, press the Run Simulation button in the Study pane of the Solver Setup window. This will open up a command window, which allows the user to specify how many processors to use.
      NOTE: The flowsolver can also be run from the command line using a Windows batch file.
      Although some simulations (i.e., those under the steady state assumption) can be run directly through CRIMSON on a local Windows desktop computer, pulsatile simulations with a mesh consisting of many tetrahedral elements (>200,000) will require a more computer high-performance computing (HPC) cluster with a Linux operating system.
4. Use the CRIMSON Navier-Stokes flowsolver to perform computations with 72–108 cores on an HPC cluster. If performing simulations on an HPC cluster, transfer all presolver files to the cluster.
   NOTE: The process of transferring files to an HPC cluster will differ for each individual and institution based on the technology and software that are available to them.
5. When the solver starts running, observe that an output file named "histor.dat" is printed in the command line. The simulation output files will be saved in a new directory called "n-procs-case" where "n" is the number of processors for the simulation.
   1. Use the linux prompt: tail -f histor.dat to view the "histor.dat" file in real time. The histor.dat file consists of multiple columns; however, the first four columns are the most important.
      1. Observe that the first column is the current time step, which may appear multiple times because within each step, the Navier-Stokes equations are solved multiple times to increase the accuracy of the numerical solution before proceeding to the next step (i.e., approaching the specified residual).
      2. Observe that the second column is the elapsed simulation time in seconds.
      3. Note that the third column is the non-linear residual, which is a measure of the quality of the current solution (a lower number indicates an improved solution).
      4. Observe that the fourth column is the logarithmic value of the current residual compared to the initial residual at the start of the simulation, which provides a measure of the current residual relative to the starting points.

7. Computational Fluid Dynamics Modeling: Post-Processing

1. Check for convergence after the simulation has completed (i.e., did the simulation approach succeed in meeting the specified residual). Use the information contained in the "histor.dat" file to plot and/or visualize the residuals.
2. To visualize detailed simulation results, post-processing is required. Navigate to the "n-procs-case" folder and perform the postsolver and multipostsolver executables (both of which can be found in the CRIMSON flowsolver install files).
   1. Use the postsolver executable (postsolver -sn <last time step> -td -ph -ybar) to generate a "ybar" file, which contains a measure of errors for each node in the mesh.
   2. Use the multipostsolver executable (multipostsolver <first time step> <last time step> <increment> <foldername>) to combine the restart files by retrieving the results at the specified increment between the first and last specified time steps.
3. Verify mass conservation by inspecting the "FlowHist.dat" file, which contains the flow waveforms for inflow of CCA, and outflow of ECA and ICA.
4. Inspect the pressure, including maximum pressure (SBP), minimum pressure (DBP), MAP, and pulse pressure (Pulse Pressure = SBP-DBP), by looking at the "PressHist.dat" file.
5. Adjust the RCR parameters to ensure agreement with patient-specific information such as blood pressure. Specifically, the resistance and capacitance are adjusted such that the simulated pulse pressure at the CCA outlet is within 5% of the patient's pulse pressure and 10% of MAP (from cuff measurement).
   NOTE: Increase the resistance to increase pressure (SBP, DBP, and MAP) and increase the compliance to decrease the pulse pressure (vice versa). Tuning the RCR is an iterative process, which is often referred to as fixed-point iteration.

8. Computational Fluid Dynamics Modeling: Data Analysis

1. After a simulation has passed the designated tuning (i.e., simulated pulse pressure within 5% of the patient's pulse pressure), export, visualize, and analyze the data.
2. Identify the "view.pht" file in the folder that was created after the multipostsolver executable was run [i.e. (multipostsolver <first time step> <last time step> <increment> <foldername>)] and import it into Paraview.
3. Calculate and visualize the following variables in Paraview.
   1. Velocity (and Flow): CRIMSON reports velocity in mm/s, however in DUS imaging velocity is reported in cm/s. Convert velocity to cm/s using a calculator in Paraview.
      NOTE: Because a no-slip or rigid wall boundary condition was used the velocity at the wall will be zero. Thus, it is best to visualize velocity using a volume render technique.
      1. Capture the velocity profile of a specific portion of the model (i.e. the point of maximum stenosis) using a clip or slice function in Paraview.
   2. Pressure (and Pressure Ratio)
      NOTE: CRIMSON reports pressure in Pascal (Pa); however, clinically, pressure is reported in mmHg. Convert pressure to mmHg using a calculator by dividing the pressure (in Pa) by 133.33.
      1. Use a clip or slice to capture the pressure proximal to and distal to the ICA stenosis. Use the "Plot Data Overt Time" filter in Paraview to obtain a pressure waveform over time (analogous to the waveform one would obtain with an invasive pressure measurement).
      2. Calculate the pressure ratio by dividing the mean distal pressure by the mean proximal pressure.
   3. Wall Shear Stress (WSS): Calculate the time-averaged WSS by first selecting the area of interest (the ICA stenosis), using a calculator to obtain the Magnitude of the WSS, and using the "Temporal Statistics Filter."
   4. Calculate the Oscillatory Shear Index (OSI) in Paraview after the time-averaged WSS is calculated (see above).
      NOTE: OSI is a measurement of how much the WSS changes direction and magnitude during a cardiac cycle. OSI values range from 0 to 0.5, where 0 indicates unidirectional WSS, and 0.5 indicates WSS with a time average of zero.

9. Analysis of plaque morphology using qMatch MRI

1. Run raw qMatch image data through MATLAB image reconstruction program to obtain post-processed images, including dark blood, T1-weighted, T2-weighted, MRA, qMatch T1 Map, and qMatch T2 Map images.
2. Use a DICOM viewer to visualize the post-processed qMatch MRI images and assess the plaque composition.
   NOTE: qMatch can identify plaque components, including calcium, intraplaque hemorrhage (IPH), lipid-rich necrotic cores (LRNC), and fibrous cap thickness and its status.
   1. In general, each component will have the following characteristics on qMatch datasets (Table 1).
      1. Calcium: Hypo-intense on dark blood, T1-weighted, and T2-weighted images.
      2. Recent IPH: Hyper-intense on T1-weighted and hyper- to iso-intense on T2-weighted images.
      3. Old IPH: Hyper-intense on T1-weighted and hypo- to iso-intense on T2-weighted images.
      4. LRNC: Hyper-intense on T1-weighted and Hypo-intense on T2-weighted images.
      5. FC: Hyper- to iso-intense on T2-weighted images.
3. Grade plaques, based on their components, using the modified American Heart Association²⁷ and/or Plaque-RADS (Reporting and Data System)²⁸ classification systems.

Results

Using this MRI-informed CFD workflow paired with qMatch MRI allows for identification of the hemodynamic loads across and ICA stenosis and the specific components of the plaque. We first begin with ensuring that we have a high-quality mesh to allow for accurate representation of flow features in critical areas. A final mesh should contain an adequate number of mesh elements with low aspect ratios (Figure 1A). A coarse mesh with high aspect ratios will likely lead to inaccurate simulation results. We then move forward with the specification of our boundary conditions (Figure 1B). After successful completion of the simulation and appropriate boundary condition tuning, non-invasive and patient-specific hemodynamics can be collected.

Specific hemodynamic metrics that can be measured, including but not limited to velocity, flow, pressure (including pressure ratios and pressure gradients), WSS, and OSI. Figure 2 shows a representative velocity profile across the carotid bifurcation and ICA stenosis. Visualization of the maximum velocity profile throughout the cardiac cycle can serve as a surrogate for a DUS-derived velocity waveform. Thus, both PSV and end-diastolic velocity (EDV) can be approximated. Figure 3 shows two representative examples of the pressure (mmHg) across the carotid bifurcation and ICA stenosis. A pressure gradient can be measured by collecting pressure waveforms proximal and distal to the stenosis.

In Figure 3A, there is minimal to no difference in pressure proximal to (red line) and distal to (blue line) the stenosis. However, in Figure 3B, there is a large difference in pressure proximal to (red line) and distal to (blue line) the stenosis. Figure 4 shows two representative examples of the WSS (Pa) mapped across the carotid bifurcation and ICA stenosis. In Figure 4A, there is a low WSS across the stenosis, whereas in Figure 4B, there is a large WSS across the stenosis. Figure 5 shows a comparison of OSI mapped across the carotid bifurcation before (Figure 5A: pre-operative) and after (Figure 5B: post-operative) CEA. Post-operative maps depict areas of higher OSI compared to pre-operative.

After appropriate postprocessing of the qMatch images, a dataset with six sets of DICOMs will be generated, including dark blood, T1-weighted, T2-weighted, MRA, qMatch T1 Map, and qMatch T2 Map sequences. Using these data sets, plaque components including calcium, IPH, LRNC, and fibrous cap thickness and/or rupture can be visualized and quantified (using the T1 map and T2 map sequences). Table 1 depicts the general characteristics of each plaque component on qMatch datasets. Figure 6 shows a representative qMatch dataset from a patient with IPH. The outline of the ICA is depicted with a solid white line, while the flow lumen is depicted with the dashed white line, and the plaque is depicted with the dashed yellow line. Features of IPH (solid red line) demonstrated by hyperintense signal in the T1-weighted image and lowered T1 measurement in the T1 map. Figure 7 shows a representative qMatch dataset from a patient with heavily calcified plaque. The outline of the ICA is depicted with solid white line, while the flow lumen is depicted with the dashed white line. Calcified portion of the plaque (dashed orange line) demonstrated by hypo-intense signal in the dark blood, T1-weighted, and T2-weighted images.

Figure 1: Overview of computational fluid dynamics modeling method. (A) Creation of patient-specific geometry and meshing as well as (B) specification of boundary conditions. (A) De-identified DICOM image data from CTA is imported into CRIMSON, and the anatomy of interest (including the CCA, ICA, and ECA) is determined. Centerline points are placed along the length of each vessel within the anatomy of interest. The boundaries of the vessel wall are specified by adding contours. Vessel branches are lofted, then combined with a fillet operation. The final geometric model is then discretized into a mesh, consisting of multiple tetrahedral elements with local mesh refinement at the level of the stenosis. (B) A 3-element Windkessel is prescribed to the ICA outlet to allow for variations in pressure and velocity. 2D cardiac-gated PC-MRI is obtained at the level of the CCA at C5 (red circle and ellipse) and above the carotid bifurcation at the proximal ECA (orange circle and ellipse) and mid ICA distal to the lesion (blue circle and ellipse) to measure volumetric blood flow waveforms. A flow waveform is prescribed to the CCA inlet and the ECA outlet. Abbreviations: CTA = computed tomography angiography; CCA = common carotid artery; ICA = internal carotid artery; ECA = external carotid artery; PC = Phase Contrast. Please click here to view a larger version of this figure.

Figure 2: Velocity information from CFD workflow. Right) Velocity (cm/s) mapped to a model of a carotid bifurcation including the CCA, ECA, and ICA with a severe stenosis in the anterior view. Left) Maximum velocity over time for one cardiac cycle can be visualized, serving as a surrogate for duplex ultrasound. Abbreviations: CCA = common carotid artery; ECA = external carotid artery; ICA = internal carotid artery. Please click here to view a larger version of this figure.

Figure 3: Representative example of pressure (mmHg) mapped across the carotid bifurcation for two cases in the anterior view. Pressure is mapped to geometric models of the CCA, ECA, and ICA. (A) Case with minimal to no difference in pressure proximal to (red line, red pressure waveform) and distal to (blue line, blue pressure waveform) the ICA stenosis. (B) Case with large difference in pressure proximal to (red line, red pressure waveform) and distal to (blue line, blue pressure waveform) the ICA stenosis. Abbreviations: CCA = common carotid artery; ECA = external carotid artery; ICA = internal carotid artery. Please click here to view a larger version of this figure.

Figure 4: Representative example of wall shear stress (Pa) mapped across the carotid bifurcation for two cases in the anterior view. WSS is mapped to geometric models of the CCA, ECA, and ICA. (A) Case with low WSS across the ICA stenosis. (B) Case with large WSS across the ICA stenosis. Abbreviations: WSS = wall shear stress; CCA = common carotid artery; ECA = external carotid artery; ICA = internal carotid artery. Please click here to view a larger version of this figure.

Figure 5: Comparison of oscillatory shear index before (preoperative) and after (postoperative) carotid endarterectomy, including both anterior and posterior views. OSI is mapped to geometric models of the CCA, ECA, and ICA. Lesion and repaired lesion (segments where OSI are compared) are highlighted. Post-operatively maps depict areas of higher OSI compared to pre-operatively. Abbreviations: OSI = oscillatory shear index; CCA = common carotid artery; ECA = external carotid artery; ICA = internal carotid artery; CEA = carotid endarterectomy. Please click here to view a larger version of this figure.

Figure 6: A representative qMatch dataset from a patient with intraplaque hemorrhage. (A) Dark blood, (B) T1-weighted, (C) T2-weighted, (D) MRA, (E) qMatch T1 Map, and (F) qMatch T2 Map sequences. The outline of the ICA is depicted with a solid white line, while the flow lumen is depicted with the dashed white line, and the plaque is depicted with the dashed yellow line. Features of IPH (solid red line) demonstrated by hyperintense signal in the T1-weighted image and lowered T1 measurement in the T1 map. Abbreviations: IPH = intraplaque hemorrhage. Please click here to view a larger version of this figure.

Figure 7: A representative qMatch dataset from a patient with calcified plaque. (A) Dark blood, (B) T1-weighted, (C) T2-weighted, (D) MRA, (E) qMatch T1 Map, and (F) qMatch T2 Map sequences. The outline of the ICA is depicted with solid white line, while the flow lumen is depicted with the dashed white line. Calcified portion of the plaque (dashed orange line) demonstrated by hypo-intense signal in the dark blood, T1-weighted, and T2-weighted images. Please click here to view a larger version of this figure.

Plaque Component	MRA	Dark Blood	T1w	T2w	T1-map	T2-map
IPH	+		+		Used for Quantification	Used for Quantification
Calcium		-	-	-	Used for Quantification	Used for Quantification
LRNC	=			-	Used for Quantification	Used for Quantification
Fibrous Cap	-/=		-/=	-	Used for Quantification	Used for Quantification

Table 1: Characteristics of plaque components on qMatch datasets. Abbreviations: MRA = Magnetic resonance angiography; T1w = T1 weighted; T2w = T2 weighted; IPH = intraplaque hemorrhage; LRNC = lipid-rich necrotic core; + = hyper-intense; - = hypo-intense; (=) iso-intense.

Discussion

Here, we presented a protocol to non-invasively characterize the hemodynamic loads and plaque composition across an ICA stenosis, thus providing a more comprehensive assessment of embolic potential than current diagnostic modalities which assess percentage stenosis alone. We begin by obtaining patient imaging and pressure data in both a retrospective and prospective manner, including CTA, PC-MRI, and blood pressure cuff data to inform our CFD models. Moreover, we tune the boundary conditions in our model, specifically the Windkessel model, to correspond to known patient data. As such, this protocol allows for the collection of precise and patient-specific data pertaining to physiologically relevant risk factors for plaque embolism and stroke.

Informing models and boundary conditions with physiologically accurate and patient-specific data is critical for accurate simulation results. Some computational simulations in the cerebrovascular space rely on DUS, numerical methods, or non-patient specific assumptions to derive inflow waveforms^21,29,30,31. The use of DUS is appealing given that it is widely available, frequently used in clinical settings, has lower cost, and is readily accessible. However, PC-MRI is generally considered to be a more accurate method for measuring flow^32,33,34. PC-MRI can directly quantify velocity at multiple locations within the lumen, thus accommodating for asymmetries within the flow field within a vessel and thus provides a more comprehensive depiction of flow dynamics^32,33. PC-MRI is also not subject to operator-specific biases introduced by DUS such as angle of interrogation and measurement selection location. On the other hand, DUS is often operator-dependent and is less precise in capturing vessel area and complex flow patterns often leading to inaccurate flows. Nonetheless, PC-MRI flow measurements are not perfect with an approximate 10% error^35,36. Special attention should be made to ensure appropriate vessel encoding, maintaining an image plane that is orthogonal to the vessel axial, appropriate temporal and special resolutions, and minimization of phase offset errors³⁷. Lastly, MRI may overestimate stenosis compared to CTA which should be considered when assessing patient geometry³⁸. Future work, focused on comparing hemodynamic outputs from CFD models informed by DUS-flow waveforms and those informed by PC-MRI flow waveforms are warranted.

The choice of outflow boundary conditions can have a significant influence on velocity and pressure fields in CFD simulations of blood flow. In our approach, we elected to impose a parabolic outflow waveform to the ECA and couple the ICA to a three-element Windkessel model. This approach for boundary condition specification allows for robust enforcement of conservation of mass between inlet and outlets, while also enabling accurate matching of the patient's blood pressure³⁹. Thus, we felt that this would provide the most accurate depiction of ICA hemodynamics. However, given that we are imposing a flow waveform to one of our model outlets (i.e., the ECA) is important to ensure that the outflow waveform is synchronized to the CCA inflow waveform³⁹. In our approach this was made possible by the collection of our flow data from 2D-cardiac gated PC-MRI. However, in cases where obtaining such data is impractical, a different boundary condition approach may be advantageous (i.e., coupling both the ECA and ICA to three-element Windkessel models) so that assumptions do not have to be made in the temporal alignment of inflow and outflow waveforms³⁹.

There are important limitations of this protocol to keep in mind. First, because this modeling approach only consists of the ipsilateral carotid bifurcation, it does not include the circle of Willis and/or important factors that impact cerebral hemodynamics, such as the presence of collaterals or the extent of contralateral ICA stenosis. Patients with incomplete collateral pathways in the circle of Willis have been demonstrated to have higher rates of severe stroke and worse prognosis after stroke^40,41,42. Furthermore, the presence of patent collaterals has been associated with a reduced risk of stroke and transient ischemic attack^9,43,44. Additionally, several studies have demonstrated that the presence of a contralateral ICA stenosis (or occlusion) impacts the ipsilateral ICA velocities^45,46,47,48. Moreover, our group has recently demonstrated that severe contralateral ICA stenoses and occlusions impact ipsilateral ICA WSS and pressures⁴⁹. However, modeling of the entire circle of Willis is resource-intensive and limits the clinical utility of our current protocol.

An additional limitation of our model is that we did not allow changes in resistance and compliance at the ICA outlet and thus do not account for cerebral autoregulation which may affect the distribution of blood flow with different stenosis severities. Furthermore, we modeled the vessel walls as rigid, rather than deformable. However, because carotid artery stenosis is associated with increased vessel stiffness, we feel that a rigid wall assumption is reasonable. In addition, there are not well-defined thresholds of WSS and PG for differing levels of ICA stenosis and associations with stroke risk are not yet defined, thus in our current model we have not validated clinical translatability and cannot yet estimate a patient's stroke risk. Finally, the qMatch MRI sequence is not readily available on standard MRI scanners. qMatch requires a 3T MRI machine and requires manual importing of the sequence as it is not a standard clinical MRI sequence. Moreover, as our protocol specifies qMatch requires complex post-processing in MATLAB which may further limit its generalizability to widespread clinical use.

New metrics to define and assess the hemodynamic impact of ICA stenosis, and better stratify individualized stroke risk are warranted, as evidenced by the current top research priority from the Society of Vascular Surgery: to develop diagnostic tools, imaging techniques, and selection strategies aimed at identifying patients who would benefit from treatment of asymptomatic ICA stenosis⁵⁰. This protocol is well equipped to non-invasively characterize both the hemodynamic loads and plaque composition across an ICA stenosis, thus providing a more comprehensive assessment of ICA plaque embolic potential than current diagnostic modalities. In our future work, we seek to better define the association of hemodynamics metrics (such as WSS and PG) with risk of ICA plaque embolism and stroke.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
CRIMSON	N/A	N/A	Open source online software
Horos	Horos	N/A	Open source online software
MATLAB version 14	Mathworks	N/A
Paraview	N/A	N/A	Open source online software
Siemens 3T VIDA MRI scanner	Siemens Healthineers	N/A