Computer simulation of three-dimensional plaque formation and progression in the carotid artery

Abstract

Atherosclerosis is becoming the number one cause of death worldwide. In this study, three-dimensional computer model of plaque formation and development for human carotid artery is developed. The three-dimensional blood flow is described by the Navier–Stokes equation, together with the continuity equation. Mass transfer within the blood lumen and through the arterial wall is coupled with the blood flow and is modeled by a convection–diffusion equation. The low-density lipoproteins transports in lumen of the vessel and through the vessel tissue are coupled by Kedem–Katchalsky equations. The inflammatory process is modeled using three additional reaction–diffusion partial differential equations. Fluid–structure interaction is used to estimate effective wall stress analysis. Plaque growth functions for volume progression are correlated with shear stress and effective wall stress distribution. We choose two specific patients from MRI study with significant plaque progression. Plaque volume progression using three time points for baseline, 3- and 12-month follow up is fitted. Our results for plaque localization correspond to low shear stress zone and we fitted parameters from our model using nonlinear least-square method. Determination of plaque location and composition, and computer simulation of progression in time for a specific patient shows a potential benefit for the prediction of disease progression. The proof of validity of three-dimensional computer modeling in the evaluation of atherosclerotic plaque burden may shift the clinical information of MRI from morphological assessment toward a functional tool. Understanding and prediction of the evolution of atherosclerotic plaques either into vulnerable or stable plaques are major tasks for the medical community.

Appendix: Growth functions and the fitting procedure

For simulation of the plaque growth, the following procedure was used.

Step 1: Start from the original in geometry at T1;

Step 2: Set

$$ \begin{gathered} f_{t1\_0} \left( {i,j} \right) = f_{T1} \left( {i,j} \right) \hfill \\ f_{t1\_1} \left( {i,j} \right) = f_{T1} \left( {i,j} \right) + \left( {f_{T2} \left( {i,j} \right) - f_{T1} \left( {i,j} \right)} \right)/m \hfill \\ \tau_{t1\_0} \left( {i,j} \right) = \tau_{T1} \left( {i,j} \right) \hfill \\ \tau_{t1\_1} \left( {i,j} \right) = \tau_{T1} \left( {i,j} \right) + \left( {\tau_{T2} \left( {i,j} \right) - \tau_{T1} \left( {i,j} \right)} \right)/m \hfill \\ \sigma_{t1\_0} \left( {i,j} \right) = \sigma_{T1} \left( {i,j} \right) \hfill \\ \sigma_{t1\_1} \left( {i,j} \right) = \sigma_{T1} \left( {i,j} \right) + \left( {\sigma_{T2} \left( {i,j} \right) - \sigma_{T1} \left( {i,j} \right)} \right)/m \hfill \\ \end{gathered} $$

Step 3: We use m time steps to go from T1 to T2 and n time steps to go from T2 to T3. This means that we use total $n + m$ time steps to go from T1 to T3. For $k = 1, \ldots n + m$ do the following:

• GF1:

  $$ f_{t1\_(k + 1)} \left( {i,j} \right) = a_{0} \left( j \right) + a_{1} \left( j \right) \cdot \left( {w(j) \cdot f_{t1\_k} \left( {i,j} \right) + \left( {1 - w(j)} \right) \cdot f_{t1\_0} \left( {i,j} \right)} \right) + a_{2} \left( j \right) \cdot \left. {\frac{{{\text{d}}f}}{{{\text{d}}t}}} \right|_{{T_{1\_k} }} (i,j) \cdot \Updelta t_{k} $$

