Optical Coherence Tomography Based Biomechanical Fluid-Structure Interaction Analysis of Coronary Atherosclerosis Progression

Summary

There is a need to determine which atherosclerotic lesions will progress in the coronary vasculature to guide intervention before myocardial infarction occurs. This article outlines the biomechanical modeling of arteries from Optical Coherence Tomography using fluid-structure interaction techniques in a commercial finite element solver to help predict this progression.

Abstract

In this paper, we present a complete workflow for the biomechanical analysis of atherosclerotic plaque in the coronary vasculature. With atherosclerosis as one of the leading causes of global death, morbidity and economic burden, novel ways of analyzing and predicting its progression are needed. One such computational method is the use of fluid-structure interaction (FSI) to analyze the interaction between the blood flow and artery/plaque domains. Coupled with in vivo imaging, this approach could be tailored to each patient, assisting in differentiating between stable and unstable plaques. We outline the three-dimensional reconstruction process, making use of intravascular Optical Coherence Tomography (OCT) and invasive coronary angiography (ICA). The extraction of boundary conditions for the simulation, including replicating the three-dimensional motion of the artery, is discussed before the setup and analysis is conducted in a commercial finite element solver. The procedure for describing the highly nonlinear hyperelastic properties of the artery wall and the pulsatile blood velocity/pressure is outlined along with setting up the system coupling between the two domains. We demonstrate the procedure by analyzing a non-culprit, mildly stenotic, lipid-rich plaque in a patient following myocardial infarction. Established and emerging markers related to atherosclerotic plaque progression, such as wall shear stress and local normalized helicity, respectively, are discussed and related to the structural response in the artery wall and plaque. Finally, we translate the results to potential clinical relevance, discuss limitations, and outline areas for further development. The method described in this paper shows promise for aiding in the determination of sites at risk of atherosclerotic progression and, hence, could assist in managing the significant death, morbidity, and economic burden of atherosclerosis.

Introduction

Coronary artery disease (CAD) is the most common type of heart disease and one of the leading causes of death and economic burden globally^1,2. In the United States, roughly one in every eight deaths is attributed to CAD^3,4, while most global deaths from CAD are now seen in low- and middle-income countries⁵. Atherosclerosis is the predominant driver of these deaths, with plaque rupture or erosion leading to coronary artery occlusion and acute myocardial infarction (AMI)⁶. Even after revascularization of culprit coronary lesions, patients have substantial risk of recurrent major adverse cardiovascular events (MACE) after AMI, largely due to the concomitant presence of other non-culprit plaques that are also vulnerable to rupture⁷. Intracoronary imaging provides an opportunity to detect these high-risk plaques⁸. Although intravascular ultrasound (IVUS) is the gold standard for evaluating plaque volume, it has limited resolution to identify microstructural features of vulnerable plaque in contrast to the high resolution (10-20 µm) of optical coherence tomography (OCT). A thin and inflamed fibrous cap overlying a large lipid pool has been demonstrated to be the most important signature of a vulnerable plaque⁹ and is best identified and measured by OCT among currently available intracoronary imaging modalities¹⁰. Importantly, OCT is also able to assess other high-risk plaque features, including: lipid arc; macrophage infiltration; presence of thin cap fibroatheroma (TCFA), which is defined as lipid-rich core with overlying thin fibrous cap (<65 µm); spotty calcification; and plaque microchannels. OCT detection of these high-risk features in non-culprit plaques post-AMI has been associated with up to a 6-fold increased risk of future MACE¹¹. However, despite this, the ability of angiography and OCT imaging to predict which coronary plaques will progress and ultimately rupture or erode is limited, with positive predictive values of only 20%-30%⁸. This limited predictive ability hinders clinical decision-making around which non-culprit plaques to treat (e.g., by stenting)^7,12.

In addition to patient factors and the biological characteristics of plaque, biomechanical forces in the coronary arteries are also important determinants of plaque progression and instability¹³. One technique that shows promise for helping to comprehensively evaluate these forces is fluid-structure interaction (FSI)¹⁴ simulation. Wall shear stress (WSS), also called endothelial shear stress, has been a traditional focal point for coronary biomechanics research¹⁵, with a general understanding that WSS plays an etiological role in atherosclerosis formation¹⁶. Predominantly simulated using computational fluid dynamics (CFD) techniques, low WSS regions have been associated with intimal thickening¹⁷, vascular remodeling¹⁸ and the prediction of lesion progression¹⁹ and future MACE²⁰. Recent advances in these analyses suggest the underlying WSS vector field topology²¹, and its multidirectional characteristics²², as a better predictor of atherosclerosis risk than WSS magnitude alone. However, WSS only captures a glimpse of the overall biomechanical system at the lumen wall, and much like imaging modalities, no one biomechanical metric can reliably discern high risk atherosclerotic features.

Further metrics are emerging as potentially important in atherosclerosis formation. Intraluminal flow characteristics²³ are one such example, with helical flow, quantified through various indices²⁴, suggested as playing an atheroprotective role by suppressing disturbed flow patterns^25,26. While CFD techniques can analyze these flow characteristics and present a wide range of useful results, they do not consider the underlying interactions between the blood flow, artery structure and general heart motion. This simplification of the dynamic system to a rigid wall misses potentially critical results such as fibrous cap stress. While the debate both for and against the need for FSI over CFD continues^27,28,29, many comparisons neglect to include the impact of ventricle function. This limitation can be overcome with FSI, which has shown that dynamic bending and compression exerted on the artery through the influence of the ventricle function can significantly impact plaque and artery structural stress as well as flow metrics such as WSS^30,31,32. This is important as structural stresses are also a key metric for analyzing and predicting plaque rupture^33,34 and have been suggested to co-locate with regions of plaque increase^14,35. Capturing these interactions allows for a more realistic representation of the coronary environment and the potential mechanisms of disease progression.

