MRI-based patient-specific human carotid atherosclerotic vessel material property variations in patients, vessel location and long-term follow up

Background Image-based computational models are widely used to determine atherosclerotic plaque stress/strain conditions and investigate their association with plaque progression and rupture. However, patient-specific vessel material properties are in general lacking in those models, limiting the accuracy of their stress/strain measurements. A noninvasive approach of combining in vivo 3D multi-contrast and Cine magnetic resonance imaging (MRI) and computational modeling was introduced to quantify patient-specific carotid plaque material properties for potential plaque model improvements. Vessel material property variation in patients, along vessel segment, and between baseline and follow up were investigated. Methods In vivo 3D multi-contrast and Cine MRI carotid plaque data were acquired from 8 patients with follow-up (18 months) with written informed consent obtained. 3D thin-layer models and an established iterative procedure were used to determine parameter values of the Mooney-Rivlin models for the 81slices from 16 plaque samples. Effective Young’s Modulus (YM) values were calculated for comparison and analysis. Results Average Effective Young’s Modulus (YM) and circumferential shrinkage rate (C-Shrink) value of the 81 slices was 411kPa and 5.62%, respectively. Slice YM value varied from 70 kPa (softest) to 1284 kPa (stiffest), a 1734% difference. Average slice YM values by vessel varied from 109 kPa (softest) to 922 kPa (stiffest), a 746% difference. Location-wise, the maximum slice YM variation rate within a vessel was 311% (149 kPa vs. 613 kPa). The average slice YM variation rate for the 16 vessels was 134%. The average variation of YM values for all patients from baseline to follow up was 61.0%. The range of the variation of YM values was [-28.4%, 215%]. For plaque progression study, YM at follow-up showed negative correlation with plaque progression measured by wall thickness increase (WTI) (r = -0.7764, p = 0.0235). Wall thickness at baseline correlated with WTI negatively, with r = -0.5253 (p = 0.1813). Plaque burden at baseline correlated with YM change between baseline and follow-up, with r = 0.5939 (p = 0.1205). Conclusion In vivo carotid vessel material properties have large variations from patient to patient, along the diseased segment within a patient, and with time. The use of patient-specific, location specific and time-specific material properties in plaque models could potentially improve the accuracy of model stress/strain calculations.

