The influence of constitutive law choice used to characterise atherosclerotic tissue material properties on computing stress values in human carotid plaques

Calculating high stress concentration within carotid atherosclerotic plaques has been shown to be complementary to anatomical features in assessing vulnerability. Reliability of stress calculation may depend on the constitutive laws/strain energy density functions (SEDFs) used to characterize tissue material properties. Different SEDFs, including neo-Hookean, one-/two-term Ogden, Yeoh, 5-parameter Mooney–Rivlin, Demiray and modified Mooney–Rivlin, have been used to describe atherosclerotic tissue behavior. However, the capacity of SEDFs to fit experimental data and the difference in the stress calculation remains unexplored. In this study, seven SEDFs were used to fit the stress–stretch data points of media, fibrous cap, lipid and intraplaque hemorrhage/thrombus obtained from 21 human carotid plaques. Semi-analytic solution, 2D structure-only and 3D fully coupled fluid-structure interaction (FSI) analyses were used to quantify stress using different SEDFs and the related material stability examined. Results show that, except for neo-Hookean, all other six SEDFs fitted the experimental points well, with vessel stress distribution in the circumferential and radial directions being similar. 2D structural-only analysis was successful for all seven SEDFs, but 3D FSI were only possible with neo-Hookean, Demiray and modified Mooney–Rivlin models. Stresses calculated using Demiray and modified Mooney–Rivlin models were nearly identical. Further analyses indicated that the energy contours of one-/two-term Ogden and 5-parameter Mooney–Rivlin models were not strictly convex and the material stability indictors under homogeneous deformations were not always positive. In conclusion, considering the capacity in characterizing material properties and stabilities, Demiray and modified Mooney–Rivlin SEDF appear practical choices for mechanical analyses to predict the critical mechanical conditions within carotid atherosclerotic plaques.


