Simulation of atherosclerotic plaque growth using computational biomechanics and patient-specific data

Atherosclerosis is the one of the major causes of mortality worldwide, urging the need for prevention strategies. In this work, a novel computational model is developed, which is used for simulation of plaque growth to 94 realistic 3D reconstructed coronary arteries. This model considers several factors of the atherosclerotic process even mechanical factors such as the effect of endothelial shear stress, responsible for the initiation of atherosclerosis, and biological factors such as the accumulation of low and high density lipoproteins (LDL and HDL), monocytes, macrophages, cytokines, nitric oxide and formation of foams cells or proliferation of contractile and synthetic smooth muscle cells (SMCs). The model is validated using the serial imaging of CTCA comparing the simulated geometries with the real follow-up arteries. Additionally, we examine the predictive capability of the model to identify regions prone of disease progression. The results presented good correlation between the simulated lumen area (P < 0.0001), plaque area (P < 0.0001) and plaque burden (P < 0.0001) with the realistic ones. Finally, disease progression is achieved with 80% accuracy with many of the computational results being independent predictors.


Results
We performed the plaque growth simulation at 94 coronary arteries, which resulted to 900 3 mm sub-segments. The simulations resulted to the deformed arterial geometries due to plaque growth. The simulation was performed assuming a time period same with the interscan period of each patient applying timesteps of one months. A case example of an arterial segment with increased lumen reduction and plaque increase is shown in Fig. 1. Table 1 presents the mean and Standard Deviation of the lumen and wall area and plaque burden change for the simulated and real arteries as found in the CTCA images.
Among the results, all species' distributions within the arterial wall were extracted in all time steps of the simulations, in order to further demonstrate each patient's arterial wall composition over time. The comparison of the simulated and the real follow-up arterial lumen area showed a statistical significant correlation (r = 0.608, P < 0.0001). Similar correlation is also found between the simulated wall area with the corresponding  www.nature.com/scientificreports/ real follow-up area (r = 0.494, P < 0.0001). Finally, the simulated plaque burden is also associated with the real follow-up plaque burden, presenting, however, lower correlation coefficient (r = 0.332, P < 0.0001). These correlations are also presented in Fig. 2 in the form of scatter plot. Going one step further, we hypothesize that the plaque growth model and its variables could be further used for the prediction of disease progression. For this purpose, we estimate the disease progression between the follow-up and the baseline for each 3 mm sub-segment. Initially, linear regression analysis was performed to identify correlations between the computational variables and the lumen area, plaque area and plaque burden change. We found a strong relation of ESS with disease progression, since regions of low ESS present increased lumen area reduction (P < 0.0001) and an increase of the plaque area and plaque burden at the follow-up (P < 0.0001). This is explained by the fact that in areas of high shear stress, blood velocity is increased preventing several species from entering in the arterial wall, as they drift along with blood. This mainly occurs in cells which due to their large size maintain their momentum, drifting along with blood. The arterial wall concentrations of foam cells,   Figure 3 presents the relation of some computational variables with the change of lumen area, plaque area and plaque burden. A multi-variate linear regression model was also built. Collinear variables (> 0.8 collinearity) were not included in the statistical model. Regarding the lumen area change, the included variables were the baseline plaque burden (B = 0.048, P < 0.0001), the collagen concentration (B = − 79.26, P < 0.0001) and the HDL concentration (B = 1538.00, P < 0.0001). The associated variables regarding the prediction of wall area change are the baseline lumen area (B = 0.16, P < 0.0001) and baseline wall area (B = − 0.52, P < 0.0001), the monocytes concentration (B = 1.99 × 10 -9 , P = 0.003), the macrophages concentration (B = − 1.07 × 10 -11 , P < 0.0001), the oxidized LDL concentration (B = 2441.94, P < 0.0001), the HDL concentration (B = − 3411.139, P < 0.0001) and the ESS (B = − 0.16, P = 0.014). Finally, regarding the plaque burden change, the associated variables are baseline plaque burden (B = − 0.45, P < 0.0001), the monocytes concentration (B = 2.48 × 10 -8 , P < 0.0001), the cytokines concentration (B = − 1.97, P < 0.0001), the oxidized LDL concentration (B = 11,967.47, P < 0.0001), the HDL concentration (B = − 0.11, P = 0.002), baseline lumen area (B = − 0.62, P = 0.001) and the collagen concentration (B = 195.03, P = 0.036).
Finally, the prognostic value of the plaque growth model in identifying segments of disease progression was achieved performing binary logistic regression modeling. To this purpose, we first performed a receiver-operator characteristics analysis to identify cut-off points to make binary the associated variables identified in the linear multi-variate analysis. Lumen area change, plaque area change, and plaque burden change were transformed to binary assuming 20% change. The accuracy to predict lumen reduction, plaque increase, and plaque burden increase is 83%, 80% and 77%, respectively.