Considerable efforts have been made by several research groups to quantify mechanical material properties of atherosclerotic vessels. Holzapfel et al. used cyclic quasistatic uniaxial tension tests in axial and circumferential directions for different components of atherosclerotic lesions from human atherosclerotic iliac arteries [27]. Anisotropic and highly nonlinear tissue components properties as well as considerable interspecimen differences were identified by the experimental data of individual samples. The circumferential direction of the fibrous cap demonstrated the lowest fracture stress (254.8 +/-79.8 kPa at stretch 1.182 +/-0.1) of all intimal tissues. The adventitia exhibited the highest and the nondiseased media the lowest mechanical strength on average [27]. Different layers displayed different direction-dependent mechanical behaviors, which are of vital importance to the realistic computational models and accurate stress/strain prediction [16,[27][28]. Williamson et al. reported that artery wall stresses have low sensitivities for material properties variability [29]. In a histology-based model, a +/-50% variation in elastic modulus was found to lead to less than a 10% change in artery wall stress at the site of rupture [29]. Tang et al. investigated the effects of plaque structure and material properties on stress behaviors in human atherosclerotic plaques by using 3D FSI models, and reported that the softer materials result in higher stress values [30]. Teng et al. performed uniaxial tests using different plaque components, axial and circumferential oriented adventitia, media and intact specimens prepared from human carotid arteries [31][32]. Most of the material properties investigations in the literature used ex vivo specimens and in vitro experimental techniques. Nieuwstadt et al. performed a numerical feasibility study for carotid plaque elasticity estimation using ultrasound elastography, MRI, and inverse FEA [24]. They were able to estimate material coefficients for carotid intima and lipid-rich necrotic core (lipid) with histology validations. Smoljkić et al. proposed a non-invasive, energy-based assessment of patient-specific material properties of arterial tissue [25]. Their results showed that imposing conditions on strain energy can provide a good estimation of carotid material properties from the non-invasively measured pressure and diameter data. Czernuszewicz et al. performed some preliminary study of non-invasive in vivo characterization of human carotid plaques with acoustic radiation force impulse ultrasound. Their method was able to differentiate soft tissues from stiffer tissues with histological validations [26]. In vivo vessels material properties and patient-specific material properties are still scarce in the current literature.
For plaque models based on in vivo image data, shrinkage rate (axial shrinkage and circumferential shrinkage) should be determined so that the in vivo vessel geometry could be shrunk to its "no-load" shape from which loaded stress/strain conditions could be obtained. Tang et al. and Yang et al. introduced a shrink-stretch process in their 3D fluid-structure interaction (FSI) models with in vivo plaque data [16,[33][34][35]. Their shrink-stretch process include: a) axial and circumferential shrink the vessel to get the numerical starting geometry; and b) axial stretch and circumferential pressurization were applied to recover the vessel in vivo shape. The axial shrinkage was 9% so that the in vivo vessel length could be recovered when a 10% axial stretch were applied. The Circumferential shrinkages of vessel lumen and outer wall were determined so that: a) the vessel volume was conserved; b) vessel shape after pressurization and 10% axial stretch recovered the original in vivo shape. Speelman et al. and Gee et al. also showed the importance and necessity of the pre-shrink process with their computed tomography (CT)-based simulations for abdominal aortic aneurysm [36][37]. Huang et al. reported a non-uniform shrink-stretch process which had better match with the in vivo vessel geometries [38]. Liu et al. introduced a non-invasive approach to quantify patient-specific vessel material properties and plaque circumferential shrinkage rate between vessel in vivo and "no-load" geometries [39]. Their material properties and circumferential shrinkage rate were calculated by 2D plaque models. Their results showed that effective Young's Modulus (YM) from the 12 human carotid arteries varied from 137 kPa to 1435 kPa and vessel circumferential shrinkage to "no-load" condition varied from 6% to 32%. Overall, quantified patient-specific shrinkage rate using in vivo data are rare in the current literature.
In this paper, a non-invasive approach [39] of combining 3D multi-contrast MRI, in vivo Cine MRI and computational 3D thin-layer model [40] was used to quantify patient-specific carotid plaque material properties and circumferential shrinkage rates. The variation of the material properties among different subjects, within each diseased artery, and with time was investigated.

In vivo serial MRI data acquisition and segmentation
Serial MRI data of carotid atherosclerotic plaques from 8 patients (5 male, 3 female; age: 62-83, mean = 71; see Table 1 for details) were acquired at the University of Washington (UW), Seattle by the Vascular Imaging Laboratory (VIL) using protocols approved by the UW Institutional Review Board and with written informed consent obtained. For each patient, MRI slices at baseline (Time 1, T1) and follow-up (Time 2, T2, Scan time intervals were about 18 months) were matched up using vessel bifurcation, stenosis features and with careful review by the MRI group. For simplicity, the vessel segment assembled using the selected slices is called the plaque. Cuff systolic and diastolic arm pressure was recorded for modeling use. In vivo Cine and three-dimensional (3D) multi-contrast MR images of the carotid arteries were acquired using a 3.0T whole-body scanner (Philips Achieva, R2.6.1, Best, The Netherlands) and a dedicated 8-channel, phased array carotid coil. The carotid bifurcation was located on two-dimensional (2D) TOF (Time of Flight) and oblique black blood MR images. A 3.5cm region centered on the carotid bifurcation was imaged by high-resolution axial bright and black blood imaging. Detailed data acquisition and segmentation procedures were published before and are omitted here [10,39]. For each patient, locations with Cine sequence and nearly-circular slices from 3D MRI were selected for calculating the material parameter values in the modified Mooney-Rivlin model used in our previous publications.   16 re-constructed carotid plaques based on in vivo MRI data from 8 patients, each with 2 time points. Fig 2 gives a sample plaque with 3D MRI slices and matching Cine lumen contours at multiple locations corresponding to minimum and maximum pressure conditions. The Cine and 3D MRI sequences were designed with equal slice distance and matched using longitudinal position from the bifurcation.