Introduction
Carotid atherosclerotic disease is responsible for around 15-20% of all ischemic strokes (Brott et al., 2011), with the luminal stenosis being the only validated diagnostic criterion for patient risk stratification. However, this criterion becomes less reliable in patients with mild to moderate carotid stenoses (Barnett et al., 1998). Increasing evidence has suggested that both the physical characteristics of atherosclerotic plaques and the mechanical loading within the structure may allow greater potential to predict clinical progression than luminal stenosis alone. A vulnerable carotid atherosclerotic plaque is characterized by the presence of intraplaque hemorrhage (IPH) and a large lipid-rich necrotic core, with symptomatic plaques also showing evidence of fibrous cap (FC) rupture. These features have been shown to predict future events in both symptomatic (Altaf et al., 2008;Eliasziw et al., 1994) and asymptomatic (Singh et al., 2009;Takaya et al., 2006) patients. As plaques are continually subject to mechanical loading due to pulsatile blood pressure and flow, FC rupture is thought to occur when loading exceeds its material strength (Richardson et al., 1989;Tang et al., 2009a). FC stress can differentiate symptomatic from asymptomatic patients Zhu et al., 2010) and both plaque deformation and FC stress have been found to be Contents lists available at ScienceDirect journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com associated with subsequent cerebrovascular ischemic events in symptomatic patients Teng et al., 2011Teng et al., , 2013. There is therefore a need to integrate both plaque morphological and compositional features with the critical mechanical conditions for improved patient risk stratification. However, the reliability of re-predicting the critical mechanical conditions is largely dependent on the accuracy of plaque geometry and the material properties of each atherosclerotic component, including intra-plaque hemorrhage or thrombus (IPH/T), lipid and FC. Accurate reconstruction of plaque geometry is critically dependent on the imaging technique utilized, limited by resolution and tissue discrimination. The behavior of material properties is determined from experimental measurements, but also influenced by constitutive laws. Several potential strain energy density functions (SEDFs) can be used to characterize the material, such as neo-Hookean (Akyildiz et al., 2011;Caille et al., 2002;Lee et al., 1996;Ohayon and Tracqui, 2005), one-term Ogden (Barrett et al., 2009), two-term Ogden (Li et al., 2006(Li et al., , 2007Tang et al., 2008;Versluis et al., 2006), Yeoh (Cunnane et al., 2015;Lawlor et al., 2011), fiveparameter Mooney-Rivlin (Gao and Long, 2008;Maher et al., 2009), Demiray (Chau et al., 2004;Delfino et al., 1997;Kaazempur-Mofrad et al., 2003a), and modified Mooney-Rivlin SEDF (Tang et al., 2009a(Tang et al., , 2013(Tang et al., , 2014b. These seven SEDFs have been used in numerous studies to model the mechanical behavior of carotid atherosclerotic plaques. However, the effectiveness of each SEDF at characterizing the material properties of carotid atherosclerotic tissue and the resulting variance in predicted critical mechanical conditions within the plaque structure using these different SEDFs remain unexplored. In this study, the seven selected SEDFs are used to fit the experimental data obtained from uniaxial extension tests performed on human carotid atherosclerotic tissues. The accuracy in computing the mechanical stress within the plaque structure is assessed by using analytical solutions, idealized 2D structure-only and 3D fully coupled fluid-structure interaction (FSI) simulations and the related material stability is discussed.

Strain energy density functions
Human carotid atherosclerotic tissues exhibit non-linear stress-strain behavior at low stretch levels (Maher et al., 2009;Mulvihill et al., 2013;Teng et al., 2014c). These complexities need to be accommodated by specific SEDFs/hyperelastic material models. Several SEDFs were adopted/developed for this purpose, including neo-Hookean, one-and two-term Ogden, Yeoh, 5-parameter Mooney-Rivlin, Demiray and modified Mooney-Rivlin models, with details as follows: neo-Hookean model

5-parameter Mooney-Rivlin model
modified Mooney-Rivlin model Ogden material models are expressed in terms of principal stretches, λ j j ¼ 1; 2; 3 ð Þ , while the others are expressed in terms of invariants of Cauchy-Green deformation tensor. I 1 ¼ J À 2=3 I 1 and I 2 ¼ J À 4=3 I 2 with I 1 and I 2 being the first and second invariant of the unimodular component of the left Cauchy-Green deformation tensor, Þis the principal stretch. J ¼ det F ð Þ and F is the deformation gradient. κ is the Lagrangian multiplier for the incompressibility and the remainder are material constants which will be determined by fitting experimental measurements.

Material testing and data fitting
Endarterectomy carotid plaque samples from 21 symptomatic patients were collected during surgery. The local ethics committee approved the study protocol and all patients gave written informed consent. Details of tissue preparation, testing protocol and equipment used have been described previously (Teng et al., 2014c). In total, stress-stretch curves were obtained successfully from 65 media strips from 17 samples, 59 FC strips from 14 samples, 38 lipid strips from 11 samples and 21 IPH/T strips from 11 samples. An energy-based average strategy (Teng et al., 2014c) was used to obtain the representative stress-stretch curve for each atherosclerotic tissue as shown in Fig. 1 and Fig. S1 in the Supplemental material.
Cauchy stress in terms of principal stretches can be obtained from each SEDF, where W is the part in SEDFs without the incompressible term, κ J À 1 ð Þ. In the case of simple uniaxial extension with an incompressible tissue strip, The stress-stretch relationship can be therefore derived, and material constants can be obtained by minimizing the following objective function, The relative error is used to assess the fitting quality, in which σ 11 and σ e 11 are the predicted and measured stress, respectively; and N is the number of data points. In this study, all material constants were constrained to be positive as one or some negative material constants might lead to unphysical phenomena, e.g., an increased stretch leads to a decreased stress.

Material stability
The material stability should be taken into account when a SEDF is used to describe the material properties and to calculate mechanical conditions within the plaque. The material stability is material constant-and loading-dependent (Adina, 2013;Ogden, 2003).
Convexity is one of the criteria for assessing the material stability to some extent defined as, implying the stress being a monotonic increasing function of the stretch ratio. If Eq.
(4) holds for all λ i 40; W is globally strictly convex. However, convexity depends on the measures being employed, such as stretch ratio, true strain and Green strain. Convexity in one measure does not necessarily guarantee material stability, but failure of convexity may have undesirable consequences for the development of numerical schemes (Ogden, 2003). The material stability can also be partially characterized by using stability curves under certain loading conditions (Adina, 2013). Considering an incompressible solid that undergoes homogeneous deformations, the equilibrium requires the equality of the external and internal virtual work as, in which R i is the external loading. For a stable material, the work done by external loadings should always be positive, which implies the minimum eigenvalue of the matrix below must be positive, where δ ij is the Kronecker delta. The material stability can be represented by the curve of the minimum eigenvalue versus the stretch ratio or strain. Due to the involvement of external loading in the matrix, the stability curve can only be drawn under specific situations of homogeneous deformation, including uniaxial and biaxial extension and pure shear deformation.

Stress calculation
The difference in stress calculation using different SEDFs is quantified under the scenario of a long circular tube with a thick wall, using 2D idealized structureonly plaque models and 3D fully coupled FSI simulation with an idealized plaque geometry. The following considerations apply: there is a semi-analytic solution for calculating stress in a statically pressurized long circular tube so the result can always be obtained precisely regardless of material stability; and 2D (Kaazempur-Mofrad et al., 2003b;Li et al., 2006Li et al., , 2007Sadat et al., 2011Sadat et al., , 2010Zhu et al., 2010) and 3D FSI (Bluestein et al., 2008;Gao et al., 2011;Huang et al., 2014a;Leach et al., 2010;Tang et al., 2009a) analyses have been widely used in re-predicting critical mechanical conditions within carotid plaques.

Long circular tube with a thick wall
Under static pressure, the deformation of a long tube can be solved with a plane strain assumption. The equilibrium equation and boundary conditions are, where r i and r o are the inner and outer radii, respectively, σrr and σ θθ are the radial and circumferential stresses, respectively and p i is the applied internal pressure. The pressure can be expressed in terms of stress by integrating the equilibrium equation, Therefore r i can be obtained by solving the above equation numerically with the assumption of incompressibility, and stress in both circumferential and radial directions can be further computed.   2.4.2. Idealized plaque models As shown in Fig. 2, the model represents a typical atherosclerotic plaque of 50% stenosis, composed of FC, lipid, and IPH (the rest was assumed to be media), with FC thickness being 600 mm. The stenosis of 50% was chosen to represent lesions with moderate luminal stenosis (30%-69%) for which the optimal treatment strategy remains unclear. FC thickness of 600 mm was the mean minimum value quantified by in vivo high-resolution magnetic resonance imaging in symptomatic patient groups (Sadat et al., 2009). The plaque length was set to be 20 mm and the lengths of proximal and distal section were 120 mm to avoid a potential entrance effect when FSI simulations were performed. The cross section at the most stenotic site ( Fig. 2A) was used for the 2D structure-only analysis.

Media
For 2D models, the plaque structure was meshed using quadrilateral elements and the symmetry in geometry was considered when the fixity was applied. The pressure waveform at the inlet shown in Fig. 3 was used as the loading condition. For 3D fully coupled FSI analyses, a volume curve-fitting technique was employed, in which, the 3D plaque domain was divided into hundreds of small "volumes" to curve-fit the irregular plaque geometry with plaque component inclusions (Tang et al., 2009b). The entire plaque domain and fluid volume were meshed using hexahedral elements with 136,080 elements for the structure and 108,864 elements for the fluid, respectively. The symmetry in geometry was also considered when the fixity was applied (Fig. 2B). The pressure waveforms at both inlet and outlet are shown in Fig. 3, with systolic and diastolic pressure at the inlet being 120 and 80 mmHg, respectively. The number of time steps was set to be 200. In 3D FSI analyses, the blood flow was assumed to be Newtonian, viscous and incompressible. FSI simulations were performed using Adina 9.0.3 (Adina Inc., MA, USA). For the fluid domain, the slipping line was the central line of the fluid domain and leader-follower pairs were specified. The energy convergence criterion was used for solid domain during equilibrium iterations with the relative energy tolerance being 0.05 and relative force and moment tolerance being 0.01. For the fluid domain, the relative tolerances for velocities, pressure and displacements were set to be 0.06 to control the equilibrium. The fluid-structure coupling was solved iteratively. Both the displacement and velocities at the fluid-structure interface and the forces on the structure due to the viscous fluid were checked for convergence. Relative displacement/velocity and force tolerances were both set to be 0.06. The maximum principal stress (Stress-P 1 ) over the diseased region at systole was used to characterize the critical mechanical conditions. The influence of using different SEDFs on Stress-P 1 was subsequently analyzed.

Stress-stretch curves fitting
As shown in Fig. 1, except for the neo-Hookean model (Fig. 1A), all SEDFs could accurately characterize the experimental data points. The detailed fitted constants for each SEDF are listed in Table 1. In addition, apart from neo-Hookean, one-term Ogden and Demiray material models, fitting results of the rest SEDFs were sensitive to initial guessed values. For example, apart from the sets of constants of the modified Mooney-Rivlin listed in Table 1, those listed in Table 2 could also well characterize the stress-stretch relationships of FC, media, Lipid and IPH/T as shown in Fig. S2 in the Supplemental material.

Stress distribution in the wall of a long circular tube
From a tube with the inner and outer radii of 4 and 5 mm, respectively, with internal pressure of 16 kPa, the inner and outer radii under deformed configuration can be obtained by solving Eq. (6) numerically under the assumption of plane strain. The stress distribution across the wall thickness can be computed by using Eq. (1). As shown in Fig. 4, except for the neo-Hookean model, the other six models have a similar prediction of stress in both circumferential and radial directions. The circumferential stress at the inner boundary is 80.79 for neo-Hookean, 130.36 for one-term Ogden, 130.33 for two-term-Ogden, 131.40 for Yeoh, 127.70 for 5parameter Mooney-Rivlin, 114.33 for Demiray and 114.33 for modified Mooney-Rivlin (unit: kPa). When neo-Hookean is excluded, the variation of circumferential stress at the inner boundary was 8.3%.

Idealized plaque models
Differences in stress prediction became more prominent when 2D structure-only analysis with idealized plaque models were performed, as shown in Fig. 5. The stress value predicted by using neo-Hookean was the lowest (243.27 kPa), while those obtained from others were much higher (one-term Ogden: 820.11 kPa; twoterm Ogden: 819.23 kPa; Yeoh: 662.43 kPa; 5-parameter Mooney-Rivlin: 703.35 kPa; Demiray: 732.88 kPa and modified Mooney-Rivlin: 732.34 kPa) with a variation of 11.1% (the value from neo-Hookean is excluded).
Consistent with previous reports, we again observed that 2D structure-only analysis overestimated the stress prediction, when compared with 3D FSI analysis (Huang et al., 2014b). As shown in Fig. 6, in 3D FSI analysis the high stress concentration in the shoulder region was much lower than that observed in 2D modeling (Fig. 5). Successful simulations were only achieved when neo-Hookean, Demiray and modified Mooney-Rivlin SEDFs were used, while the solution procedure was interrupted when a certain internal pressure loading level was reached when one-term Ogden, two-term Ogden, Yeoh and 5-parameter Mooney-Rivlin models were used.

Material stability
The material stability is determined by SEDF type (Adina, 2009), the relationship among material constants in each SEDF (Adina, 2009;Ogden, 2003) and the initial-/boundary-value problems (Adina, 2013;Zheng, 2008). Considering these complexities, it is nearly impossible to validate the stability comprehensively. However, material stability can, in part, be examined by analyzing the energy contour convexity and stability curves under homogeneous deformations.
As an example, the energy contour of FC is shown in Fig. 7 (those of media, lipid and IPH/T can be found in Figs. S3-S5 in the Supplemental material). The contours of neo-Hookean, Yeoh, Demiray and modified Mooney-Rivlin models are strictly convex, while those of one-term Ogden, two-term Ogden and 5-parameter Mooney-Rivlin are not. The stability curve of FC of each SEDF is shown in Fig. 8 (those of media, lipid and IPH/T are presented in Figs. S6-S8 in the Supplemental material). In certain stretch levels, the minimum eigenvalue of the matrix, shown in Eq. (5), becomes negative or approaches to zero in SEDFs of one-term Ogden, twoterm Ogden and 5-parameter Mooney-Rivlin.

Discussion
Calculation of structural stress has been shown to be complementary to anatomic determinants in assessing the vulnerability of atherosclerotic Teng et al., 2014a) and aneurysmal (Khosla et al., 2014) lesions. However, the accuracy and reliability of the calculation depend on a number of factors  including the resolution and tissue discrimination of the imaging modality, modeling strategy (Huang et al., 2014b), the boundary/ loading conditions and the constitutive laws used to describe tissue material behavior. We believe that this is the first study to: (1) assess the difference between constitutive laws (SEDFs) in characterizing experimental data obtained from uniaxial extension tests with human carotid atherosclerotic tissue; (2) quantify the difference in stress concentrations within the plaque structure using different SEDFs; and (3) characterize the material stability of different SEDFs. In total, seven SEDFs which had been previously used in calculating structural stress in the carotid artery were tested in this study. Results obtained indicate that, except for neo-Hooken model, all other six models could fit the experimental data appropriately (Fig. 1). However, the fitted results were sensitive to the initial guessed values for two-term Ogden, Yeoh, 5-parameter Mooney-Rivlin and modified Mooney-Rivlin models (Tables 1 and  2  (2) it is difficult to interpret the physical meaning of each material constant in these SEDFs. However, such non-uniqueness appears not to have any perceptible effect on the final stress calculations. For the case of modified Mooney-Rivlin model, when the constants listed in Table 2 rather than Table 1 were used, the computed high stress concentration was 208.46 kPa, which is nearly identical as the one shown in Fig. 6C computed using constants listed in Table 1.
As listed in Table 1, except for neo-Hookean model, the relative error was similar for the same type of tissue when different SEDF was used. However, subtle differences were observed when the fitted curve and experimental data points were compared, in particular, at a low stretch range (Fig. 1). These subtle differences may result in discrepant values when predicting the stress distribution across the vessel wall when using a semi-analytic approach (Fig. 4) and the high stress concentration in the shoulder region in 2D structure-only analyses (Fig. 5). Further analyses indicated that SEDFs from the same family (Ogden: oneterm and two-term Ogden models; polynomial Mooney-Rivlin: Yeoh and 5-parameter Mooney-Rivlin models; and Mooney-Rivlin with exponential term: Demiray and modified Mooney-Rivlin models) had a similar capacity in characterizing the stress-stretch behavior of each tissue type and an overall similar calculation of stress (Figs. 4 and 5). Theoretically, SEDF with any combination of stretch ratios can be developed, but not all of them can be incorporated into numerical schemes leading to a successful solution to solve initial-/boundary-value problems. In this study, successful 3D FSI simulations were obtained only with neo-Hookean, Demiray and modified Mooney-Rivlin models. The failure of one-term Ogden, two-term Ogden and 5 parameter Mooney-Rivlin models can be explained partially by the non-strictly convex energy contours (Fig. 7) and stability curves (Fig. 8). However, it is worth noting that convexity of energy contours (strain measure-dependent) and positive stability indicators (loading-dependent) do not guarantee successful 3D analyses, as exampled by Yeoh model for the idealized geometry used in this study (Fig. 2). A successful 3D fully coupled FSI analysis depends on many factors, including mesh quality, settings of time steps and interaction boundaries, etc., and it may not be reasonable to conclude that the material instability definitely accounts for the failure of the analysis with Yeoh model. However in this study, for all 3D analyses except for the difference in SEDF, all of these other factors that may influence stability were unchanged. A successful 3D structure-only analysis was achieved using Yeoh model when the pressure at the inlet (Fig. 3) was used as the loading condition applying on the entire inner surface, whilst the pressure gradient along the model was ignored. It is, therefore, possible that under a certain deformation condition, Yeoh model with material constants listed in Table 1 becomes unstable. 3D structure-only analyses with one-term Ogden, twoterm Ogden and 5 parameter Mooney-Rivlin models were not successful and it was successful with neo-Hookean, Demiray and modified Mooney-Rivlin models. For those that failed in 3D analyses, the failure persisted despite efforts to adjust mesh density, time function, time steps and leader-follower settings. It should be noted that material stability under a general loading condition is complicated (Ogden, 2003) and the energy convexity shown in Fig. 7 and stability curves shown in Fig. 8 are insufficient to fully describe this characteristic.
Considering the poor capacity of neo-Hookean models to characterize the material properties of atherosclerotic tissues and the large deviations in predicting stress concentrations, this model may not be appropriate for mechanical analysis for carotid atherosclerotic plaques. Instead, both the Demiray and modified Mooney-Rivlin models may be more appropriate choices considering their capacity to characterize material properties of atherosclerotic tissues and the convergence in 2D structure-only and 3D FSI analyses. However, compared with modified Mooney-Rivlin model, Demiray model appeared to have a relatively poorer capacity in fitting experimental data from aortic tissues (Teng et al., 2015). The efficiency of modified Mooney-Rivlin models in calculating critical mechanical conditions have been validated in numerous patient-specific 2D structure-only , 3D structure-only (Huang et al., 2014b;Teng et al., 2011), 3D oneway FSI (Huang et al., 2014b) and 3D full coupled FSI (Tang et al., 2009a(Tang et al., , 2014b analyses. In this study, only the first invariant of the deformation gradient, I 1 , was included. It is always possible to include the second invariant, I 2 , although the involvement of I 2 did not improve the fitting quality as listed in Table 3. Moreover, the predicted stress values with I 2 involvement were nearly identical in 3D FSI analyses (Fig. 6D).
In this study, Stress-P 1 was used as the stress measure. The von Mises stress has also been widely used to assess the critical mechanical conditions within atherosclerotic plaques, as summarized in a recent review (Holzapfel et al., 2014). The von Mises stress, σ v , is an important measure charactering material yielding due to excessive shear stresses, in which σ i (i¼1, 2, 3) stands for the principal stress in the ith direction and Stress-P 1 ¼max [σ i (i¼1, 2, 3)]. In the idealized plaque model used in this study, the distribution of Stress-P 1 and von Mises stress was similar with less than 5% difference in calculating the high stress concentration in the shoulder region ( Fig. 9). In general, Stress-P 1 is often used to assess the failure of brittle materials and may better govern plaque structure failure, as tensile and compressive stresses frequently coexist due to geometrical complexity ( Fig. 9A and C). However, debate continues regarding which stress measure is optimal to characterize the critical mechanical conditions within an atherosclerotic plaque (Holzapfel et al., 2014).
There are some limitations of the current study, (1) atherosclerotic tissues are fiber-oriented and their anisotropic material properties were not considered. This may result in the stress levels reported in this study being underestimated; (2) no attempt was made to model residual stresses and this may lead to an overestimation of stress levels (Delfino et al., 1997); (3) calcium was not included in idealized plaque models; (4) the convergence of FSI depends on geometry, mesh, settings of loading steps and interaction boundaries, amongst other parameters. Thus, it is virtually impossible to ensure that failure was due to material instability alone, excluding all other factors; and (5) different commercial numerical packages may use different strategies to solve linear equations and handle the iteration and convergence. The conclusions obtained in this study are based on Adina 9.0.3 (Adina Inc., MA, USA) and may not be valid when other numerical packages are used. Table 3 The fitted material parameters when I 2 was included in the modified Mooney-Rivlin model.

Disclosure
The authors do not have any conflict of interest to be declared.