Addressing this, here we outline the process of developing a patient-specific geometry from OCT imaging³⁶ and the setting up and running of an artery FSI simulation using a commercial finite element solver. The process to manually extract the lumen, lipid and outer artery wall is detailed before the three-dimensional computational reconstruction of the patient's artery. We outline the simulation set-up, coupling and the process of comparing baseline, and follow-up OCT imaging parameters to determine lesion progression. Finally, we discuss the post-processing of numerical results and how these data may have clinical relevance by comparing the biomechanical results with lesion progression/regression. The overall method is demonstrated on non-culprit, mildly stenotic, lipid-rich plaques in the right coronary artery (RCA) of a 58-year-old Caucasian male patient who presented with an acute non-ST elevation myocardial infarction in the setting of hypertension, type 2 diabetes mellitus, obesity (BMI 32.6) and a family history of premature CAD. Coronary angiography and OCT imaging were performed during his initial admission, and then 12 months later as part of an ongoing clinical trial (COCOMO-ACS trial ACTRN12618000809235). We anticipate that this technique can be further refined and used for identifying coronary plaques that are at high risk of progressing.

Protocol

The following deidentified data was analyzed from a patient recruited into the ongoing COCOMO-ACS randomized-controlled trial (ACTRN12618000809235; Royal Adelaide Hospital HREC reference number: HREC/17/RAH/366), with additional ethics approval granted by Central Adelaide Local Health Network (CALHN) Research Services for the purpose of biomechanical simulation (CALHN Reference Number 14179). Figure 1 summarizes the complete workflow outlined in the following protocol, which can be applied to any FSI capable software or codes.

1. Image evaluation

1. Match baseline and OCT follow-up images using anatomical landmarks such as bifurcations and using images immediately proximal to the distal bifurcation and distal to the most proximal bifurcation. The matched images between these landmarks are to be analyzed, as described in Figure 2A.
2. OCT lumen cross-section
   1. Load each OCT image into the image digitizer and click to mark points at the catheter center point and limits of the scale (Figure 2B). Export these points to be used later.
   2. Manually mark the edge of the lumen, starting at the same location in every image, being sure to capture the curves of the lumen as accurately as possible. Leave a gap at the catheter artefact as the reconstruction process will interpolate across this region at a later stage. Export these files in .dat format and repeat this for every image.
3. OCT outer wall and lipids
   1. In the DICOM software, extract the outer wall in high attenuation regions by using visible parts of the outer elastic membrane to manually fit an ellipse to estimate the outer wall location, as described in Figure 3. Click and drag the left mouse button to define the ellipse and position appropriately.
   2. Manually define the lipid arc, calculated to the lumen centroid, and fibrous cap thickness, as described in Figure 3, by clicking and dragging angle and distance measures, respectively. These will be used to analyze lesion progression along with lumen area.
   3. Import these overlaid images into the image digitizer and manually select the outer wall points, using the fitted ellipse as a guide in regions of high attenuation where the outer elastic membrane is not visible. Repeat step 1.2.2 to select and export the points to a .dat format.
   4. Similarly for the lipids, manually select the lipid surface, starting from the same end of the lipid in every case. Use the outer wall ellipsoidal guide (step 1.3.1) for a consistent backside arc. Export points to a .dat file and repeat for all images with lipids present, leaving a gap across the guidewire artefact as described in step 1.2.2.
      NOTE: Lesion progression is analyzed by comparing three metrics, namely, lumen area, lipid arc and fibrous cap thickness, which can be assessed directly from the DICOM viewer. The technique to extract the outer wall and lipid backside is required due to OCT's limited penetration depth. OCT was used in this investigation due to the focus on the relationship between plaque composition and biomechanical forces.
4. Angiography based centerline
   1. Load the first angiographic image in the image digitizer³⁷. Select the edges of the catheter to scale the image in later steps, and then manually mark the catheter centerline beginning with the proximal marker and moving distally, with evenly spaced points, as shown in Figure 4A. Export the data to .dat format and repeat for the second angiographic plane.
      NOTE: Generally, planes with an angle greater that 20° between them improves the three-dimensional centerline reconstruction robustness. The catheter and OCT guidewire should be visible in each image.

2. Three-dimensional reconstruction

1. Angiography projections
   1. Load the data files that were exported in step 1.4. Use the first two points to scale the data to millimeters (the first two points are used with the known catheter specifications, 6F in this case). Subtract the proximal data point from the remaining points in each dataset so that the curve begins at the origin of the coordinate system.
   2. Generate rotation matrices for each angiographic view, where θ and Φ represent the RAO/LAO and CAU/CRA angles, respectively. We use LAO and CRA angles as negative. The two rotation matrices in the x (Rot_x) and y (Rot_y) directions, respectively, are:
           (1)
   3. Multiply the rotation matrices together, and then multiply them with the coordinates of each point from step 2.1.1. The resultant equation:
            (2)
      gives the three-dimensional location of the catheter point on its respective angiogram plane (Pt3D) by rotating the two-dimensional points that were specified from each angiographic image.
   4. Calculate the normal vector to each angiographic plane by multiplying the x and y rotation matrices by the unit vector in the z direction. From proximal to distal location, project each point normal to its respective plane and calculate the midpoint of the shortest distance between the projections. This results in the three-dimensional point on the OCT guidewire in space.
   5. Using the 'interparc' function, available from the MATLAB central file exchange³⁸, divide the three-dimensional centerline into equally spaced points. The spacing between points should be equal to the spacing between OCT images, which is determined by the pullback speed. These are the locations where the OCT cross-sections will be placed.
2. OCT cross-section rotation
   1. Using the data file containing the catheter center and scale, convert each cross-section from pixels to mm using the second and third points in the scaling file. To center the cross-section about the catheter location, subtract the first point in the scaling file (the catheter center) from all cross-section points. Calculate the normal vector to the cross-section (parallel to the catheter in the artery) by subtracting the three-dimensional centerline point from the next distal point along the catheter curve.
   2. Rotate the OCT cross-section to align perpendicular to the catheter centerline by multiplying the scaled data points by the rotation matrix:
           (3)
      where
           (4)
      and N_X, N_Y, and N_Z are the x, y, and z components, respectively, of the normal vector calculated in section 2.1. Add the three-dimensional centerline point to all the rotated points in the cross-section, resulting in the cross-section location in three-dimensional space (Figure 4B).
   3. Repeat steps 2.2.1-2.2.2 for every cross-section (lumen, artery, and lipid). Export the cross-sections to a text file, which can be imported into the computer aided design (CAD) software for final solid body creation.