The 3D thin-layer model to determine material parameter values
A 3D thin-layer modeling approach introduced by Huang et al. [40] was used to determine material parameter values in our selected material model. For every slice that Cine data was available, a thin slice thickness was added to make a 3D thin-layer model (see Fig 3). Then our established 3D model construction and mesh generation procedures with axial and circumferential shrink were applied in the material value determination process [40]. The carotid artery was assumed to be hyperelastic, isotropic, incompressible and homogeneous. The nonlinear modified Mooney-Rivlin (M-R) model was selected to describe the material properties of the vessel wall [41][42]. The strain energy function was given by: where C = [C ij ] = X T X is the right Cauchy-Green deformation tensor; I 1 and I 2 are the invari- h i is the deformation gradient; c 1 , c 2 , D 1 and D 2 form the material parameter set. The modified Mooney-Rivlin model was selected because it was able to fit carotid artery vessel properties measured by uniaxial and biaxial mechanical testing data and good agreement was obtained [9,43]. With our limited pressure data (only two pressure data points: maximum and minimum pressures), material parameter values could not be uniquely determined. According to our previous experiences [10,16], we set c 2 = 0 and D 2 = 2. The c 1 / D 1 value was adjusted. For lipid cores, we used parameter values in our previous publications: c 1 = 0.5 kPa, D 1 = 0.5 kPa, D 2 = 0.5 [10]. For each 3D thin-layer model, a 10% axial shrinkage rate was applied. Then an iterative procedure [39] was followed to adjust the parameter values in the modified M-R model and the circumferential shrinkage rate to match both maximum and minimum Cine lumen circumferences corresponding to systolic and diastolic pressures, respectively (see Fig 4 for

details).
Material parameters values c 1 (36.6 kPa) and D 1 (14.4kPa) from our previous paper [10] were used as initial estimates for each slice. The initial lumen circumferential shrinkage rate S 1 was set to be 10%. The details of the iteration procedure are given in

The 2D models for comparison purpose
To compare our 3D thin-layer model approach with 2D model approach, parameter values for Modified Mooney-Rivlin (MR) models and circumferential shrinkage values for the 81 slices were also calculated using 2D modeling method [29] with the iterative procedure given in Fig  4. The main difference between 2D model and 3D thin-layer model is that 2D model does not have axial stretch while 3D thin-layer model does which makes the 3D thin-layer model closer to full 3D models. Their differences are reported in Section 3.5.

Using an effective Young's modulus as a stiffness indicator for material properties
The stress-stretch relationship for the Mooney-Rivlin model is given by: where σ is Cauchy stress, and λ is stretch ratio. In order to facilitate comparison, it is easier to use a single parameter to compare vessel stiffness from different patients or slices. The effective Young's modulus (YM) E for the stretch ratio interval [1.0, 1.3] is defined as: The least-squares technique was used to calculate the YM values that best fit the M-R model.

Statistical method
For all correlation study related variables in Section 3.4, Shapiro-Wilk normality test showed no evidence of non-normal distribution. Standard student t-test was performed for possible correlations between plaque progression and material stiffness, plaque burden and wall thickness using our patient follow-up data. We also applied the Kolmogorov-Smirnov test to check the normality. A p-value 0.81 was obtained, still indicating no evidence of non-normality. In case the normality assumption was a concern, we calculated the correlation p-values under the non-parametric Spearman rank test (see 3.4), which does not require the normality assumption. The results were consistent with those by the parametric t-test.  Quantifying in vivo carotid material property variations in patients, vessel location and time lumen circumferences (Cir -Max and Cir -Min ), relative lumen circumferences variation rates (δ -Cir ), maximum and minimum lumen pressure (P), effective Young's modulus (YM) and circumferential shrinkage rate (C-Shrink) were provided by Table 2. The relative lumen circumferences variation rate (δ -Cir ) is defined as:

