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.

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 ...

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 methodolog...

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