3. 3D solid model creation
   1. In a 3D modeling software, import and generate the cross sections one file at a time. Import the text files containing the cross-sections into the 3D modeling software by clicking on the concept drop-down box (Figure 5A-1) and selecting 3D curve (Figure 5A-2). Click on Generate.
   2. To create a solid component, select all the curves in order and loft them together (Figure 5A-3), ensuring add frozen is selected to generate a new solid. Carry out these steps for the lumen, lipids, and outer wall to create separate solids, ensuring to enable merge topology.
      NOTE: It may be necessary to skip a curve if problematic geometry arises. In this reconstruction, omit a small mid-section lipid due to its size and the added computational cost and numerical complexity associated with its inclusion.
   3. To subtract the lumen and lipids from the artery wall, create a Boolean operation from the create drop-down list and choose the target body as the wall and the lipids/lumen as the tool bodies to subtract the lumen and lipids from the artery wall (Figure 5A-4).
   4. Share topology between the wall and lipids to ensure mesh nodes are shared in future steps. To do this, manually highlight the wall and lipids and right click to form a new part (Figure 5A-5).
      NOTE: This step ensures mesh nodes are shared between the surfaces preventing improper contact regions or mesh penetration between the two layers, greatly assisting in the solution phase. The final geometry of the catheter centerline, lipids, lumen and artery wall is visualized in Figure 5B.