Overview of material parameter values, effective YM and circumferential shrinkage rates
The averaged YM value of the 81 slices was 411kPa, which was consistent with the current literature [10-11, 28, 32]. The averaged C-Shrink and δ -Cir value of the 81 slices was 5.62%,

Plaque material properties have large patient-to-patient variations
Stress-Stretch Ratio curves from Mooney-Rivlin models for the 16 plaque samples using average parameter values of their slices are presented in Fig 6. The average YM values and circumferential shrinkage (C-Shrink) values from 16 plaque samples were given in Table 3. The average YM values for the stiffest plaque sample (P16) was 922 kPa, 746% higher than that for the softest plaque (P13, YM = 109 kPa). This showed that plaque material properties have large variations from patient to patient and patient-specific material properties should be used in plaque models. It should be noted that pressure difference had a big effect on the material properties of vessels. Average C-shrink value from the 16 samples was 6.29%. The softest sample had 21.7% C-shrink, while the stiffest sample had a negative C-Shrink value (-1.78%). Negative C-Shrink means to obtain the zero-load geometry of the 3D thin-layer model, the in vivo slice lumen needed to expand slightly so that it could regain the in vivo circumference when 10% axial stretch and pressure were applied. Axial stretch makes the vessel to shrink in radial direction.

Plaque material properties change over time
We investigated whether vessel stiffness changed over time, and how material stiffness changes were associated with plaque progression and other risk factors. The effective Young's modulus (YM), wall thickness (WT) [44], plaque burden (PB) and circumferential shrinkage (C-shrink) values at baseline (T1) and follow-up (T2) for the 8 patients are given in Table 5. The variation of YM over time (δ YM-Time ), the difference of PB (δ PB ) and the WT increase (WTI) are given by: where YM T1 , PB T1 and WT T1 are the YM, PB and WT values at T1, and YM T2 , PB T2 and WT T2 are the YM, PB and WT values at T2, respectively. The average δ YM-Time , δ PB and WTI values in Table 5 Table 5, it has a Shapiro-Wilk p-value 0.17, indicating no evidence of non-normality. We also applied the Kolmogorov-Smirnov test to check the normality. We got p-value 0.81, still indicating no evidence of non-normality. When the Kolmogorov test is applied, data have to be re-scaled based on its mean and standard deviation. If the YM data were not re-scaled, the Kolmogorov test would give a small p-value. Since our data size is small and the normality assumption could be a concern, we also reported the correlation p-values under the nonparametric Spearman rank test, which does not require the normality assumption. The results are consistent with those by the parametric t-test.

The differences between 2D model and 3D thin-layer model
To compare our 3D thin-layer model approach with 2D model approach, parameter values for Modified Mooney-Rivlin (MR) models and circumferential shrinkage values for the 16 plaque samples from 8 patients were also calculated using 2D modeling method [29] with the iterative procedure given in Fig 4. The average effective YM and C-Shrink values for the 16 plaque samples calculated by the 2D model and 3D thin-layer model are given in Table 6. The variation of YM (δ -YM-Model ) and the difference of C-Shrink (δ C-Shrink ) between the 2 models are given by: where YM 2D and C-Shrink 2D are the YM and circumferential shrinkage values, respectively, calculated by the 2D model, and YM 3D and C-Shrink 3D are the YM and circumferential shrinkage values, respectively, calculated by the 3D thin-layer model. YM values from the 3D thin-layer models were 7.98% higher than that from the 2D models. The range of δ -YM-Model values was [5.19%, 9.04%]. C-Shrink values from the 3D thin-layer models were 4.85% lower than that from the 2D models. The range of δ -C-Shrink values was [-5.44%, -3.29%]. These differences showed that using 3D thin-layer model to replace 2D model would lead to improved accuracy of material parameter value estimations.

