*In vivo* serial MRI data acquisition and segmentation

After informed consent, serial MRI data of carotid atherosclerotic plaques from 16 patients (15 male, 1 female; age: 59–81, mean=71.9; two patients had L/R carotid MRI usable) were acquired 2–4 times (scan interval: 18 months; 4 scans: 7; 3 scans 8; 2 scans: 1) by the University of Washington (UW) Vascular Imaging Laboratory (VIL) using protocols approved by the UW Institutional Review Board. A total of 38 scan pairs (baseline and follow-up) were formed for progression analysis. MRI scans were conducted on a GE SIGNA 1.5-T whole body scanner using an established protocol (Yuan and Kerwin, 2004). Multi-contrast images in T1, T2, proton density (PD), time-of-flight (TOF), and contrast-enhanced (CE) T1 weighted images of atherosclerotic plaques were generated to characterize plaque tissue composition and luminal and vessel wall morphology [14–16]. Segmentation was performed using a custom-designed computer package CASCADE (Computer-Aided System for Cardiovascular Disease Evaluation) developed by the Vascular Imaging Laboratory at the University of Washington. The slice thickness was 2 mm. Field of view = 160 mm × 160 mm. Matrix size 512 × 512 (the real matrix size was 256 × 256. Images were machine interpolated to 512 × 512). After interpolation, the in-plane resolution was 0.31 × 0.31 mm^{2}. Figure 1 gives two examples re-constructed from MRI data showing plaque progression and regression, respectively.

### 3D geometry re-construction and mesh generation

Under in vivo conditions, arteries are axially stretched and pressurized. Therefore, in vivo MRI plaque geometry needs to be shrunk axially and circumferentially a priori to obtain the no-load starting geometry (see Figure 2) for computational simulations as our model starts from the no-load geometry with zero stress/strain, zero pressure and no-flow conditions. The shrinkage in axial direction was 9% so that the vessel would regain its *in vivo* length with a 10% axial stretch. Circumferential shrinkage for lumen (about 8-12%) and outer wall (about 2-5%)was determined by trial-and-error so that: 1) total mass of the vessel was conserved; 2) the loaded plaque geometry after 10% axial stretch and pressurization had the best match with the original *in vivo* geometry.

Because advanced plaques have complex irregular geometries and 3D FSI models involve large deformation and large strain, a geometry-fitting mesh generation technique was developed to generate mesh for these models [17, 18]. Using this technique, the 3D plaque and fluid domains were divided into hundreds of small “volumes” to curve-fit the irregular plaque geometries (see Figure 3). Computational meshes for these volumes were then generated automatically by ADINA (ADINA R & D, Inc., Watertown, MA, USA), a commercial finite element software used to solve these FSI models.

### 3D fluid–structure interaction plaque model and solution methods

3D FSI models were constructed using well-established procedures for the 38 pairs (a total of 54 plaque models) to obtain plaque wall stress (PWS) and flow shear stress (FSS) for correlation analysis [17, 18]. The artery wall was assumed to be hyperelastic, isotropic, incompressible and homogeneous. The nonlinear modified Mooney-Rivlin model was used to describe the material properties of the vessel wall [19, 20]. The strain energy function was given by,

\mathrm{W}={\mathrm{c}}_{1}\left(\phantom{\rule{0.25em}{0ex}}{\mathrm{I}}_{1}\u20133\right)+{\mathrm{c}}_{2}\left(\phantom{\rule{0.25em}{0ex}}{\mathrm{I}}_{2}\u20133\right)+{\mathrm{D}}_{1}\left[\mathrm{exp}\left({\mathrm{D}}_{2}\left(\phantom{\rule{0.25em}{0ex}}{\mathrm{I}}_{1}\u20133\right)\right)\u20131\right],

(1)

{\mathrm{I}}_{1}={\displaystyle \sum {\mathrm{C}}_{\mathrm{ii}},\phantom{\rule{1em}{0ex}}{\mathrm{I}}_{2}}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\phantom{\rule{0.12em}{0ex}}\left[{\mathrm{I}}_{1}{2}^{}\u2013{\mathrm{C}}_{\mathrm{ij}}{\mathrm{C}}_{\mathrm{ij}}\right],

(2)

