Are two-dimensional images sufficient to assess the atherosclerotic plaque vulnerability: a viscoelastic and anisotropic finite element model

The correlation between plaque morphology as well as composition and plaque vulnerability has been the motivation for many recent studies. In a generic point of view, instability of atherosclerotic plaques is known to be the result of a thin fibrous cap and a large and highly compliant necrotic core area. There have been numerous two-dimensional (2D) and three-dimensional (3D) computational models mostly based on the finite element method (FEM) to assess the plaque vulnerability. It is well-known that 2D models are not reliable as they do not provide a consistent assessment on the vulnerability of plaques and are highly erroneous. 3D models offer a more effective evaluation but creating 3D models to be further assessed by computational means such as the finite element method is time-consuming. However, 2D models are easier to develop and are less time-consuming to assess. In this study, we propose a novel computational platform by which the plaque vulnerability is assessed using only 2D plaque models. We develop idealistic 2D models and their corresponding idealistic 3D models. The idealistic 3D models resemble the worstand best-case scenarios for each 2D model. Using these 3D idealistic models, a standard error (SE) is estimated and then added to the peak stress values calculated earlier using 2D models. These SEs are also used to assess the probability of plaque stability. In this platform, the effect of viscoelasticity and anisotropy of the plaque composition is taken into consideration and the transmural pressure considered is similar to that of physiological conditions (dynamic pressure). The current study may suggest a more realistic insight to the prediction of atherosclerotic plaques rupture using 2D images.