• GF2:

  $$ \begin{aligned} f_{t1\_(k + 1)} \left( {i,j} \right) = & a_{0} \left( j \right) + a_{1} \left( j \right) \cdot \left( {w(j) \cdot f_{t1\_k} \left( {i,j} \right) + \left( {1 - w(j)} \right) \cdot f_{t1\_0} \left( {i,j} \right)} \right) \\ + a_{2} \left( j \right) \cdot \left. {\frac{{{\text{d}}f}}{{{\text{d}}t}}} \right|_{{T_{1\_k} }} \left( {i,j} \right) \cdot \Updelta t_{k} + a_{3} \left( j \right) \cdot \tau_{t1\_k} \left( {i,j} \right) + a_{4} \left( j \right) \\ \cdot \left. {\frac{{{\text{d}}\tau }}{{{\text{d}}t}}} \right|_{{T_{1\_k} }} \left( {i,j} \right) \cdot \Updelta t_{k} \\ \end{aligned} $$

• GF3:

  $$ \begin{aligned} f_{t1\_(k + 1)} \left( {i,j} \right) = & a_{0} \left( j \right) + a_{1} \left( j \right) \cdot \left( {w(j) \cdot f_{t1\_k} \left( {i,j} \right) + \left( {1 - w(j)} \right) \cdot f_{t1\_0} \left( {i,j} \right)} \right) \\ + a_{2} \left( j \right) \cdot \left. {\frac{{{\text{d}}f}}{{{\text{d}}t}}} \right|_{{T_{1\_k} }} \left( {i,j} \right) \cdot \Updelta t_{k} + a_{3} \left( j \right) \cdot \tau_{t1\_k} \left( {i,j} \right) + a_{4} \left( j \right) \\ \cdot \left. {\frac{{{\text{d}}\tau }}{{{\text{d}}t}}} \right|_{{T_{1\_k} }} \left( {i,j} \right) \cdot \Updelta t_{k} + a_{5} \left( j \right) \cdot \sigma_{{t1_{k} }} \left( {i,j} \right) + a_{6} \left( j \right) \cdot \left. {\frac{{{\text{d}}\sigma }}{{{\text{d}}t}}} \right|_{{T_{1\_k} }} \left( {i,j} \right) \\ \cdot \Updelta t_{k} \\ \end{aligned} $$

  where $\left. {\frac{{{\text{d}}f}}{{{\text{d}}t}}} \right|_{{T_{1\_k} }} \left( {i,j} \right) = \frac{{f_{t1\_k} (i,j) - f_{t1\_(k - 1)} (i,j)}}{{t_{k} - t_{k - 1} }},\,\,\left. {\frac{{{\text{d}}\tau }}{{{\text{d}}t}}} \right|_{{T_{1\_k} }} \left( {i,j} \right) = \frac{{\tau_{t1\_k} (i,j) - \tau_{t1\_(k - 1)} (i,j)}}{{t_{k} - t_{k - 1} }}$ , and $\left. {\frac{{{\text{d}}\sigma }}{{{\text{d}}t}}} \right|_{{T_{1\_k} }} \left( {i,j} \right) = \frac{{\sigma_{t1\_k} (i,j) - \sigma_{t1\_(k - 1)} (i,j)}}{{t_{k} - t_{k - 1} }}$are derivatives of displacement, shear stress and solid stress, respectively, $\Updelta t_{k} = t_{k + 1} - t_{k} = \frac{T3 - T1}{m + n}$ is a time step, $f$ are $x$ and $y$ coordinates, $\tau$ are WSS values, $\sigma$ are solid stress values of nodal points, $j = 1,2, \ldots ,24$ is the slice number and $i$ is the index for the points on each slice. $a_{0} (j),a_{1} (j),a_{2} (j),a_{3} (j),a_{4} (j),a_{5} (j),a_{6} (j)$ and $w(j)$ are coefficients of growth functions GF1, GF2 and GF3 to be determined in such a way to obtain the best match of calculated geometries and experimental geometries at times T2 and T3. Since we use $m$ time steps to go from T1 to T2 and $n$ time steps to go from T2 to T3, we compared $f_{t1\_m}$ with experimental geometry at time T2 and $f_{t1\_ (m + n)}$ with experimental geometry at time T3. The previous formulas of growth functions are very similar with the formulas that Yang used in his paper [21].