where I_{1} and I_{2} are the first and second strain invariants, **C** =[C_{ij}] = **X**^{T}**X** is the right Cauchy-Green deformation tensor, c_{i} and D_{i} are material parameters chosen to match experimental measurements and the current literature [21, 22]. Parameter values used in this paper: vessel/fibrous cap, *c*_{1}=36.8 KPa, *D*_{1}=14.4 KPa, *D*_{2}=2; calcification, c_{1}=368 KPa, D_{1}=144 KPa, D_{2}=2.0; lipid-rich necrotic core, *c*_{1}=2 KPa, *D*_{1}=2 KPa, *D*_{2}=1.5; loose matrix, *c*_{1}=18.4 KPa, *D*_{1}=7.2 KPa; *D*_{2}=1.5. *c*_{2} = 0 was set for all materials [17].

Blood flow was assumed to be laminar, Newtonian, viscous and incompressible. The incompressible Navier–Stokes equations with arbitrary Lagrangian–Eulerian (ALE) formulation were used as the governing equations. A no-slip condition, natural traction equilibrium boundary condition and continuity of displacement were assumed on the interface between solid and fluid. Inlet and outlet were fixed (after initial pre-stretch) in the longitudinal (axial) direction, but allowed to expand/contract with flow otherwise. Patient-specific systolic and diastolic pressure conditions from the last hospital admission were used as the maximum and minimum of the imposed pulsatile pressure waveforms at the inlet and outlet of the artery. Details of the FSI model were given in Tang, et al. (2004, 2009) [17, 20].

The 3D FSI models were solved by ADINA, using unstructured finite element methods for both fluid and solid domains. Nonlinear incremental iterative procedures were used to handle fluid–structure interactions. The governing finite element equations for both solid and fluid models were solved by Newton–Raphson iteration method. More details of the computational models and solution methods can be found in Tang et al. (2004, 2009) and Bathe (2002) [18–20]. Plaque wall stress and flow shear stress data corresponding to peak systolic pressure were recorded for analysis.

### Plaque progression measurements and data extraction for correlation analysis

For each scan pair, slices from the baseline (Time 1, or T1) and follow-up (Time 2, or T2; baseline-follow up pairs can also be T2-T3 and T3-T4) scans were matched using the carotid bifurcation as the registration reference. Only matched common carotid artery (CCA) and internal carotid artery (ICA) slices were chosen for analysis. For each matched slice, 100 evenly-spaced points from the lumen were selected and vessel wall thickness, PWS, and FSS from 3D FSI model solutions at each point for baseline and follow-up were obtained for analysis. For the 38 pairs, 400–1000 matched data points for each plaque were obtained for correlation analysis. Plaque progression at each data point was expressed by vessel wall thickness increase (WTI) defined as

\mathrm{WTI}=\mathrm{Wall}\phantom{\rule{0.25em}{0ex}}\mathrm{Thickness}\phantom{\rule{0.25em}{0ex}}\mathrm{at}\phantom{\rule{0.25em}{0ex}}\mathrm{follow}-\mathrm{up}\u2013\mathrm{Wall}\phantom{\rule{0.25em}{0ex}}\mathrm{Thickness}\phantom{\rule{0.25em}{0ex}}\mathrm{at}\phantom{\rule{0.25em}{0ex}}\mathrm{baseline}.

(3)

### Statistical analysis

The Linear Mixed-Effects (LME) models [23] were used to study the correlation between WTI and the predictors (PWS and FSS) at the initial time point of each time pair. The measured points are nodes, and it seems reasonable to capture the dependence among nodes based on their 3-dimensional spatial locations on the vessel. The models are specified as follows.

For individual patients, the LME model was defined as

{y}_{\mathit{ij}}={\beta}_{0}+{\beta}_{1}{x}_{\mathit{ij}}+{\in}_{\mathit{ij}},

(4)

where *y*_{
ij
} is the WTI value at the *i* th node on the *j* th slice, *x*_{
ij
} is the corresponding value of FSS (or PWS). *β*_{0} and *β*_{1} are the fixed effects of the predictor for the baseline and the changing rate of WTI, respectively. The spatial dependence structure among *y*_{
ij
} is accounted for by the exponential isotropic variogram model in 3D space, which is analogous to the autoregressive model in 1D space [23]. Specifically, the vector of random error terms (*∈*_{
ij
}) follows a joint Gaussian distribution with mean 0, and the correlation between {\in}_{{i}_{1}}{{j}_{1}}_{} and {\in}_{{i}_{2}}{{j}_{2}}_{} is assumed an exponential function *ϕ*^{s}, where *s* is the Euclidean distance between the 3-dimensional spatial locations of the two nodes on the vessel, *ϕ* is the correlation parameter to be estimated in the model fitting by restricted maximum likelihood (REML) algorithm. The null hypothesis that no correlation exists between WTI and FSS (or PWS) indicates *β*_{1} = 0. A small p-value of the Student’s t-test for the coefficient [23] provides a strong statistical evidence to reject the null and accept the existence of correlation.