Significance of in vivo patient-specific and location-specific vessel material properties
Most of the research on determining arterial wall material properties has been performed using ex vivo specimens and in vitro experimental techniques. In vivo estimation of patientspecific material properties is scarce, which is a serious limitation for patient-specific plaque models. A noninvasive approach of combining in vivo Cine and 3D MRI and simple 3D thinlayer modeling was introduced to quantify patient-specific and location-specific vessel material properties and improve model prediction accuracies. Our results from 16 plaques and 81 slices showed that slice YM values could vary from 70 kPa to 1284 kPa, 16 times of the lowest YM value. YM values from slices of the same vessel could vary by 311% (Table 4). Future studies should render 3D plaque models using patient-specific and location-specific material properties to quantify their impact on stress/strain calculations.

Significance of the 3D thin-layer model
Full 3D computational models are essential for plaque stress and strain calculations. They should also be used in the iteration process to determine the material parameter values that fit the variation in luminal circumference with the cardiac cycle, as provided by imaging techniques such as MR Cine data. However, full 3D model construction is time consuming. The material parameter adjusting process takes 20-50 iterations. The labor cost is unbearable if full 3D models were used. Previously, 2D models were used as an approximation [39]. We proposed a 3D thin-layer modeling method as an approximation with much less computational cost than full 3D models. 3D axial stretch was applied to the 3D thin-layer models, mimicking full 3D model procedures. The 3D thin-layer model was closer to the real in vivo conditions than used 2D models. This 3D thin-layer approach needs less computational cost than a full 3D model to achieve convergence for both circumferential shrinkage rate and material parameter value. Compared to 2D model results, the YM values from the 3D thin-layer models were 7.98% higher. The difference showed that using 3D thin-layer model to replace 2D model would lead to improved accuracy of material parameter value estimations. Note: c 2 = 0 and D 2 = 2 for all cases. Due to axial shrink applied to the 3D thin-layer model, some C-Shrink values in 3D thin-layer model were negative. Both averages by patients and by slices are given. https://doi.org/10.1371/journal.pone.0180829.t006

Vessel material properties change over time and plaque component material properties
Plaque material properties changed over time with a 61.0% average variation of YM among all patients. It is clearly of great importance to be able to quantify patient-specific, location-specific and time-specific vessel material parameter values under in vivo conditions. However, there are some limitations as the available Cine MR data did not provide plaque component size change information needed to calculate plaque component material parameter data. In this study, lipid core parameter values were extracted from previous publications [10]. New imaging protocols or in vitro experimental techniques could be used to obtain patient-specific plaque component properties.

Impact of hypertension
It is well-accepted that hypertension is a common thread in most cardiovascular diseases. Table 2 shows that the six patients who had higher than 140 mmHg systolic pressure had higher YM values. The average YM value for the six patients was 612 kPa, whereas the average YM value for the remaining patients was 263 kPa. Since we do not have direct pressure measurements at present time, we plan to investigate the impact of hypertension when better pressure data becomes available. And the number of patient is the important limitation for the correlation prediction between hypertension and vessel material.

Model limitations
Cine MRI was used to determine vessel material parameter values, matching in vivo plaque geometries under both systolic and diastolic pressure conditions. Cine MRI is widely accepted to acquire time-dependent vessel motion and deformation. Multi-layer structure and anisotropic material properties of arteries were not considered since MRI does not provide layer information. Cine data provided only circumference variations under cardiac pressure. Another limitation was that location-specific pressure measurement was not available. Currently, arm cuff pressure values are used in most image-based studies. Noninvasive acquisition of intraplaque pressure data remains a challenge.