Discussion
In this work, we have developed a new plaque growth model which simulates the blood flow dynamics, the species transport in the arterial wall, the oxidation of LDL, the inflammation, the formation of foam cells and finally the development and growth of plaque consisted of smooth muscle cells, collagen and foam cells. Specifically, this model introduces each one of these features in one of its three modelling levels: (i) blood flow modelling and evaluation of ESS distribution, (ii) lipoprotein (LDL, and HDL) transport within the arterial wall, inflammation modelling and plaque volume evaluation, (iii) wall thickening modelling. It is the first time that a complex model has been applied to a population with serial imaging-based assessment to validate its outcomes. More interestingly, the selected population is considered low-risk, since no invasive angiography is available for these patients. Moreover, for the first time, prediction of disease progression is achieved with 80% accuracy using computational results other than the ESS and the LDL concentration.
A significant characteristic of this model is that the endothelial dysfunction was integrated in the endothelial permeability model. Endothelial dysfunction is the culprit of atherosclerosis initiation, while some major influential factors leading to endothelial dysfunction are the ESS magnitude and the endothelial nitric oxide concentration level, where both of them are included in this model. More specifically, low ESS is experimentally found to promote the endothelial fluxes of lipoproteins, a fact that is explained due to the increase of the endothelium intercellular space.
Increased lipoproteins accumulation within the arterial wall triggers the inflammatory process initiated with the oxidation of LDL caused by free radicals. Subsequently, the presence of OxLDL within the arterial wall causes cytokine production, which is an inflammatory signalling molecule for monocyte attraction. Therefore, monocytes accumulate within the affected area and differentiate into macrophages to uptake the OxLDL. The presence of clusters with foam cells in areas containing abnormal lipoprotein concentrations was established by histological analyses, ultimately exposing type I lesions. Increased concentrations of lipid-laden smooth muscle cells were detected in advanced lesions, leading to the classification of type II lesions, namely fatty streaks. Extracellular lipid particles can be encountered in further advanced lesions, preventing the coherence of smooth muscle cells www.nature.com/scientificreports/ in the intima layer, eventually leading to type III and IV lesions, distinguished by the size of the extracellular lipid droplets. Type V lesions, distinguished by the existence of thick layers of fibrous connective tissues in the lesions, prevail due to additional advancement of the latter lesion classes, typically occurring during the fourth decade of life. Lastly, manifestation of fissure, hematoma or thrombus within the lesion, led to their referral type IV lesions. Our model considers the dynamics of LDL, foam cells and smooth muscle cells within the arterial wall and therefore it constitutes a mathematical tool for predicting type I and II lesions. A novelty of this work is that the wall thickening is enabled in the atherosclerotic areas and it is based on the total calculated plaque volume, which consists of foam cells, collagen and smooth muscle cells-the main components of the fatty streak plaques. These components that can identify both type I and II lesions, are used to evaluate the volumetric growth of the arterial wall, by considering the volumetric growth in each element of the arterial wall. This is important for development of prevention strategies. In fact, the selected population of SMARTool project consists of low-risk patients without events. For this reason, the developed plaque growth model can be used to predict the arterial progress of patients with at least 20% disease progression in a period of six years. A major novelty of this work is the use of a low-risk population with a large follow-up. This enables the examination of disease progression to patients who are not susceptible to disease progression compared to interventional studies, which include high-risk patients. Similar models were presented previously 6,9 , but only as proof-of-concept approaches or in ideal geometries. Our work, however, can be compared against other predictive studies of disease progression, which are mainly using ESS and LDL concentration accumulation as predictors 5,10,11 . In fact, in these studies, maximum accuracy of previous studies to predict disease progression is less than 65%. Our results show that in our population we have about 80% accuracy of disease progression defined as lumen reduction, plaque increase or plaque burden increase. This contributes to the conclusion that this complex computational model may describe better the pathophysiology of atherosclerosis and worth applied to even larger populations.
However, the main outcome of this work is the dynamics of this proof of concept study. Compared to ESS based predictive models, where only blood flow modeling is implemented, using a plaque growth model the pathophysiology of the disease is simulated. In this concept, it is possible to increase the complexity of the model by the inclusion of additional biological pathways. Considering large initiatives for data collection including omics data, a plaque growth model can be improved considerably adding the effect of genetic phenotype or implementing machine learning approaches for the prediction of disease progression.
The current work has however some limitations. First, the population even if it is the largest population used for the validation of such models, it is considered still small to make safe conclusions. However, the results demonstrate the trend that this model should be applied to larger populations and potentially can be used as a preventive tool. Also, the different stiffness of the plaque regions was neglected, since almost all patients presented low disease progression and small plaque regions. Finally, in this work, we didn't consider the plaque composition and the relation of the computational variables with specific plaque types. This task requires the use of intravascular imaging either ultrasound or optical coherence tomography. We aim to utilize such datasets in a future work.