Introduction
Although plenty of studies have been performed on the detection and therapy of instable plaques, thrombus mediated ischemic cardiovascular disease still remains the leading cause of death in developed countries [1]. The rupture of unstable atherosclerotic plaque in coronary artery regularly leads to a significant number of ischemic cardiovascular events. Conventional procedures for the prediction of rupture based on imaging of the plaque morphology and composition [2,3] still provides rather inaccurate and insufficient predictors of risk [4]. The challenge for imaging methods is that prediction of the coronary plaque rupture requires not only an accurate quantification of fibrous cap thickness [5] and necrotic core morphology [6], but also a precise knowledge of the mechanical properties of the arterial wall and plaque components at any given stage of the atherosclerotic progression [7][8][9][10].
Several reports indicated that such vulnerable plaques can be detected clinically by various techniques, including Intravascular Ultrasound (IVUS) [11,12], Optical Coherence Tomography (OCT) [13,14] Computed Tomography (CT) [15] doi: 10.7243/2052-4358-3-3 and Magnetic Resonance Imaging (MRI) [16]. Both IVUS and OCT are used to make measurements for lesion length and lumen size, but OCT is shown to be more accurate. With intravascular imaging, OCT offers clear, photographic quality images, as opposed to grainy, lower-resolution IVUS images. Also, it is known that OCT has a resolution 10 times greater than IVUS (OCT being 10 μm vs. IVUS, 100 μm) [13]. Moreover, OCT can only penetrate about 2-3 mm into the vessel wall, where as IVUS's penetration is 4-8 mm. If a vessel has significantly remodeled due to plaque burden, the outline of the true lumen disappears on the OCT image and the vessel outer boundary on the IVUS images.
Finite element methods (FEM) have been extensively used as a strong tool to assess the effects of mechanical factors on instability of plaque rupture [10,[17][18][19][20][21]. Depending on the assumption made, all of the models provided are accurate to some extent. This is because there are some major issues associated with the applied models such as: (1) the constitutive material models applied are either isotropic linear, isotropic nonlinear, or anisotropic, (2) viscous properties of atherosclerotic plaque is either overlooked or somewhat oversimplified, and (3) mostly the modeling is implemented for steady-state conditions. The main aims of this study are to develop a novel computational platform which improves the estimation of mechanical stresses applying on the fibrous cap under actual physiological loading conditions. In this study, the fibrosis is assumed to be viscoelastic and the transmural pressure is considered to be dynamic. The solver used in this study is ANSYS, R15.0, a commercial FEM code.

Methods
This study has two phases: phase I: working on idealistic models and phase II: working on real data or data obtained from patients.

Phase I
(1) 3 idealistic 2D models are developed to comprehensively define all types of plaque morphologies. As the 3D structure of these 2D models is unknown, for each idealistic 2D model, (2) 2 idealistic 3D models are developed. These 3D models are designed to provide the worst and best case scenarios for stress distribution if the 2D model were hypothetically to be considered 3D. (3) The 2D models are solved for stress analysis. (4) The 3D models are used to define lower and upper limits of the stress values obtained earlier. These limits are implemented to evaluate the plaque stability for the 2D model.

Phase II
(1) The 2D models are obtained from patients are solved for stress analysis. (2) The similar idealistic 2D model to the real 2D model is chosen. The idealistic model is chosen based on the geometry of the plaques where the fibrous cap is minimal, i.e., the most critical region. (3) The lower and upper limits obtained by the idealistic 3D models are implemented for the real model and then the probability of plaque stability is evaluated. Each stage is described below in more detail.

Idealistic 2D models
Using a comprehensive range of clinical data, a variety of authentic morphologies for plaque composition was carefully investigated. The objective of this investigation was to define principal and idealistic morphology to represent the geometry of all types of necrotic cores, fibro-athermatous caps, and arterial remodelling indexes [16,17]. Four representative IVUS images were carefully chosen based on the geometry of the thin fibrous cap which was then categorized to: (1) Nodal ( Figure 1A), linear (Figures 1B and 1C), and curve-linear ( Figure 1D) [14]. This classification was the foundation for the generation of 3 idealistic 2D models shown in Figure 1E (nodal), Figure 1F (linear) and Figure 1G (curve-linear). Any other types of 2D geometries of plaque can be created by one or a combination of two of these idealistic models.

Idealistic 3D models
In this section, each of the idealistic 2D models is considered separately, and then for each model, 2 further idealistic 3D models are developed such that the cross section of the plaque in the 3D models are identical. These 2 idealistic 3D models resemble the worst and best case scenarios for stress distribution, if the 2D model were hypothetically to be considered 3D. The stress distribution values obtained from the 2D model are between the lower and upper limits obtained from the idealistic 3D models. Using these lower and upper limits, a standard error (SE) is defined and added to the results obtained from 2D models. Figures 2A-3D

Actual 2D imaging data obtained from real patients
In this section, 4 images (2D models) obtained from patients diagnosed with atherosclerosis are considered the subject of the study. These images are developed by two well-established imaging techniques, IVUS and OCT. Given that the IVUS technique provides images deep in the vessel wall but with less resolution compared to the OCT technique which provides a higher resolution but with less penetration, a co-registered image of these two techniques may provide a realistic geometry of the plaque composition including the fibrous cap area, the necrotic core area, and calcific nodules (if present).

Constitutive model and boundary conditions
A cylindrical coordinate is defined in which z denotes the axial direction of the artery and r and θ indicate the radial and circumferential directions, respectively, of the cross section of a typical plaque (Figure 1). Given that the fibrous cap is considered anisotropic, it is assumed that mechanical properties of fibrosis are identical in the axial (z) and radial (r) directions but different in the circumferential (θ) direction (much stronger) (Table 1) [23]. Er is the modulus of elasticity  Figure 1E; (B,E) are the idealistic 3D models for Figure 1F; and (C,F) are the idealistic 3D models for Figure 1G. This classification is based on the geometry of the fibrous cap which is considered to be: nodal (A), linear (B,C), and curve-linear (D), and the corresponding idealistic models are defined as: E (nodal), F (linear) and G (curve linear) [22].
Fibrosis 100 1000 0.01 0.27 500 Table 1. Mechanical properties of fibrosis [15].  [23], and (D) The stress distribution and the peak cap stress obtained in the current study. For the current model, the PCS is 378 kPa which is 2 kPa larger than the computed PCS in that study. This small discrepancy may be because of the differences in the mesh models.
in the radial direction, Ez is the modulus of elasticity in the axial direction, Eθ is the modulus of elasticity in the circum ferencial direction, v rθ , v rz , and v θz are the Poisson's ratios. Necrotic core or lipid pool was assumed to be isotropic with elastic modulus of 1kPa [18] and the Poisson's ratio of 0.49 (in compressibility) [23]. In order to define viscoelastic behavior for fibrosis, a time Prony series with its known constants was implemented ( Table 2) [24]. This time Prony series has five spring-damper elements i.e; where relaxation times spanning five decades (τ1=0.001 s, τ2=0.01 s, τ3=0.1 s, τ4=1.0 s, and τ5=10.0 s) were selected. These numerical values were chosen in order to define completely the time scales of the dynamic studies. The relaxation times, τi, where i=1, …, 5, are associated with the damping coefficient η and spring constant E in each Maxwell element which is defined as τi=ηi/Ei where i=1, …, 5. The CAD models of these three plaque morphologies were created in Solid Works 2014, and then these models were converted to ANSYS design files through ANSYS Design Module 2015. Material properties and suitable boundary conditions (fixed all around) were assigned to the models. In this study, the transmural pressure was assumed to be dynamic following the physiological conditions ( Table 3) [25]. For comparison 5 ( / ) between 2D and 3D models, von-Mises stress distribution (with assumption of plane strain state for 2D models) at t=4.2 s has been considered. Preliminary results show that the maximum stress occurs at t=0.2 s but since heart beats almost 100,000 times a day, the maximum value at t=0.2 s seems to be numerical artifact, thus, t=4.2 s (cycle 4) is considered for the peak stress values in which the peak dynamic stress remained unchanged (more detail is provided in the result section). For FE computation, models were meshed by ANSYS Mechanical 2015. As the stress map in the fibrous cap area is important, the mesh density in this area was selected to be high. The solver implemented in this study was Mechanical APDL 2015 running on an Intel® Core™ 2 Due T6670 @ 2.2GHz and 2.00 GB of RAM.

Results and discussion Validation
In order to validate our proposed computational platform, one of the results in the study by Ohayon et al., [23] was reproduced and resolved. In this model, the peak cap stress (PCS) ofthe model with the same morphology, mechanical properties, blood pressure, etc., was calculated and compared. Results indicate that the PCS obtained by our solver (378 kPa) and that of Ohayon et al., (376 kPa) shows 2 kPa discrepancy(less than 1% error) (Figure 3). This small discrepancy could be because of small differences between the two mesh models.

Mesh independency study
Mesh independency of our results is of particular significance in this study which is done by performing further computation for each plaque models with higher mesh density. The mesh size is decreased until the point where by increasing the mesh density the results are not improved. The result of mesh independency study is outlined in Table 4.

Case study I-effect of viscoelasticity of fibrous cap
The idealistic model shown in Figure 1E is considered. The fibrous cap thickness was set to 95 µm (unlike 70 µm [22]), the diameter of artery is assumed to be 3.3 mm, and the lumen area was considered elliptic with diameters of 250 μmx200 μm, representing a typical plaque.
In Figure 4, the PCS values for 2D models (Figure 1E) are shown at different times. As observed, the maximum PCS occurs at t=0.

Case study II-the assessment of the PCS in the idealistic 2D models
The idealistic 2D models represented in Figures 1E-1G are solved for the value of PCS in the fibrous cap area (Figure 5). Results indicate that the maximum PCS (336 kPa) occurs in model 1E in which the fibrous cap geometry is focal. It follows with the model with a linear fibrous cap area ( Figure 1F) in which the PCS is estimated to be 299 kPa (11% decrease) and then the curve-linear fibrous cap geometry (Figure 1G) with the PCS being 216 kPa (36% decrease). The area with the maximum PCS in all models are located in the vicinity of the NC and not the lumen, which is consistent with previous studies [22]. Also, numerical values for the PCS are consistent with those obtained by Mohammadi et al.,

Case study III-the assessment of the PCS in the idealistic 3D models
In this study, two idealistic 3D models (Figure 2) are developed for each idealistic 2D model (Figures 1E-1G). Both these 3D models have the same geometry in the cross sectional area to that of the 2D model. The objective for emerging these models is that although they share the same cross sectional area of the 2D models, they are designed to yield the minimum and maximum PCS in the fibrous cap area. The stress distribution in the plaque model for these 6idealistic 3D models is shown in Figure 6.
The range of the PCS for the 3D models corresponding to Figure 1E is 99 kPa to 412 kPa, whereas the PCS in the 2D model is 336 kPa. Since the range of 300 kPa-99 kPa=201 kPa (lower limit) is greater than the range of 412 kPa-300 kPa=112 kPa (upper limit), the 2D plaque is more likely stable. The range of PCS for the 3D models corresponding to Figure 1F is 176 kPa to 385 kPa, whereas the PCS in the 2D model is 299 kPa. Since the range of 300 kPa-176 kPa=124 kPa (lower limit) is greater than the range of 385 kPa-300 kPa=85 kPa (upper limit), the 2D plaque is more likely stable. The range of PCS for the 3D models corresponding to Figure 1G is 144 kPa to 351 kPa, whereas the PCS in the 2D model is 216 kPa. Since the range of 300 kPa -144 kPa=156 kPa (lower limit) is significantly greater than the range of 351kPa-300 kPa=51 kPa (upper limit), the 2D plaque is likely stable. Results are outlined in Table 5, in which the discrepancy of the lower and higher limits of PCSs is compared with that of obtained from 2D models. The traditional way of evaluation of the vulnerability of plaques which is based on 2D models, easily  Figures 1E -1G recommend of stability of plaques shown in Figures 1F and 1G and instability of plaque shown in Figure 1E. Figure 7 shows the range of safe value (PCS<300 kPa) for PCS (green bar) and range of PCS (yellow bar) obtained using viscoelastic material model for all models. In order to assess the risk of plaque rupture, the length of the yellow bar in and out of the green zone is implemented such that:

Case study IV-effect of viscoelastic fibrosis on the PCS
In this step, all nine 3D and 2D models are solved for PCS using elastic and viscoelastic material models listed in Tables 1 and 2, respectively. Results outlined in Table 6 show that in all models the PCS obtained using elastic material model is higher than   Table 6. Results of PCS computation in models with elastic and viscoelastic material. As it can be seen in the table, PCS in models with elastic material is higher than PCS in models with Viscoelastic material. PCS at t=4.2s has been considered for this comparison. *S for sphere-like models and C for cylinder-like models.
that of obtained using viscoelastic material model. For all 2D models this discrepancy is between 30-32%, whereas for 3D models this discrepancy falls within the range of 23-52%.

Case study V-2D data obtained from patients
In this section, 2D geometries of the atherosclerotic plaque obtained from 4patients were developed using the combination of OCT and IVUS techniques (co-registered image) as shown in Figures 8A-8D. In these 2D models, a clear geometry of the fibrous cap and the necrotic core is available. These 2D models are solved using the proposed computational platform for the PCS. Figure 8A1 shows the stress distribution in the entire tissue of the model shown in Figure 8A. In this model, the cap thick-ness is measured to be 0.118 mm and the PCS is estimated to be 355 kPa. By comparing this 2D model and the idealistic 2D models developed earlier, it is seen that this plaque has morphology similar to model 1G ( Figure 1G) in general, and in particular has a similar morphology to model 1E (Figure 1E) where the cap thickness is minimal (0.118 mm). In order to choose similar idealistic models to an intended actual data obtained from patients, the area in which the cap thickness is minimal is considered. Following the idealistic models and their corresponding 3D models, the probable range for the PCS for such geometry is 106 kPa to 437 kPa and the probability of the plaque stability is calculated to be =59%.

Application on data obtained from patients
The conventional assessment applied for this plaque suggests that the plaque will definitely fail, but our proposed computational platform suggests that since the 3D structure 306-106 437-106 doi: 10.7243/2052-4358-3-3 Figure 7. The range of safe value for PCS (green bar) and the total range of PCS (yellow bar). Comparison between the length of yellow bar outside and inside of green bar provides practical information about probability of rupture in a plaque. In all of these models, plaque is more likely to be stable since the probability of stability is greater than that of instability. of the plaque is not known, there is a 59% chance that the plaque may stay stable. Figure 8B1 shows the stress distribution in the entire tissue regarding the model shown in Figure 8B. In this model, the cap thickness is measured to be 0.067 mm and the PCS is estimated to be 493 kPa. Comparing this 2D model and the idealistic 2D models developed earlier, this plaque has morphology similar to model 1F (Figure 1F). Following the idealistic models and their corresponding 3D models, the probable range for the PCS for such geometry is 0.6 PCS-1.28 PCS which is 296 kPa-631 kPa, and the probability of the plaque stability is calculated to be =1%.
The conventional assessment applied for this plaque suggests that the plaque will definitely fail, and our proposed computational platform reinforces this regardless what the 3D structure of the plaque might be. Figure 8C1 shows the stress distribution in the entire tissue for the model shown in Figure 8C. In this model, the cap thickness is measured to be 0.141 mm and the PCS is estimated to be 270 kPa. After comparing this 2D model and the idealistic 2D models developed earlier, this plaque has morphology similar to model 1E ( Figure 1E) idealistic models and their corresponding 3D models, the probable range for the PCS for such geometry is 81.1 kPa to 332 kPa and the probability of the plaque stability is calculated to be =87%. The conventional assessment suggests that the plaque is most likely unstable, but our proposed computational platform suggests that since the 3D structure of the plaque is not known to us, there is a chance of 13% that the plaque may fail. Figure 8D1 shows the stress distribution in the entire tissue regarding the model shown in Figure 8D. In this model, the cap thickness is measured to be 0.106 mm and the PCS is estimated to be 178 kPa. When comparing this 2D model and the idealistic 2D models developed earlier, it is seen that this plaque has morphology similar to model 1G ( Figure 1G). Following the idealistic models and their corresponding 3D models, the probable range for the PCS for such geometry is 119 kPa to 290kPa and the probability of the plaque stability is calculated to be >0%. The conventional assessment applied for this plaque suggests that the plaque is definitely stable and our proposed computational platform proves the stability of this plaque.

Conclusion
In this study, a novel computational platform was developed by which the vulnerability risk of an atherosclerotic plaque is evaluated by using only 2D models. We developed idealistic 2D models and corresponding 3D models (two 3D models for each) which provide the worst and best case scenarios for each 2D model. This is because the 3D structure of the plaque is not known to us, so defining the worse and best case scenarios may be helpful. The PCS is assessed in the 2D model first (denoted as PCS 2D ) and then the idealistic 3D models provide its corresponding lower and upper limits. This lower and higher limits are considered here as the standard errors (SE). Using the obtained SE, the probability of the plaque stability is assessed. In a real application and for a given 2D model obtained from patients, (1) The similar idealistic 2D model to the 2D model obtained from patients is found, (2) The lower and upper limits (SE) corresponding to the idealistic 2D model is assigned, (3) the 2D model is solved for its PCS 2D , 0.59xPCS 2D 1.28xPCS 2D

Figure 1G
0.67xPCS 2D 1.63xPCS 2D Table 7. The platform by which the PCS 2D is used for plaque stability risk assessment. The lower and upper limits are calculated using idealistic 3D models explained earlier.
models developed in this study, the lower and upper limits are defined as outlined in Table 7.
The computational platform seems to be useful. For example, in model shown in Figure 8B, the cap thickness is measured to be 0.067 mm and the PCS is estimated to be 493 kPa. Conventional evaluation models would suggest that this plaque is unstable, which is completely consistent with our results (99%, unstable). For the plaque composition, a suitable viscoelastic and anisotropic material model was considered and the transmural pressure was defined to be similar to that of physiological condition (dynamic condition). Results strongly suggest that the viscoelastic properties of the plaque composition are of particular importance for the assessment of the PCS. If elastic material model is implemented instead, for 2D models it may contribute to up to 32% error and for 3D models this error may be as high as 52%. Even when viscoelastic material model is used, the difference between the PCSs obtained in the first and the fourth cycles (when PCS is settled) is almost 15%.

Limitations
It should be noted that without information such as spatial resolution along all three directions of IVUS/OCT or other images, it is not possible to estimate how accurate the plaque models are. These parameters are available in the Digital Imaging and Communications in Medicine(DICOM) headers of images.
In order to improve the current work, more data analysis is required. More data obtained from patients would be helpful to be analyzed. Developing more realistic 2D and 3D models and comparing them with results obtained from 3D imaging data obtained from patients could help more strengthen the computational platform developed in this study.