To study the correlation between WTI and FSS (or PWS) based on the cohort of patients, the patients were grouped together according to the first (T1-T2), the second (T2-T3), and the third (T3-T4) time-pairs. For all available patients at a particular time pair, the LME model is defined as

{y}_{\mathit{ijk}}={\beta}_{0}+{\beta}_{1}{x}_{\mathit{ijk}}+\phantom{\rule{0.5em}{0ex}}{b}_{k}+{\in}_{\mathit{ijk}},

(5)

where *y*_{
ijk
} is the WTI value at the *i* th node on the *j* th slice of the *k* th patient, *x*_{
ijk
} is the corresponding value of FSS (or PWS). *β*_{0} and *β*_{1} have the same meanings as those in equation (4). For the *k* th patient, random effect *b*_{
k
} follows a Gaussian distribution with mean 0, which models the clustering dependence of WTI values within this patient. The vector of error terms (*∈*_{
ijk
}) follows a joint Gaussian distribution with mean 0. The patients are assumed to be independent so the correlation between the error terms from different patients is 0. The correlation between {\in}_{{i}_{1}}{{{j}_{1}}_{}}_{k} and {\in}_{{i}_{2}}{{j}_{2}k}_{} from the same *k* th patients is assumed an exponential function {\varphi}_{k}^{s}, where *s* is the Euclidean distance between the 3D spatial locations, *ϕ*_{
k
} is the patient-specific correlation parameter to be estimated in model fitting. To test the correlation, we again use the p-value for testing *β*_{1} = 0.

To assess the above 3D spatial LME models, we randomly permuted the response WTI values and then calculated the p-values. Such a permutation breaks potential correlations between WTI and FSS (or PWS), so the corresponding p-value is expected to follow a Uniform (0, 1) distribution. Our results verified that the empirical distributions of these p-values after permutations were close to such expectation, which evidenced that our models have the type I error well controlled at the level of individual hypothesis test. That is, these models are appropriate in fitting the correlations between WTI and FSS (or PWS), and the observed small p-values are reliable to indicate the significance of the correlations. We controlled the type I error rate at 0.05.

To measure the direction and magnitude of the linear correlations between WTI and FSS (or PWS), we applied a dependence-adjusted correlation coefficient *r* based on the LME models to account for the spatial dependence structure among the nodes by the exponential isotropic variogram model in 3D space [23], which is described above for models (4) and (5). This *r* is analogous to the Pearson’s correlation coefficient that is valid only for independent observations. Specifically, note that the Pearson’s correlation coefficient between *x* and *y* can be written as:

\widehat{\mathit{cor}}\left(x,y\right)=\widehat{\beta}\sqrt{\frac{\widehat{var}\left(x\right)}{\widehat{var}\left(y\right)}},

(6)

where \widehat{\mathrm{var}}\left(x\right) and \widehat{\mathrm{var}}\left(y\right) are the sample variances, and \widehat{\beta} is the estimated slope coefficient by fitting *x* to *y* with a simple regression that requires independence among the observations. Now, since the observations of FSS (or PWS) are dependent, Pearson’s correlation coefficient cannot be directly used for measuring the correlation between FSS (or PWS) and WTI. However, following the same idea as Pearson’s correlation coefficient, we can apply a dependence-adjusted correlation coefficient *r* by

r={\widehat{\beta}}_{1}\sqrt{\frac{\widehat{\mathit{var}}\left(x\right)}{\widehat{\mathit{var}}\left(y\right)}},

(7)

where {\widehat{\beta}}_{1} is the estimated slope coefficient by fitting FSS (or PWS) to WTI with the LME model (4) or (5). Because {\widehat{\beta}}_{1} is estimated in the models that have adjusted the 3D spatial dependence structure among the observations of FSS (or PWS), *r* is the dependence-adjusted correlation coefficient.