Conclusions
We have developed and validated a multi-level plaque growth model using serial CTCA imaging from 94 patients. The plaque growth model simulates the main mechanisms of disease progression. The results show that the simulated and generated arterial geometries are correlated well with the real follow-up geometries. Additionally, we investigated the role of the plaque growth model to predict disease progression. Our results show that prediction of lumen area reduction and plaque area increase can be achieved with 80% accuracy.

Methods
Population, CTCA analysis and 3D reconstruction. CTCA imaging was acquired from the SMARTool clinical study, which is a multi-center EU funded (GA number: 689068) project aiming to the development of decision support systems for the management of CHD in terms of risk stratification, diagnosis and prognosis. Each patient's dataset included baseline clinical and biohumoral data and CTCA imaging at two time points. In the current analysis, we have used only the LDL and HDL concentration as boundary condition, the blood pressure applied in the equations and the interscan period of each patient to define the maximum simulation period. The clinical characteristics of our population is presented in Table 2. For the present analysis 94 patients were selected. More specifically, a population of 275 patients was created in SMARTool project, from which 12 patients did not perform follow-up CTCA and 76 were excluded due to stenting after the baseline examination or low image quality. This results at 187 patients, from which we randomly selected 50% of them (94 patients) from various clinical centers. Full explanation of the investigational nature of the study was provided to all participants and written consensus obtained. Ethical approval was provided by each participating center (National Research Council, University of Turku, University of Zurich, Fondazione Toscana Gabriele Monasterio, Warsaw National Institute of Cardiology) through the approval of the clinical study by the Ethics Committee Vast Area Northwest of Tuscany (CEAVNO), Pisa, Italy, and all subjects gave written informed consent. Our clinical study follows the declaration of Helsinki.
3D reconstruction of the coronary lumen and outer vessel wall was performed using an in-house software, which provided measurements of lumen area, plaque volume and plaque burden as previously described and validated 12,13 . Finally, baseline and follow-up CTCA scans were co-registered using landmarks such as the bifurcations and calcified objects.  15 . The fluid dynamics of blood is based on the Navier-Stokes equations for laminar, incompressible and Newtonian flows, that account for momentum and continuity conservation, respectively. A steady blood flow was assumed, since this model considers a time duration of years and therefore a steady flow profile is assumed that does not affect the plaque growth (1-3).