4. Pre-processing: Boundary conditions
   NOTE: Before setting up the simulation, patient specific Boundary Conditions (BC's) are needed. Here displacement extracted from the angiography was used, which is applied to the inlet and outlet of the simulation and blood flow velocity/pressure measured from human patients and described in the literature³⁹.
   1. Displacement
      1. Repeat steps 1.4 and 2.1, but choosing only the distal and proximal markers, beginning with the angiographic image immediately preceding end-diastole. Do this for all angiographic images over one cardiac cycle.
      2. Fit smoothing splines to the x, y, and z coordinates of the two sets of points. This results in the displacement of the inlet and outlet regions. Representative results for the patient displacements are shown in Figure 6A.
         NOTE: Displacement analysis was begun at the image preceding end-diastole to best match phases between the extracted displacement and the applied pressure and velocity profiles found in section 3.1.2, whose systolic phase begins at 0.1 s (corresponding to the spacing between angiographic images). When extracting motion, ensure that there is no table panning/image movement throughout the image set.
   2. Blood velocity/pressure
      1. Create profiles that describe the pulsatile blood velocity and pressure by compiling User Defined Functions (UDF). Here transient profiles measured from human patients in the literature were applied ³⁹, modeled as a Fourier series, mathematically described by:
         ,     (5)
         where t is the time, w₀ is the frequency, T is the signal period, n is the number of terms, and a_0-11, b_1-11 are coefficients fitted to profiles described in the literature. In this case, we are using the first 11 terms.
      2. NOTE: These profiles are described in Figure 6B and should be written to in a C formatted file in an integrated development environment such as Microsoft Visual Studio. Outlet pressure is a flat profile and inlet velocity is applied as a fully developed, parabolic profile, described as sufficient to reproduce realistic conditions⁴⁰. Further development of this procedure could include measuring patient blood velocity (such as by doppler echocardiography⁴¹) and pressure (using pressure wires) as more realistic boundary conditions. Furthermore, concurrently measuring displacement, blood velocity and pressure would ensure their phases are accurately matched.

3. Artery/structural

1. To set the material properties for the artery and lipid, enter engineering data and add a new material called artery. Drag density and the 5-parameter Mooney-Rivlin model to the new material and set their parameters. Enter a density of 1,000 kg/m³ and the hyperelastic coefficients described in Table 1, based on intima⁴² and lipid⁴³ properties in the literature. Repeat this for the lipid.
   NOTE: The Mooney-Rivlin model is described by⁴⁴:
        (6)
   Where c₁₀, c₀₁, c₂₀, c₁₁, and c₀₂ are material constants and d is the incompressibility parameter (zero for incompressible material in this case). Here I_x is the x^th invariant of the strain tensor and J is the elastic deformation gradient determinant.
2. Enter the model component, suppress the lumen/fluid component by right clicking on the Lumen/Fluid and selecting Suppress (Figure 7A). Assign the previously defined materials to the artery and lipid solids by selecting them from the material drop down list, checking whether the units are appropriate.
3. The geometry now needs to be meshed. Click on mesh (Figure 7B), set the physics preference to nonlinear mechanical and specify the mesh sizing. Here adaptive meshing with a target size of 0.14 mm was used. Adjust the mesh preferences as needed to obtain reasonable mesh skewness values and aim for at least two to three mesh elements across gaps such as the fibrous cap. Generating the mesh may take some time due to the complex geometry.
   NOTE: A mesh independence study must be conducted to ensure results are not impacted by mesh characteristics. Gradually decrease mesh size and compare results until variation is less than a set limit; in this case, we use 2%⁴⁵ (measured at the fibrous cap of the third plaque). Furthermore, to ensure mesh quality, check mesh skewness; high mesh skewness will result in numerical difficulties during convergence or inaccurate results. To lower skewness, try decreasing the mesh size or adjust the growth rate, maximum size and/or curvature angle. Results for our mesh independence test are outlined in Table 2, with percentage variation in results compared to the medium mesh sizing, which was used throughout this analysis.
4. Click on Analysis settings (Figure 7C). For FSI simulations, turn automatic time stepping off and set the number of substeps to one (system coupling will control substeps), set the simulation end time, in this case 0.8 s (patient heart rate of 75 bpm). System coupling will control the time and substeps.
5. In the solver controls drop-down list, set the solver type to program controlled to use either the direct or iterative method. Direct methods are more robust but use a significant amount more memory. Set the Newton-Raphson method to full. (Due to the complexity of the geometry and nonlinearity in the simulation, the direct method and the full Newton-Raphson iterative method may be needed; however, these significantly increase the computation cost.)
6. Specify the system coupling domain as the inner wall of the artery by inserting a fluid-solid interface. Do this by right clicking and inset a fluid-solid interface under the Transient tab (Figure 7D). Select the inside of the artery wall for the interface. This will pass data between the structure and fluid at this location.
7. The displacement boundary conditions can be entered as a displacement function in the x, y, and z direction applied at the inlet and outlets. Do this by right clicking under the Transient tab and inserting displacements (Figure 7E). Duplicate the displacement for the x, y, and z directions. In the direction drop-down list, select function and copy the displacements extracted in step 2.4.1.
   NOTE: Displacement can be input as a function or as a table of points depending on preferences.
8. To assist in troubleshooting errors, under the Solution tab, insert four Newton-Raphson residuals. These can be viewed if errors arise to find the troublesome geometry or mesh locations.
   NOTE: To insert post-processing options, such as maximum principal stress, right click on the Solution tab and insert the appropriate results (Figure 7F).

4. Blood/fluid

1. Enter the Model tab, check the units, and suppress the artery and lipid part, leaving the fluid domain, in similar fashion to step 3.2.
2. Specify the mesh metrics and generate the mesh, checking skewness and adjusting if necessary (we applied a 0.14 mm mesh size with a maximum wall size of 0.12 mm). It is good practice to use similar mesh size and shape, as done in the structural part, on the areas where the fluid-solid interaction is occurring.
   NOTE: As with step 3.3, a mesh independence test should be conducted to ensure results are independent of mesh properties, as shown in Table 2. Check the quality of the mesh and adjust the element size, growth rate, refinement or curvature as needed to ensure skewness remains low and that mesh independence is reached.
3. Create named selections for the inlet, outlet, and wall before entering the fluid setup, by right clicking on the respective surface and selecting insert named selection.
4. Enter the Setup tab and ensure double precision is enabled. Set the Solver type to Pressure Based and ensure that the Time is set to Transient by checking their respective tick boxes (Figure 8A).
5. Enable the k-omega Viscous Turbulence model and enable Shear Stress Transport and Low Reynolds Corrections by entering the Viscous Models tab (Figure 8B) and checking their respective tick boxes.
6. To enable nonlinear viscosity models with turbulence, enter the command '/define/models/viscous/turbulence-expert/turb-non-newtonian?' in the command console (Figure 8C) and enter 'yes' when prompted.
7. Under Materials (Figure 8D), define the blood properties by entering density and selecting non-Newtonian power law in the viscosity drop down list. Do this by renaming the fluid as blood, set a density of 1,050 kg/m³, and set the Power-Law non-Newtonian consistency index, k, to 0.035, the power law index, n, to 0.6.
   NOTE: The Power Law non-Newtonian viscosity model was chosen based on literature to describe the nonlinear blood viscosity⁴⁶, η, in terms of the fluid strain rate, , as:
        (7)
   Various non-Newtonian blood viscosity models exist to capture the shear-thinning nature of blood. Several publications^46,47,48,49 have investigated the efficacy of various viscosity models and their coefficients, which should be consulted for further information when choosing the appropriate model.
8. Compile our user defined function, previously described in step 2.4.2, containing the transient blood velocity and pressure, checking the command lines for any errors (Figure 8C). Now load the UDF by entering the User Defined tab (Figure 8E), selecting Compiled and navigating to the directory of the UDF before importing it and clicking on Build, and then on Load.
   NOTE: Text will appear in the console (Figure 8C). Check this carefully to ensure no errors or warnings appear. If the UDF loads correctly, the UDF names will appear in the console, (highlighted in Figure 8C).
9. These can be applied to the inlet and outlet. To do this, select the Boundary Conditions tab. Double click on Inlet (Figure 8F) and choose the inlet UDF from the profile drop-down list. Repeat this step to also define the outlet pressure.
10. Enable the dynamic mesh (by checking the tick box under the Dynamic Mesh tab shown in Figure 8G), including smoothing, remeshing, and 6° of freedom solver tick boxes, setting the diffusion parameter to 1.5 and the appropriate maximum and minimum scales for your mesh.
11. Ensure that the maximum and minimum mesh scales are within the limits of the mesh zone and that the target skewness is set to 0.7. The mesh properties can be shown by clicking on the Mesh Properties tab.
12. Create a new dynamic mesh zone by clicking on the Create button, specify the wall of the lumen in the Region drop down list and select System Coupling. This is the interface to pass data to the artery component of the simulation.
13. Create deforming mesh zones for the inlet, outlet, and interior lumen with appropriate values for the mesh scale. Do this by clicking on Create in the Dynamic Mesh tab and choosing Deforming. Enable remeshing and smoothing and set the mesh scales based on the limits of each zone. Often, negative cell volume errors are associated with this dynamic mesh, so check carefully and adjust the mesh scales if needed for each region.
14. Ensure the pressure-velocity coupling is set to coupled and set the transient formulation and spatial discretization schemes to second order by entering the Methods tab (Figure 8H) and making selections from the respective drop-down lists.
15. In controls (Figure 8H), enter a courant number of two, and set the residual convergence criteria in the Monitors tab (Figure 8I). We have used a value of 1e-5 for continuity and 1e-6 for the remainder.
    NOTE: The Courant number can be estimated based on the mesh size, dx, time step size, dt, and velocity of the blood, v, using:
         (8)
    Enter this number under the courant number section in the Controls tab (Figure 8H). Here we apply a Courant number of two. The Courant number is generally less than one; however, as coupled pressure-velocity solver with implicit solution methods is used, the result is inherently more stable and less sensitive to this value; hence, two is considered as acceptable.
16. To define a custom function for results such as local normalized helicity (LNH), select custom field functions under the Parameters and Customization tab (Figure 8J) and insert a new function by right clicking and selecting New. Use the popup window to define as necessary. Enter the formula using the drop-down list of solver variables. As a representative result, we use LNH^50,51, a measure of the alignment between the velocity, , and vorticity, ω, vectors, as a custom function described by:
          (9)
    NOTE: Other custom variables should be defined at this step, such as the oscillatory shear index (OSI)^52,53, a measure of flow reversal.
17. In the Run Calculation tab (Figure 8K), set the number of timesteps to 160 (a step size of 0.005 s and end time of 0.8 s), time step size of 5 ms and the number of iterations to 300 to ensure the result is time independent.
    NOTE: Depending on the complexity of the simulation, greater iterations per step may be required. Multiple cardiac cycles may be required for complete numerical convergence, something we note as a limitation; however, this is often applied in coronary biomechanics simulations due to the computational cost associated with these simulations.
18. Check whether the Data Sampling for Time Statistics tick box is enabled and ensure that Wall Statistics and Flow Shear Stresses are selected, as well as the previously defined custom function.
19. Create the data export in the Calculation Activities and Autosave tab (Figure 8L), selecting the CFD-Post Compatible option for post processing. If one wishes to process results in a separate software, adjust the export type as necessary. Select all regions (wall, interior mesh, inlet, outlet) and the results to be exported.
20. Finally, initialize the simulation with the hybrid scheme by entering the Initialization tab (Figure 8M), selecting the Hybrid scheme, clicking on Settings, and increasing the number of iterations to 20. Click on Initialize.

5. System coupling

1. Make sure both the structural and fluid setups are connected to system coupling and updated. Do this by clicking and dragging the structural and fluid setup to system coupling to link them, as shown in Figure 9A, ensuring both setups are updated by right clicking and selecting Update.
2. In System Coupling, set the end time to 0.8 s and the timestep to 0.005 s. Do this by selecting Analysis Settings (Figure 9B-1) and entering the end time and time step size. Set the maximum iterations to 10.
   NOTE: Generally, between 10 and 15 iterations is sufficient if both the structural and fluid components are converging well.
3. Select the wall and solid interface from the fluid and structural components, respectively, and add a data transfer by holding Ctrl and selecting the two fluid-structure interfaces (Figure 9B-2); right click and create a data transfer between the fluid and structural components (Figure 9B-3). Adjust the under-relaxation or ramping of the force being transferred from fluid to structure to assist in convergence.
   NOTE: Depending on the complexity of the model, boundary conditions and material properties, data transfer ramping or under-relaxation may be required for numerical convergence. These can be applied to the fluid data transfer (i.e., the force being transferred from the fluid component to the artery wall). These options are available within the created data transfers (Figure 9B-2).
4. When ready to run, click on Update. Simulation data such as structural and fluid convergence and their respective data transfer convergence are printed in the console.
   NOTE: Note that FSI simulations are computationally expensive, with this simulation taking 11 days on a 16-core machine (2.6 GHz Intel Xeon Gold using 180 Gb of physical memory (RAM)), with further variation in simulation times depending on hardware setup and model complexity. Representative data transfer residuals are shown in the chart (Figure 9B-4) and solution data is printed in the console (Figure 9B-5). Over the first few iterations, convergence of the data transfer residuals may not be completely obtained until an equilibrium state is reached. This is described in more detail in the caption for Figure 9B.
5. When the simulation is complete, the results can be postprocessed within the commercial software or in a separate software, depending on your data export type described in step 4.19.

Representative Results

Representative results are presented for both established and emerging biomechanical markers of atherosclerosis progression. Established metrics such as WSS and WSS-derived results (including time averaged wall shear stress (TAWSS) and oscillatory shear index (OSI)) are visualized in Figure 10. The wall shear stress over the cardiac cycle is largely driven by the blood velocity, however, artery geometry and its motion/contraction play a significant role in its spatial distribution. This can be seen in the TAWSS and OSI contours, with OSI, a measure of flow recirculation, having been linked to atherosclerosis formation. Such variations are representative of narrowings/bends and the complex nature of artery geometry. Results near the inlet and outlet should be disregarded as they will be significantly impacted by boundary conditions. While WSS as a scalar value has been well studied, the underlying vector field, which presents both a scalar value and direction has received far less attention. Here the wall shear stress vector field is shown in Figure 11, with areas highlighted to show attraction and expansion regions, which could influence near wall transport processes. The distal region in particular shows a significant attraction region, representative of sudden changes in lumen area, which could suggest a risk site for atherogenesis. The attraction and expansion regions have been visually assessed here; however, this result could be extended through further post-processing to extract the underlying topological structure and fixed points⁵⁴ numerically. Such attraction/expansion regions are representative of changes in lumen area and bifurcation regions.

Intraluminal flow characteristics (away from the lumen wall) could also play an important role in atherosclerosis progression. Local normalized helicity (LNH), a measure of alignment between blood velocity and vorticity vectors, is useful for visualizing counter-rotating intraluminal flow characteristics²⁶ and is shown in Figure 12. A distinct difference can be seen between the proximal and distal artery sections. If the developed simulation was sub-optimal (i.e., boundary conditions were very close to the region being analyzed or proximal geometry such as bends are omitted to focus on a specific lesion), this result could be significantly affected, hence, impacting the association with changes in plaque morphology. This is highlighted inset in Figure 12. Finally, the stresses within the artery wall are analyzed using the Von Mises effective stress and are presented in Figure 13. Note the high stress region at the outlet (distal) due to the boundary condition, which should be ignored. Von Mises stress (VMS) is impacted by a combination of fibrous cap thickness, artery geometry, material properties, blood velocity/pressure and artery motion/contraction. Due to this, biomechanical simulation is required to determine the stresses in the wall that result from the nonlinear interaction between these factors. The thin fibrous cap in the proximal region results in the largest stress, with the transient distribution over the cardiac cycle driven by artery motion. This is representative of the thin nature of the fibrous cap. Conversely, in the distal region, the VMS resulting at the fibrous cap is predominantly driven by blood pressure. Capturing these results is only possible through fully coupled FSI simulations due to the inherent nonlinearity in the simulation. These variations could play a yet undefined role in atherosclerosis progression.

Finally, lumen area, total lipid arc and minimum fibrous cap thickness were compared between baseline and follow-up OCT imaging. These results are shown in Figure 14, with OCT images of selected regions shown for visualization of the changes. Results are generalized at the highlighted regions to understand the role WSS, intraluminal flow, and structural stresses have on lesion progression; however, further post-processing could be carried out to quantitatively compare biomechanical results at each OCT cross-section location with changes in plaque composition. The progression of the mid-artery lesion was associated with initially low TAWSS and VMS, disturbed LNH patterns, and saw a strong WSS attraction region (WSSat) and OSI. Note that this lesion was not included in the initial geometry due to its small size and the added computational expense and complexity associated with its inclusion. This is directly contrasted with the proximal region, which saw more consistent LNH patterns, high TAWSS and VMS and weaker WSS expansion region and OSI associated with lipid arc reduction and a thinner fibrous cap. The distal lesion saw the most significant progression and was associated with a moderate VMS that was driven by blood pressure, not artery motion, unlike the upstream regions. When taken together, this methodology and results could be used to develop a more complete picture of a patient's potential risk, provided a larger number of simulations are able to be produced to better understand the statistical significance of results.

	C10 (MPa)	C01 (MPa)	C20 (MPa)	C11 (MPa)	C02 (MPa)	d (Pa-1)
Artery	-0.19	2.03	11.3	-0.19	20.1	1.00E-05
Lipid	-0.17	0.21	5.02	-1.88	13.5	1.00E-05

Table 1: Material properties for the artery wall and lipid using a five-parameter Mooney-Rivlin model.

Structural
	Average mesh size (mm)	Target Elements Over Gap	Average Mesh Skewness	Number of Elements	Maximum Principal Stress Variation (%)
Coarse	0.17	2	0.25	1,266,029	4.7
Medium	0.14	3	0.25	1,657,589	-
Fine	0.11	4	0.24	3,382,733	1.28
Fluid
	Average Mesh Size (mm)	Max Face Size (mm)	Average Mesh Skewness	Number of Elements	Wall Shear Stress Variation (%)
Coarse	0.17	0.15	0.23	527,103	6.42
Medium	0.14	0.12	0.22	772,767	-
Fine	0.11	0.09	0.22	1,392,534	1.85

Table 2: Mesh independence results for three different mesh sizes tested.

Figure 1: Workflow for conducting fluid-structure interaction analyses from patient imaging. Please click here to view a larger version of this figure.

Figure 2: Outline of the process for aligning OCT images between baseline and follow-up before extracting the lumen edge. (A) Schematic of OCT cross-section alignment based on landmarks such as the proximal bifurcation (*) and distal side branch (^) to select corresponding images from the same artery segment at both baseline and follow-up. (B) Visualization of the process for selecting points for the catheter center and scale (red crosses) as well as the lumen start point based on the cross-section line (white dashed) and the lumen outline (white dots). Please click here to view a larger version of this figure.

Figure 3: Outline of the process for reconstructing the backside of lipids and the outer edge of the artery wall in regions of high attenuation. Manual selection of visible points on the outer elastic membrane (red crosses); projection of an ellipse fitted to the selected points to form the outer artery cross-section estimate (white dashed line); manual selection of inner lipid surface (blue dots); projection of the lipid backside (yellow dots) by shrinking the outer wall estimate to reach the desired backside thickness. Overview of measurements used for lesion comparisons are also noted (lumen area, fibrous cap thickness, and lipid arc). Please click here to view a larger version of this figure.

Figure 4: Outline of biplane angiogram catheter point selection and three-dimensional reconstruction. (A) Angiogram with scaling points and catheter curve highlighted. (B) OCT based cross-sections of the lumen (blue) rotated and placed along the three-dimensional catheter centerline (red points). Please click here to view a larger version of this figure.

Figure 5: Software layout for generation of three-dimensional geometry. (A) 3D solid model creation: (1) Drop-down menus for operations; (2) Imported 3D cross-sections extracted from OCT; (3) Loft between cross-sections to create solid; (4) Boolean operation for uniting/subtracting solid bodies; (5) Shared topology by creating a part with artery and lipids. (B) The final 3D solid geometry components, including the three lipids, catheter centerline, lumen, and artery wall. Note that a small mid-section lipid was not considered in the reconstruction due to its small size and the added computational cost associated with its inclusion. Please click here to view a larger version of this figure.

Figure 6: Boundary conditions applied to the simulation. (A) Plot of displacement in the x, y and z directions for the proximal (inlet) and distal (outlet) sections of the artery extracted from angiography. (B) Velocity and pressure boundary conditions for the inlet and outlet, respectively, from the literature. Please click here to view a larger version of this figure.

Figure 7: Software layout for the structural setup. (A) Solid bodies and their material property definitions; (B) Mesh settings; (C) Solver/analysis settings; (D) Fluid-solid interface applied to the inner wall of the artery; (E) Displacement boundary conditions applied to the inlet and outlet; (F) Solution post-processing. Please click here to view a larger version of this figure.

Figure 8: Software layout for the fluid setup. (A) General solver settings; (B) Turbulence model settings; (C) Console/text user interface for typed commands; (D) Fluid property setup; (E) User-defined function compiler; (F) Boundary condition settings for specified named selections; (G) Dynamic mesh settings and fluid-structure interaction zone setup; (H) Solver settings and controls; (I) Specification of residual convergence criteria; (J) Region to define custom functions such as LNH; (K) Calculation setup for timesteps and iterations; (L) Data export format settings; (M) Solution initialization. Please click here to view a larger version of this figure.

Figure 9: Outline of system coupling interface. (A) Final layout of the coupled fluid-structure interaction workflow. (B) Software layout for system coupling setup and simulation solve: (1) Analysis settings; (2) Fluid-structure interfaces for the structural and fluid components; (3) Data transfers; (4) System coupling residuals monitor; (5) Solution information. While convergence of the structural and fluid components may be obtained in every step, the data transfer residuals may not converge together over the first few time steps. By the seventh step, the data transfer residuals are tracking within 10 percent of each other, showing good convergence. For complete convergence, data transfers must also reach the target criteria, with the default set to minus two (logarithmic scale), shown in the green arrow. The start of each new time step is outlined by the blue arrows. In this case, increasing the number of system coupling iterations could allow the data transfers to converge earlier but adds computational cost to the simulation. Simulation convergence was reached from step 11 onwards. Please click here to view a larger version of this figure.

Figure 10: Representative wall shear stress results. (A) Time averaged wall shear stress, (B) oscillatory shear index and (C) wall shear stress over the cardiac cycle at the four locations defined on contour (A). The wall shear stress result variations are driven by a combination of the lumen shape, its motion/bending, and blood velocity, resulting in a peak in shear stress over the distal lipid (point four). The fluctuations seen over the initial time steps are a result of the simulation finding equilibrium and should not be considered further when interpreting results. Please click here to view a larger version of this figure.

Figure 11: Wall shear stress vector field colored by magnitude (inverse scale). Analyses of the underlying vector field shows promise for better understanding near wall transport processes, with attraction and expansion regions (highlighted inset) suggesting atheroprone and atheroprotective near wall flow conditions, respectively. Attraction regions are seen more significantly in the distal regions suggesting potential risk for lesion progression. Please click here to view a larger version of this figure.

Figure 12: Cycle averaged local normalized helicity values for right-hand (positive) and left-hand (negative) rotating flow structures. The distal locations see larger disturbances in the counter-rotating flow structures when compared to the proximal region. A suboptimal simulation (shown inset) focused on the distal plaque with too little of the proximal geometry included resulted in varied helical flow structures which alter the associations seen between biomechanical results and plaque progression. Please click here to view a larger version of this figure.

Figure 13: Representative von Mises stress distribution results. (A) Von Mises stress distribution in the artery wall. The regions of high stress are located at the minimum fibrous cap thickness sites. Plaque shoulders produce stress intensifiers. (B) Plot of Von Mises stress over time for one cardiac cycle, with the transient distribution impacted by a combination of plaque morphology, blood pressure, blood velocity/momentum, and the artery motion (cyclic bending and compression). Numbered points in (B) correspond to those shown in (A). Note the difference between points one and four; both result on a fibrous cap; however, stress at the first location is driven by artery motion/bending, while the fourth point is dominated by blood pressure. Please click here to view a larger version of this figure.

Figure 14: Representative change in measured plaque and vessel characteristics compared to biomechanical results. (A) Lumen area, (B) total sum of lipid arc angle and (C) minimum fibrous cap thickness comparisons between baseline and follow-up imaging. Baseline and follow-up images are shown corresponding to the axial distance marked by the green arrows. General behavior for time averaged wall shear stress (TAWSS), oscillatory shear index (OSI), wall shear stress attraction region (WSSat), local normalized helicity (LNH) and the Von Mises stress (VMS) is shown in the highlighted regions. Markers represent increased, decreased, and neutral values for the red arrow, green arrow, and black line, respectively. Please click here to view a larger version of this figure.

Discussion

The use of FSI methods to analyze coronary biomechanics is still a developing field from both numerical modelling and clinical result aspects. Here we have described the outline of setting up a patient specific FSI analysis, based on the finite element/finite volume methods, utilizing OCT and angiographic imaging. While the method we describe here utilizes a commercial finite element solver, the procedure can be applied to any FSI capable software. There are still several limitations to be improved upon in the methodology. Firstly, we acknowledge the limitation of only presenting representative results for a single patient; however, we present the current protocol to better outline the process of biomechanical FSI simulation of arteries in the hope the presented methodology can be further developed and applied to larger cohorts in the future. Furthermore, during the image analysis stage, much of the process described is manual, introducing potential interobserver variability when segmenting the lumen, lipids, and outer wall. Further developments can be made to automate these processes⁵⁵ or include multiple expert analyses which will be particularly useful for providing a more robust comparison between baseline and follow-up images. The OCT imaging used in this procedure was also not ECG-gated, resulting in slight misalignment of cross-sections due to cardiac motion/contraction. Taking this into account in the future could also increase comparison accuracy.

As the focus of the study was on the relationship between biomechanical forces and plaque composition, an OCT centered approach was taken to make use of OCT's high accuracy (10-20 µm). This, however, presents challenges due to the significant attenuation that results from OCT's limited tissue penetration depth⁸. Unlike purely CFD techniques, which only require detailing the lumen shape from OCT images for fluid analyses, FSI also requires detailed information on lipids and the artery wall. To overcome this limitation an estimation technique was presented in this methodology, as during normal clinical practice only one intravascular imaging modality is used. Future research should be undertaken to validate the accuracy of this technique and its impact on biomechanical results through co-registration of IVUS and OCT. For this procedure, we have assumed that this impact is negligible, as lumen shape and fibrous cap thickness^56,57 are generally more dominant drivers of stress results in the wall. Combining OCT and Intravascular Ultrasound (IVUS) could overcome this difficulty⁵⁸; the superior tissue penetration of IVUS allowing analysis of plaque characteristics and a near 10-fold increase in spatial resolution of OCT leading to highly accurate lumen representation⁵⁹. Developments in multimodal catheters present a significant opportunity to address this in the future⁵⁸. Similarly, for the angiographic centerlines, further automation to more accurately account for table panning, artery compression/relaxation, and foreshortening can be made by building on this methodology⁶⁰.

The numerical simulation can be further improved by including the anisotropic nature of coronary tissue. The isotropic hyperelastic tissue described in the methodology does not consider the impact of artery microstructure. This has been shown to be important for stress magnitude and distribution in the artery wall³¹. Tissue anisotropy could be included in the future by creating user-defined material models. Alternate finite element solvers also have built-in anisotropic hyperelastic models developed specifically for biological tissue behavior, such as the Holzapfel-Gasser-Odgen model⁶¹. Current simulations are also limited by primarily using data provided from the literature⁴² excised from cadavers and mechanically tested outside of the physiological environment. As intravascular imaging techniques improve, we also see the possibility of defining tissue properties tailored to each specific patient. Similarly, the intima, media and adventitia were simplified to a single layered structure as separating these layers from OCT presents challenges, especially when ensuring no overlap occurs in the reconstruction process. The added numerical cost is also significant, requiring further investigation to find the balance between cost/time and accuracy as the multiple layers do play a role in the overall stress response^62,63.

Furthermore, from an intravascular flow perspective, while we only present representative results of LNH, several helical flow indices have been developed in the literature, with results such as the h2 index, which quantifies helicity intensity, suggested as being atheroprotective. These indices can be added to the protocol by defining them in custom field functions (step 4.16) in the same manner as LNH and we direct interested readers to the cited literature^24,25,51 for further information.

Finally, our methodology is limited to applying velocity and pressure boundary conditions measured in the literature, rather than measured patient-specific ones. This could be improved by extracting blood velocity and pressure boundary conditions specific to each patient using doppler echocardiography⁶⁴ and quantitative/fractional flow reserve measurements⁶⁵ in the future. Such improvements would be particularly useful in ensuring that artery displacement could be exactly matched to blood velocity and pressure phases (start/end systole in particular), something we provide an estimate for in this analysis. Furthermore, while it has been suggested that fluid boundary conditions only impact the very near region of the boundary, the region of interest for assessment of plaque growth/regression of remodeling should be set around the mid-section of the simulation geometry to minimize any impact of these inlet and outlet conditions. The impact of simulating a sub-optimal artery section that is too short has been highlighted in Figure 12, where the inlet and outlet conditions greatly influence the presented results, rather than the artery geometry and dynamics. This will have direct consequences for predictive capability if not adequately accounted for. Furthermore, while we present a simulation using the Power Law non-Newtonian viscosity model based on suggestions in the literature^46,48, more recent investigations⁴⁷ have suggested the Carreau model as a more appropriate viscosity model. Different viscosity models can be chosen in step 4.7 of the protocol. By building directly from this methodology, incorporation of these advances could greatly increase result specificity and accuracy in the future.

During the development of the solid model, it is important to take steps to improve the likelihood of a high-quality mesh and reduce error prone regions. Sharing topology between the lipid and artery wall is an important first step, enabling mesh nodes to be shared between the solid surfaces when the meshing step is carried out. By sharing nodes, the possibility of errors associated with imperfect mesh alignment and mesh penetration is reduced, a not insignificant risk due to the complex geometry of the model. Following the setup of both the solid and fluid analyses components, it is also critical to check for mesh independence. This ensures that all results are not influenced by the size of the generated mesh. For mesh independence, results should remain stable despite changes in mesh size. For efficient calculations, choose the largest possible mesh (least number of elements) that ensures this independence. Furthermore, ensuring that both the Fluent based residuals converge and the system coupling data transfers converge is a critical step. Generally, it is good practice for the structural and fluid data transfer residuals to converge to within 10% of each other at the end of each time step. As the simulation finds an equilibrium state, the first few iterations of any FSI simulation will sometimes not completely converge, as previously described in Figure 7. Generally, after the first 5%-10% of simulation time is complete, convergence should be reached for every step. While we also demonstrate the methodology here over one cardiac cycle, a simplification frequently made in coronary biomechanics, multiple cardiac cycles are often required for numerical convergence. Due to the computational cost associated with coronary biomechanics simulations, however, multiple cardiac cycles are often not feasible. This is a limitation that warrants further development.

Several possible errors may also arise during simulation. Most common of these are element distortion/deformation from the structural side and negative cell volume from the fluid side. Element distortion can result from stress intensifiers such as pinch points/poor element quality areas or from the applied force exceeding material strength. To find the locations where errors are occurring, insert Newton-Raphson residuals (from under the Solution tab in transient structural). Newton-Raphson residuals will show the regions of the geometry with the largest residual error. If mesh quality is low in this region, refining the mesh may help. Viewing the force/displacement convergence from the drop-down list in the Solution tab may also help troubleshoot this error. From the fluid perspective, negative cell volume is generally associated with the dynamic mesh settings. Review these settings carefully and try increasing the margin between the mesh zone maximum/minimum sizes and the manually specified maxima/minima. Shallow/poor quality elements in the boundary region may also impact this error due to the deformation occurring during system coupling. Reviewing mesh quality in the boundary region may also assist in troubleshooting. When troubleshooting, solving the structural and fluid components separately can reduce the time spent addressing errors.

Moving forward, FSI-based simulation shows significant potential for coronary lesion analyses as it overcomes the limitations of purely structural or fluid (CFD) based approaches. The presented methodology also has further applications in coronary stenting⁶⁶, ventricle function⁶⁷, and spontaneous coronary artery dissection⁶⁸ analysis. However, there remains the challenge of correlating biomechanical results to clinical outcomes. This remains difficult due to the computation cost associated with FSI techniques, limiting the number of analyses able to be performed within realistic timeframes. This means that currently results cannot be used in near-real time (such as in the catheterization laboratory) to directly assist percutaneous coronary intervention (PCI), but rather are retrospective analyses leading to delays in data access. This timeframe may still be helpful in guiding staged PCI procedures or individualizing more aggressive pharmacotherapy or surveillance for non-culprit coronary lesions that are identified as being at high risk of progression or destabilizing due to their biomechanical profile. With multiple biomechanical markers emerging with potential links to CAD, a multifactorial approach capable of simulating and comparing results will give a clearer picture of the relationship between artery biomechanics and atheroprone sites, a possibility that FSI simulations are uniquely placed to complete. The possibility for this FSI simulation to further integrate and inform machine learning algorithms could also overcome the limitations associated with each separate approach⁶⁹. Overall, this methodology can assist predicting coronary atherosclerosis progression, and with further development could become an integral part of high-risk patient care.

Materials

Name	Company	Catalog Number	Comments
ANSYS Workbench (version 19.0)	ANSYS		Commercial finite element solver
MATLAB (version 2019b)	Mathworks		Commercial programming platform
MicroDicom/ImageJ	MicroDicom/ImageJ		Open Source DICOM reader
Visual Studio (version 2019)	Microsoft		Commercial Integrated Development Environment