Coefficients of the plaque volume growth functions (GF1, GF2 and GF3) $a_{0} ,a_{1} ,a_{2} ,a_{3} ,a_{4} ,a_{5} ,a_{6}$ and $w$ are calculated, independently for all 24 slices using simplex optimization method, the method which does not involve derivative calculations, developed by John Nelder and Roger Mead [13]. We minimized sum of the squared errors between calculated and real geometry at times T2 and T3 for each of 24 slices.

$$ \begin{aligned} {\text{ESS}}\left( j \right) = \mathop \sum \limits_{i = 1}^{{N_{j} }} \left( {\left( {x_{T2,i} \left( j \right) - \bar{x}_{T2,i} \left( j \right)} \right)^{2} + \left( {y_{T2,i} \left( j \right) - \bar{y}_{T2,i} \left( j \right)} \right)^{2} } \right) \\ + \mathop \sum \limits_{i = 1}^{{N_{j} }} \left( {\left( {x_{T3,i} \left( j \right) - \bar{x}_{T3,i} \left( j \right)} \right)^{2} + \left( {y_{T3,i} \left( j \right) - \bar{y}_{T3,i} \left( j \right)} \right)^{2} } \right) \\ \end{aligned} $$

where $N_{j}$ is the number of nodes for slice $j$, $x_{T2,i} \left( j \right)$, $y_{T2,i} \left( j \right)$, $x_{T3,i} \left( j \right)$, and $y_{T3,i} \left( j \right)$ are real x and y coordinates at time steps T2 and T3 for slice $j$, $\bar{x}_{T2,i} \left( j \right)$, $\bar{y}_{T2,i} \left( j \right)$, $\bar{x}_{T3,i} \left( j \right)$ and $\bar{y}_{T3,i} \left( j \right)$ are calculated x and y coordinates at time steps T2 and T3 for slice $j$.

The best results we obtained using growth function GF3 which takes into account wall shear and solid stress. Total squared error is calculated as:

$$ {\text{ESS}} = \mathop \sum \limits_{j = 1}^{24} {\text{ESS}}(j) $$

Total squared errors for all growth functions are GF1 = 36.02, GF2 = 29.98, GF3 = 26.31.

Total squared error does not give a picture of how our model is really accurate, it only serves to compare the results obtained with different growth function. Because of that, we calculated mean relative percent error:

$$ {\text{RE}}\left( j \right) = \frac{1}{{2N_{j} }}\left( {\mathop \sum \limits_{i = 1}^{{N_{j} }} \frac{{\Updelta P_{T2,i} \left( j \right)}}{{r_{T2,i} \left( j \right)}} + \mathop \sum \limits_{i = 1}^{{N_{j} }} \frac{{\Updelta P_{T3,i} \left( j \right)}}{{r_{T3,i} \left( j \right)}}} \right) \times 100 $$

$$ {\text{RE}} = \frac{1}{24}\mathop \sum \limits_{j = 1}^{24} {\text{RE}}\left( j \right) $$

where $\Updelta P_{T2,i} \left( j \right)$ and $\Updelta P_{T3,i} \left( j \right)$ are distances between real and predicted position of ith point of jth slice at times T2 and T3. $r_{T2,i} \left( j \right)$ and $r_{T3,i} \left( j \right)$ are distances between real position of ith point and center of gravity for jth slice at times T2 and T3.

Mean relative percent errors for all growth functions are GF1 = 2.7 %, GF2 = 2.62 %, GF3 = 2.51 %. It can be observed that model which uses growth function GF3 is the most accurate. Mean relative percent errors for GF3 is 2.51 %, which means that distance between predicted position of point and real position of point is in average only 2.51 % of distance between real position of point and slice center of gravity. This seems to be very good result.