Multi
ρ is the blood density, U is the blood velocity, P is the pressure, τ is the shear stress, μ blood is the blood's dynamic viscosity and S M accounts for momentum sources. Plasma is considered as an incompressible Newtonian fluid that presents a density value of 1000 [kg/m 3 ] 14 and a dynamic viscosity value of 0.001 [Pa s] 15 . Plasma flow in the arterial wall is governed by the modified where, τ = −µ blood ∇U + (∇U) T − 2 3 δ∇ · U . www.nature.com/scientificreports/ Navier-Stokes equations laminar, incompressible and Newtonian flows inside a porous media, since the arterial wall can be characterized as such, with a porosity coefficient 0.96 16 . These equations derive from the Navier-Stokes and continuity equations, in which the reduced flow space is considered using the area porosity tensor K while the momentum losses are considered in a source term S M . In a material of porosity γ, the area porosity tensor is calculated by (4). The momentum losses S M are formulated using the permeability K perm and loss coefficient K loss , but in our case of a low plasma velocity within the arterial wall, the term including the loss coefficient K loss is neglected as it includes a second degree of the velocity that is tiddley (8). Using the porosity γ and the area porosity tensor K, the equations for the conservation of mass (5) and momentum (6) account for the decreased flow space and the momentum losses 17,18 .
μ plasma is the plasma viscosity, K perm is the Darcian permeability, γ is the volume porosity and K is the area porosity tensor, which is a symmetric second rank tensor. Momentum losses, S M , account for the losses due to the permeability decrease and the velocity direction changes 17,18 .
Dynamics of the atherosclerotic process. In lumen, blood flows along with LDL, HDL and monocytes. The blood concentration of LDL, HDL and monocytes are patient specific and therefore are used as boundary conditions. However, after their infiltration in the arterial wall, their concentrations within the arterial wall are calculated using the convection-diffusion-reaction Eq. (9), which is also modified to account for flow in porous media as presented previously 18,19 .
The first term accounts for time-dependent concentration changes, while the second is the advection term and accounts for concentration changes due to velocity drift. The third one is the diffusion term, where D is the diffusivity. The last term accounts for concentration sources or sinks.
LDL macromolecules react with the free radicals in the arterial wall, forming oxidized LDL and therefore their concentration reduces 20 . However, HDL macromolecules, that infiltrate to the arterial wall, also react with the free radicals and reduce the free radical concentration, leading to a decrease of the LDL oxidation rate. An HDL-dependent oxidation rate of LDL was implemented as it was presented by Sakellarios et al. 6 which was based on the experimental results of Esterbauer et al. 20 , and it is referred as HDL protection .
In the LDL sink term, r LDL (= 1.4 × 10-4 s −120 ) is the LDL oxidation rate, in the absence of HDL macromolecules, while HDL protection (= − 3 × 10 −5 c LDL 2 + 5 × 10 -4 c LDL + 1.0056 6 ) corresponds to the relative decrease of the LDL oxidation rate due to the presence of HDL macromolecule.
The presence of oxidized LDL (OxLDL) within the arterial wall sets an inflammatory signalling to monocytes, which are accumulated within the arterial wall and differentiate into macrophages to uptake OxLDL, forming foam cells. k 2 is the OxLDL uptake rate from monocytes, (= 12 × 10 −19 m 3 cells −1 s −121 ) and depends on both OxLDL and macrophage concentration.
Each equation describing a cells' dynamics has a zero term for advection, due to the cell's large size in respect to the porous dimension. In most of them the diffusion term is also neglected for the same reason. However, the diffusion terms of monocyte and macrophage dynamics are considered, as we assume that these immune cells can move according to diffusion forces as previously presented by Cilla et al. 9 .
Cytokines are produced due to the presence of macrophages and OxLDL within the arterial wall and cause the differentiation of contractile smooth muscle cells (SMCs) into synthetic SMCs.
The cytokine degradation rate is d c (= 2.3148 × 10 -5 s −123 ) and its production rate is d r that depends on the OxLDL and macrophage concentration. The contractile SMCs differentiation rate into synthetic SMCs has been found experimentally to be equal to (1 + exp(− S r c cytokines /c cytokines,max )), where S r = 4.16 × 10 -8 s −124 , showing its explicit dependence to the cytokine concentration.
Synthetic SMCs secrete collagen as connective tissue, further increasing plaque volume.
Endothelial flow rates. The infiltration of species and cells through the endothelium layer can be performed by three different pathways. These are the vesicular transcytosis pathway and the pathways through leaky and normal junctions. The first one accounts for intercellular transport, which is only active in the presence of endothelial cell receptors 27 . The other two, account for transport between the endothelium cells, which are regulated only by mechanical and diffusion factors, such as pressure and concentration differences across the endothelium. The volume fluxes of vesicular transcytosis are negligible in relation to the other two, which allows the application of the Kedem-Katchalsky equations that calculate fluxes based on pressure and concentration differences across biological membranes 28 . The endothelial monocyte fluxes initiate after inflammatory signalling and are defined by an experimental equation using both the OxLDL concentration and the blood monocyte concentration, but also the endothelial ESS. In particular, high ESS retain monocyte from infiltration.
Nitric oxide concentration is considered implicitly in the LDL and HDL diffusive permeability. It is a chemical product of the endothelial nitric oxide synthases (eNOS) and its concentration depends on the partial pressure of oxygen and the eNOS concentration 29 . Following, eNOS concentration dependence of the applied endothelial WSS has been proven experimentally 30 .
Equation (21) accounts for plasma velocity through the endothelium, while Eqs. (22)(23)(24) account for LDL, HDL and monocyte fluxes, respectively. The second term of the plasma velocity equation, which corresponds to the decrease of velocity due to osmotic pressure differences, is neglected. L p is the hydraulic conductivity of the endothelium and it is a function of the endothelial ESS developed from blood flow. DP LDL and DP HDL are the LDL and HDL diffusive permeability respectively, and are a function of the endothelial nitric oxide concentration (13) ∂ ∂t γc monocytes = ∇ · D m K · ∇c monocytes − γm d c monocytes − γ d m c monocytes .   Endothelial NO concentration depends on the partial pressure of oxygen and the endothelial nitric oxide synthase concentration (eNOS), which is a function of the mean applied endothelial shear stresses 30 .
Wall thickening. To simulate the arterial wall thickening, the arterial wall is treated as a solid domain, which expands due to the volumetric strain that is caused by the generation of plaque. More specifically, the volumetric strain is calculated based on the simulated plaque volume, which enables a strain-based arterial wall thickening. Volumetric strain is defined as the change in volume divided by the original volume 32 . Although, the most used stress-strain relationship of the arterial wall is a 6-parameter Mooney-Rivlin material model 33,34 , we implemented a linear elastic material model, since the volumetric strains that are caused due to plaque growth are relatively small. According to that the arterial wall presents a young modulus of 1.06 MPa and a Poisson ratio of 0.45 7 .
Due to the low plaque progression rate of our population, we assumed that a one-way interaction between the CFD analysis of the atherosclerotic process and the structural analysis of the wall thickening is adequate. Therefore, in our analysis, the volume change of each finite element is equal to the total volume of the foam cells, collagen and synthetic muscle cells, because their initial concentration inside the arterial wall was considered zero.
where, in (11), c Foam cells , c Synthetic SMC and c collagen are the concentrations of foam cells (cells/m 3 ), synthetic SMCs (cells/m 3 ) and collagen (g/m 3 ) respectively, while V Foam cells , V Synthetic SMC and v collagen are the cellular volume of foam cells and the cellular volume of synthetic SMCs and the specific volume of collagen respectively. Sequentially, in (30), V is the volume of a cubic finite element and ε is the resulted strain after a change in volume ∂V. In our model, we consider a constrained arterial wall thickening that is available only to the direction of the centerline and therefore the directional volumetric strain ε directional can result from the function: Numerical implementation. The simulations were conducted using the ANSYS software (ANSYS, Canonsburg, PA) which integrate among other finite element and finite volume analysis solvers. Specifically, fluid dynamics are performed using the ANSYS CFX software, while the wall thickening is performed in ANSYS mechanical module.
Both lumen and the arterial wall domain of each patient were meshed using tetrahedral elements of 0.15 mm edges. Several methods were used to improve the mesh quality, such as the patch-independent technique and the inflation. Specifically, these parameters were evaluated after monitoring of the element quality statistics. Moreover, a sensitivity analysis of a normal/typical patient arterial segment was performed based on the tetrahedral element size, to further validate the mesh quality (Fig. 4).

Boundary conditions. Fluid domain of lumen.
To simulate the blood flow in lumen, inlet is conditioned by the patient blood velocity, outlet is conditioned by the patient blood pressure and the endothelium layer is considered as an impenetrable and no-slip wall, since endothelial fluxes are too small to affect the blood flow.
Fluid domain of arterial wall. Plasma flow in the arterial wall is constrained with a zero-velocity to the vertical sides of the arterial wall, a standard pressure to the adventitia layer and a Kedem-Katchalsky derived velocity at endothelium. The arterial wall vertical sides are impenetrable to any substance and species, while at the adventitia layer all substances can pass through. The adventitia layer is also constrained with patient-specific LDL and HDL concentration values.
Scientific Reports | (2020) 10:17409 | https://doi.org/10.1038/s41598-020-74583-y www.nature.com/scientificreports/ Validation. Plaque growth simulations were performed at the baseline reconstructed arteries and we compared the simulated results with the realistic follow-up examination. The validation approach is based on the rationale that the plaque growth model enables the arterial wall deformation in a way that reaches the follow-up geometry as this has been assessed by the CTCA. The validation approach requires the extraction of 0.5 mmdistanced cross-sections of both the baseline and follow-up coronary arteries as shown in Fig. 5. Considering that CTCA has also a mean slice thickness of 0.5 mm, we implemented an approach, proposed by Stone et al. 35 , which is considered good approach for the comparison of computational results with morphological findings from the imaging data or 3D reconstructed arteries. In particular, to eliminate possible registration errors, the rationale is to generate 3 mm segments by combining six sequential 0.5 mm-distanced cross-sections for all geometries (baseline, simulated and follow-up) (Fig. 5).
Parameters of the plaque growth model. The plaque growth model simulates the major mechanisms of atherosclerosis by employing many differential equations. The equations require some parameters, which in majority are found in literature in experimental studies. All the parameters of the plaque growth model are shown in Table 3.

Statistical analysis.
All continuous variables are presented as mean ± standard deviation. Wilcoxon, Fisher's exact and chi-squared tests were used for the comparison of characteristics within group and between groups. Initially, linear regression analysis was performed between all computational variables and the morphological characteristics (lumen area, plaque area, plaque burden and their changes between the baseline and follow-up examination). In the next step, a step-wise multivariate regression model was employed, in which only the associated variables (P < 0.1) were entered.