Effects of the Haemodynamic Stimulus on the Location of Carotid Plaques Based on a Patient-Specific Mechanobiological Plaque Atheroma Formation Model

In this work, we propose a mechanobiological atheroma growth model modulated by a new haemodynamic stimulus. To test this model, we analyse the development of atheroma plaques in patient-specific bifurcations of carotid arteries for a total time of 30 years. In particular, eight geometries (left or right carotid arteries) were segmented from clinical images and compared with the solutions obtained computationally to validate the model. The influence of some haemodynamical stimuli on the location and size of plaques is also studied. Plaques predicted by the mechanobiological models using the time average wall shear stress (TAWSS), the oscillatory shear index (OSI) and a new index proposed in this work are compared. The new index predicts the shape index of the endothelial cells as a combination of TAWSS and OSI values and was fitted using data from the literature. The mechanobiological model represents an evolution of the one previously proposed by the authors. This model uses Navier-Stokes equations to simulate blood flow along the lumen in the transient mode. It also employs Darcy's law and Kedem-Katchalsky equations for plasma and substance flow across the endothelium using the three-pore model. The mass balances of all the substances that have been considered in the model are implemented by convection-diffusion-reaction equations, and finally the growth of the plaques has been computed. The results show that by using the new mechanical stimulus proposed in this study, prediction of plaques is, in most cases, better than only using TAWSS or OSI with a minimal and maximal errors on stenosis ratio of 2.77 and 32.89 %, respectively. However, there are a few geometries in which haemodynamics cannot predict the location of plaques, and other biological or genetic factors would be more relevant than haemodynamics. In particular, the model predicts correctly eleven of the fourteen plaques presented in all the geometries considered. Additionally, a healthy geometry has been computed to check that plaque is not developed with the model in this case.

In this work, we propose a mechanobiological atheroma growth model modulated by a new haemodynamic stimulus. To test this model, we analyse the development of atheroma plaques in patient-specific bifurcations of carotid arteries for a total time of 30 years. In particular, eight geometries (left or right carotid arteries) were segmented from clinical images and compared with the solutions obtained computationally to validate the model. The influence of some haemodynamical stimuli on the location and size of plaques is also studied. Plaques predicted by the mechanobiological models using the time average wall shear stress (TAWSS), the oscillatory shear index (OSI) and a new index proposed in this work are compared. The new index predicts the shape index of the endothelial cells as a combination of TAWSS and OSI values and was fitted using data from the literature. The mechanobiological model represents an evolution of the one previously proposed by the authors. This model uses Navier-Stokes equations to simulate blood flow along the lumen in the transient mode. It also employs Darcy's law and Kedem-Katchalsky equations for plasma and substance flow across the endothelium using the three-pore model. The mass balances of all the substances that have been considered in the model are implemented by convection-diffusion-reaction equations, and finally the growth of the plaques has been computed. The results show that by using the new mechanical stimulus proposed in this study, prediction of plaques is, in most cases, better than only using TAWSS or OSI with a minimal and maximal errors on stenosis ratio of 2.77 and 32.89 %, respectively. However, there are a few geometries in which haemodynamics cannot predict the location of plaques, and other biological or genetic factors would be more relevant than haemodynamics. In particular, the model predicts correctly eleven of the fourteen plaques presented in all the geometries considered. Additionally, a healthy geometry has been computed to check that plaque is not developed with the model in this case.

INTRODUCTION
Atherosclerosis is a disease that causes the formation of atheroma plaques in arterial walls. The effect of atheroma plaques is that the thickness of the arterial wall increases-due to an accumulation of some substances such as low density lipoproteins (LDL) and foam cells (FC) in it-and, therefore, the lumen area decreases and blood cannot flow properly. It can derive in several events, such as heart attacks, ischaemias or strokes, and currently it is one of the main causes of mortality in developed countries (Gaziano and Gaziano, 2012). Although this pathology has been widely studied, it has not yet been completely understood. Therefore, it is relevant to study the process of formation of atheroma plaques and to foresee the locations that are susceptible to the emergence of plaques in arteries.
It has been accepted that some mechanical stimuli can cause shape changes of endothelial cells (Dai et al., 2004), and depending on it, they can induce LDL transport into the arterial wall through the endothelium, initiating the growth of atheroma plaques in the vessel. These mechanical stimuli can depend on several factors such as cyclic stretches, cardiac cycle, geometry of the arteries and oscillatory shear stress (Ohayon et al., 2011).
One of these mechanical stimuli is the wall shear stress (WSS) caused by blood flow in the endothelium. It is an index that has been widely used to predict the location of plaques; however, it has the limitation of being calculated for stationary blood flow and does not consider the cardiac cycle. It is well known that areas with physiological WSS promote endothelial cells to have an elongated shape, so pores between them are small and limit the flow of substances across the endothelium. In contrast, for areas with very low WSS, endothelial cells are more circular, so pores are larger and allow flow of substances between them, resulting in plaque emergence. The threshold of WSS below which plaques grow depends on the considered artery (Olgac et al., 2008;Filipovic et al., 2013). In the case of carotid arteries, areas with WSS lower than 2 Pa could be considered atheroprones, while areas of higher WSS are atheroprotectives (Zhao et al., 2002;Younis et al., 2004;Filipovic et al., 2013).
To avoid the limitation of not considering transient blood flow, other studies use the time averaged wall shear stress (TAWSS) instead of WSS to take into account the cardiac cycle and to improve the accuracy of the prediction (Sáez et al., 2015;Alimohammadi et al., 2017). Another index that has also been used in some studies is the oscillatory shear index (OSI). There is some evidence about the influence of this mechanical stimulus on cell shape, and therefore in the emergence of plaques, being areas of high OSI susceptible to developing plaques (Alimohammadi et al., 2017). This index also considers the complete cardiac cycle. However, most of these studies only take into account TAWSS or OSI to predict the location of plaques and do not consider the inflammatory process to reproduce the growth of plaques.
There are other indices that have been investigated recently such as Cross-flow index (CFI) (Arshad et al., 2020), Transverse Wall Shear Stress (transWSS) (Peiffer et al., 2013) and Topological Shear Variation Index (TSVI) (Morbiducci et al., 2020), but their implementation into the model was not possible due to that there are not enough experimental data correlating SI of endothelial cells with them (Morbiducci et al., 2020).
Finally, some studies use patient-specific geometries and calculate blood flow in the transient mode with growth of plaques, but they do not consider all the substances that take part in the disease progression and do not add volume to the final plaque (Filipovic et al., 2013;Díaz-Zuccarini et al., 2014;Alimohammadi et al., 2017).
In this study we use patient-specific geometries with different degrees of atheroma plaques in carotid arteries to fictitiously reconstruct healthy arteries and to computationally reproduce the inflammatory process of the emergence of plaques in them. We use transient blood flow, taking into account the cardiac cycle, and analyse plaque growth under three different hypotheses considering distinct mechanical stimuli: TAWSS, OSI, and a combination of them that we propose as a new stimulus.
The aim of this study is to analyse the predictability of the different mechanical stimuli to estimate the emergence of plaques, improving a previous model developed by the authors under an axisymmetric hypothesis (Cilla et al., 2014). Finally, we compare the plaques predicted using the computational mechanobiological model with the real plaques of patients in clinical images in order to determine which haemodynamical stimulus can better predict the location and size of plaques.

MATERIALS AND METHODS
The mechanobiological model proposed by Cilla et al. (2014) was improved by the introduction of a new mechanical stimulus and some terms were simplified in order to improve the numerical convergence, see section 2.5. The endothelium was modelled as a thin layer of endothelial cells, which change their shape as a function of the haemodynamical stimulus -TAWSS, OSI and a combination of them-by becoming rounder or elongated according to the stimulus and thus allowing more or less substance transport along the endothelium. The arterial wall was modelled as a single layer (intima-media) with a permeable membrane (endothelium).

Patient-Specific Geometries
Clinical images of four different male patients with atherosclerosis and one healthy volunteer were segmented using the software Materialise Mimics (Materialise N. V., Leuven, Belgium) to obtain eight different patient-specific geometries of carotid artery bifurcations, including common, internal and external carotid arteries (CCA, ICA and ECA, respectively). The clinical images were provided by the Hospital Clinico Universitario in Zaragoza, Spain, according to ethics guidelines of the hospital. One geometry corresponds to a healthy volunteer without pathology, and the others correspond to four patients with developed atheroma disease. Carotid images are shown in Figure 1, with their respective plaques indicated by arrows. The images in Figure 1 correspond to the real geometries of the patients, without making changes in them, to show the differences between the real carotids and the computed ones, that can be observed in Figure 7. These differences are due to the simplifications that are necessary to compute the model, e.g., the small side branches of the carotids were eliminated in the computed geometry to simplify the model and it could have some influence on the blood flow distribution.
A thresholding segmentation technique was used to reconstruct the lumen of the vessel from the clinical images. Once the thresholding was done, the plaques of the carotids were located, and plaques were digitally removed to obtain geometries that we considered healthy arteries previous to the development of the pathology. The geometry corresponding to the healthy volunteer was not modified. Finally, the arterial wall was extruded with the software Rhinoceros (Robert McNeel & Associates, Seattle, WA, United States) from the lumen to obtain a 3D geometry with variable thickness, imposing a thickness of 0.7 mm for the CCA and 0.53 mm to the ICA and ECA. Finally, the thickness in the area close to the bifurcation was progressively reduced when advancing from the CCA to de ICA and ECA according to their respective thicknesses (Sommer et al., 2010).
The different geometries were coded from "A" to "H." The patient-specific geometry "A" was used to estimate the parameters not found in the literature to computationally reproduce the location and size of the real plaque. Once done these estimations, these parameters were used for the rest of geometries. The healthy volunteer, named "H, " was used to demonstrate that the mechanobiology model can also predict a healthy case without relevant growth of plaques. The total time of the numerical simulations implemented was 30 years (Insull, 2009) and patients were supposed to have a high level of hypercholesterolemia with an LDL concentration in blood of 6.98 mol m 3 , equivalent to 270 mg dL (Goldstein and Brown, 1977).

Numerical Methods
The geometries were meshed using triangular elements. Mesh sensitivity analysis was performed for both the lumen and the arterial wall, including the number of boundary layers, to determine optimal meshes to finally compute the whole process. The final mesh had a total of 850,000 elements for the lumen, with two boundary layers near the endothelium, and 550,000 elements for the arterial wall, with boundary layers near the endothelium and the adventitia to ensure the correct calculation of fluxes across the arterial wall. The software COMSOL Multiphysics (COMSOL AB, Burlington, MA, USA) was used to computationally solve the model following four consecutive stationary and transient steps. A first transient step was used to simulate the blood flow along three cardiac cycles. Then, a second step under stationary hypothesis was performed to solve the plasma flow across the endothelium. Afterwards a third step was computed in transient mode to calculate the concentrations of all the substances in all the arterial wall during 30 years, and finally, a fourth and last stationary step was made to compute the growth of the plaques from the concentration of all the substances at 30 years. In Figure 2 there is a scheme with the followed workflow.
A direct linear solver (PARDISO) was used to solve the transient blood flow along the lumen. Another direct linear solver (MUMPS) was employed to compute the plasma flow through the endothelium, the inflammatory process of all the substances (in an iterative way, using different segregated steps for groups of substances) and finally the growth of the plaques.
The development of the mathematical model is presented in the following sections, separating the equations referring to the blood flow, plasma flow, inflammatory process and growth of plaques.

Blood Flow Model
According to Caro et al. (1978) and Perktold et al. (1991), blood was modelled as a Newtonian and incompressible fluid because the lumen diameter of the considered arteries is higher than 0.5 mm. Additionally, blood flow was considered laminar due to the Reynolds number in carotid arteries under physiological conditions. Blood is basically composed of a liquid component called plasma, but it also contains solid particles. Nevertheless, these particles are very small in comparison to the lumen diameter. Therefore, blood was considered a homogeneous fluid (Malvè et al., 2014). Blood flow in the lumen is governed by Navier-Stokes and continuity equations: where subscripts b and l refer to blood and lumen, respectively, so parameters ρ b and µ b are the density and dynamic viscosity of blood, respectively, while u l and P l are the velocity and pressure of blood flow in the lumen, respectively. Finally, F l corresponds to internal forces of the fluid, which are negligible in comparison with the friction between blood flow and the arterial wall. All parameters necessary to calculate blood flow along the lumen are shown in Table 1.
Blood flow was modelled in the transient mode. Therefore, an analysis of the number of cardiac cycles necessary to model blood flow was performed, obtaining that a total number of three cycles is sufficient to completely develop blood flow and establishing the validity of the results of flow obtained for the third cardiac cycle. At the inlet of the lumen, transient mass flow was imposed, and transient pressures were imposed at the two outlets of the lumen. Transient flow and pressures were imposed following their respective shapes along a cardiac cycle obtained from Malvè et al. (2014). Additionally, Murray's law was applied in all the geometries to establish the correct division of blood flow at the bifurcations, with an average pressure at the outlet of the ICA of 100 mmHg. In Figure 3, blood flow at the inlet of the CCA and pressures at the outlets of ICA and ECA imposed at geometry A can be seen. Finally, a non-slip condition was imposed at the endothelium.

Plasma Flow Across the Endothelium
The arterial wall is permeable; therefore, some elements contained in the blood flow can cross the arterial solid wall. In particular, there is a plasma flow through the endothelium from the lumen that can be modelled with Darcy's Law: where u w and k w are the velocity of the plasma on the arterial wall and its permeability, respectively. ∇p w is the pressure gradient in the arterial wall, and finally, µ p is the dynamic viscosity of plasma. Furthermore, continuity of plasma flow has to be accomplished: where ρ p and ǫ w are the density of plasma and the porosity of the arterial wall, respectively. The last term, denoted as J v , is the plasma flow through the endothelium. It can be calculated with Kedem-Katchalsky equations, considering three different types of pores in the endothelium by which plasma flow is allowed: normal junctions, leaky junctions and vesicular pathways. However, plasma flow through vesicular pathways is very small related to the other two, so it is negligible (Olgac et al., 2008). Total plasma flow throughout the endothelium (J v ) can be calculated as: where Jv nj , Jv lj , and Jv v are plasma flows across normal junctions, leaky junctions and vesicular pathways, respectively. Finally, Lp nj and Lp lj are the hydraulic conductivities of normal and leaky junctions, respectively. The value of Lp nj depends on the thickness of the arterial wall (and, therefore, of the artery that we consider) and of the intraluminal pressure (Tedgui and Lever, 1984). P is the pressure drop in the endothelium, which depends on the intraluminal pressure and takes values of 18 and 28 mmHg for an intraluminal pressure of 70 and 180 mmHg, respectively, Tedgui and Lever (1984). Finally, lj is the fraction of leaky junctions and is defined as the ratio of the number of leaky cells and the total number of cells (Weinbaum et al., 1985;Huang et al., 1994;Huang and Tarbell, 1997), and in our model, it depends on the haemodynamical stimulus and the number of mitotic cells, see section 2.6.
Following Weinbaum et al. (1985) and Yuan et al. (1991), hydraulic conductivity of leaky junctions can be calculated as: where Ap S is the fraction of total area occupied by leaky junctions, Lp slj is the hydraulic conductivity of a single leaky junction, R cell the radius of endothelial cells, and w l the half-width of a leaky junction (Weinbaum et al., 1985;Yuan et al., 1991;Huang et al., 1994). For more details about the derivation of (6), see Appendix A.
The total number of mitotic cells in the endothelium depends on the haemodynamical stimulus that we are using. In experimental studies, the number of mitotic cells (MC) was determined in areas of the known shape index (SI) (Chien, 2003), so it was developed in the next experimental correlation with a unit area of 0.64 mm 2 : In addition, based on experimental studies, the next correlation between the number of leaky cells (LC) and mitotic cells for the unit area of 0.64mm 2 is defined (Lin et al., 1989;Olgac et al., 2008): lj is defined as the ratio between the area of leaky cells and the area of all the cells, and it can be calculated as: taking as A unit the unit area considered in all the anterior experimental correlations of 0.64 mm 2 . By using all the experimental correlations, it is now possible to obtain lj and, therefore, A p S . Consequently, knowing the SI that it is computed by Equations [(32)-(35)] for our model, we can compute lj and the hydraulic conductivity Lp lj .
Finally, the hydraulic conductivity of a unique leaky junction, Lp slj , is defined following Olgac et al. (2008): where µ p is the dynamic viscosity of the plasma, and w l and l lj the width and the length of a leaky junction, respectively. Therefore, with this derivation, plasma flow through the endothelium is completely determined. The parameters necessary for the plasma flow calculation are depicted in Table 2.
The normal velocity of plasma flow through the endothelium is J v , which has already been defined in section 2.4. Additionally, the pressure at adventitia defined in Olgac et al. (2008) is also prescribed (17.5 mmHg).

Inflammatory Process of the Arterial Wall
Once plasma flow across the endothelium has been modelled, we can compute the inflammatory process that takes place on the arterial wall. There are many substances involved in this process, among which we consider LDL, oxidized LDL (LDLox), monocytes (m), macrophages (M), cytokines (C), contractile and synthetic smooth muscle cells (CSMC and SSMC), foam cells (FC), and collagen (G).
The behaviour of cells and substances on the arterial wall obeys convection-diffusion-reaction equations of the form: where X i is the concentration of the considered substance and D X i its diffusion coefficient on the arterial wall. K lag is the solute lag coefficient. The first term of the equation corresponds to temporal variations of the cells or the substances on the arterial wall, while second and third terms are, respectively, diffusion and convection of the cells or substances. Finally, the reaction term represents the interaction between cells and/or substances (chemotaxis, proliferation, differentiation, apoptosis, degradation or generation), and they are different for each of the considered cells or substances. All the parameters for the inflammatory process are collected in Table 3.
On the other hand, flow of substances across the arterial wall can be defined as: Initial concentrations of all the substances at the artery wall are null, except for the case of CSMCs, by which all of the arterial wall is composed at the beginning of the inflammatory process. In addition, we suppose a hypercholesterolemia level for all the patients, and their concentrations of LDL and monocytes at the lumen are C LDL,l and C m,l , respectively. The specification of these equations for each of the substances and cells of the process is as follows.

Evolution of LDL Concentration
Due to their small size, LDL molecules suffer convection due to plasma flow across the arterial wall. They additionally have diffusion. Once LDL molecules are on the arterial wall, they are oxidized, so their reaction term is: where d LDL is the degradation ratio of LDL on the arterial wall, and C LDL,w its concentration at each time. LDL flow through the endothelium can be calculated with the Kedem-Katchalsky equation (Olgac et al., 2008).
where C LDL,l is the LDL concentration at the lumen, LDL dep the quantity of LDL molecules that are deposited into the arterial wall and P app the coefficient of apparent permeability of the arterial wall, which is composed of the permeability of normal junctions, leaky junctions and vesicular pathways (P app,nj , P app,lj and P app,v , respectively): Molecule transport through endothelium occurs in different ways depending on the size of the particles. For molecules with a size lower than 2 nm, transport is allowed through all the possible ways, but for greater molecules (such as LDL, whose size is approximately 11 nm), transport across normal junctions is not allowed, so for this case, molecular transport through the endothelium only occurs by leaky junctions and vesicular pathways. According to Olgac et al. (2008), molecular transport of LDL through vesicular pathways is 0.1 of the flux through leaky junctions.
The apparent permeability of leaky junctions can be defined as: where P lj , Z lj and σ f ,lj are, respectively, diffusive permeability of leaky junctions, a factor of reduction of the concentration gradient of LDL at the entrance of flow and the solvent-drag coefficient of leaky junctions. Therefore, LDL flux across the endothelium can be written as: Diffusive permeability of leaky junctions is defined as: where χ is the difference between the total area of endothelial cells and the area of cells separated by leaky junctions, where LDL flux is allowed: where α lj is the ratio between the radius of an LDL molecule (a m ) and the half-width of a leaky junction (w l ): Finally, P slj is the permeability of a leaky junction that can be computed using the equations of Appendix B.
In addition, we have to impose as an additional boundary condition the LDL concentration at adventitia (C LDL,adv ) which is obtained from experimental data of Meyer et al. (1996) to comply with the experimental LDL distribution of LDL across the arterial wall.

Evolution of Oxidized LDL Concentration
We considered that once LDL becomes oxidized, it does not experiment convection, but it has the diffusion term. Their reaction terms are due to two factors: the first one refers to LDL that becomes oxidized in the arterial wall, and the second one refers to oxidized LDL that is absorbed by macrophages: Frontiers in Bioengineering and Biotechnology | www.frontiersin.org where C LDL ox,w and C M,w are oxidized LDL and macrophage concentration at each point of the arterial wall, respectively. LDL ox,r is the ratio of the quantity of oxidized LDL that a single macrophage absorbs.

Evolution of Monocyte Concentration
Monocytes are cells, so they do not have convection. The first reaction term of the equation corresponds to the monocytes that disappear because of their differentiation into macrophages. The second term is due to apoptosis of monocytes: where C m,w is monocyte concentration on the arterial wall, d m a parameter that represents the rate of monocytes that differentiate into macrophages, and m d monocyte rate of death.

Evolution of Macrophage Concentration
Similar to monocytes, macrophages are cells, so they do not have convection. Their reaction terms are: where the first term is the differentiation of monocytes into macrophages and the second one their apoptosis.
LDL ox,r is the constant rate of oxidized LDL taken up by macrophages, and n FC is the maximum amount of oxidized LDL that a single macrophage has to ingest to turn into a foam cell. To obtain this value, it was considered that these cells are capable of ingesting particles up to 1.44 times their radius (Cannon and Swanson, 1992), taking into account that the density and molecular weight of LDL are 1,063 kg m 3 (Ivanova et al., 2017) and 386.65 g mol (Guarino et al., 2006), respectively.

Evolution of Cytokine Concentration
Cytokines are proteins, so we did not consider convection for them. In addition, they are surrounded by macrophages, so their diffusion can be considered negligible (Cilla et al., 2014). Cytokine reaction terms are due to their degradation and production: where C c,w is cytokine concentration on the arterial wall. C r is the ratio of cytokine production due to the presence of oxidized LDL and macrophages on the arterial wall, and d c the cytokine degradation rate.

Evolution of Contractile Smooth Muscle Cell Concentration (CSMC)
Due to the size of CSMCs, they have neither convection nor diffusion. At the beginning of the inflammatory process, all the muscle cells on the arterial wall are of a contractile phenotype, but the presence of cytokines on the arterial wall make them change into a synthetic phenotype, so the reactive term of CSMCs is expressed as follows.
f C CSMC,w (C CSMC,w , C c,w ) = −C CSMC,w · S r · C c,w k c · C th c,w + C c,w (26) C CSMC,w is CSMC concentration at the arterial wall. S r is the CSMC differentiation rate due to the presence of cytokines on the arterial wall, and finally, C th c,w is the maximum cytokine concentration allowed at the arterial wall.

Evolution of Synthetic Smooth Muscle Cell Concentration (SSMC)
Analogous to CSMCs, SSMCs have neither convection nor diffusion. Their reaction terms in the arterial wall are due to differentiation of CSMCs into SSMCs, proliferation and apoptosis of SSMCs: where C SSMC,w and C th SSMC,w are SSMC concentration and their maximum allowed concentration at the arterial wall, respectively, p ss the SSMC proliferation rate and r Apop the SSMC apoptosis rate.

Evolution of Foam Cell Concentration (FC)
Foam cells have neither convection nor diffusion given that they are large cells. The reaction term is due to apoptosis of macrophages into foam cells and can be written as: All parameters in Equation (28) have already been defined.

Evolution of Collagen Fibres
Finally, collagen fibres are composed of many molecules, so they cannot move between arterial wall pores. Therefore, collagen fibres have neither convection nor diffusion. Reaction terms of collagen fibres are due to collagen segregation due to SSMC presence on the arterial wall and by collagen degradation.
where G r and d G are the collagen secretion and degradation rate, respectively, and C G,w its concentration on the arterial wall.

Haemodynamical Stimuli to Initiate the Inflammatory Process
Three different mechanical stimuli were analysed in this study as potential triggers for the inflammatory process and predictors to foresee the position and growth of atheroma plaques. The first one is TAWSS, which is defined as: where T is the period of a cardiac cycle and |τ (t)| the magnitude of WSS dependent on time, with WSS defined as: where τ x , τ y and τ z are components of the tangential stress vector appearing in the lumen-wall interface of the model. The shape index of endothelial cells directly depends on TAWSS, being proximal to 1 in the case of low TAWSS, meaning that endothelial cells are almost circular. To determine the behaviour of endothelial cells with TAWSS, we propose a numerical correlation based on the experimental results of Levesque et al. (1986). This correlation is shown in Figure 4. The endothelial shape index (SI) is: It is well accepted that areas of low TAWSS are atheroprones. In particular, for carotid arteries, areas below 2 Pa are susceptible to the emergence of atheroma plaques (Zhao et al., 2002;Younis et al., 2004;Filipovic et al., 2013). The values of the parameters k 1 , k 2 , k 3 , and k 4 are shown in Table 4.
The second mechanical stimulus that we considered is OSI: SI can be considered directly dependent on OSI; therefore, to determine this behaviour, we propose the next correlation obtained from the experimental data of Levesque et al. (1986). The graphical correlation is shown in Figure 4.
Areas with high OSI are areas in which atheroma plaques are more likely to appear. To estimate the OSI threshold, we first calculated the value of SI corresponding to a value of TAWSS of 2 Pa [Equation (32)], assuming that for this value of SI, the LDL molecules can pass through the endothelium. Therefore, by replacing this SI value in Equation (34), we can estimate that atheroma plaques will grow in areas of OSI higher than 0.1910. The values of the parameters k 5 , k 6 , and k 7 are shown in Table 4. Finally, we proposed a new index to calculate the growth of plaques as a combination of TAWSS and OSI to take into account the effect of both stimuli. For that, we used pseudo-experimental data from Sáez et al. (2015) to approximate the variable SI as a function of TAWSS and OSI, obtaining: SI = k 8 · e k 9 ·OSI + k 10 · e k 11 ·(TAWSS) 2 , using thresholds obtained before for TAWSS and OSI to determine areas of plaque growth with this new index. The approximation surface is shown in Figure 5. All the adjustment parameters used in the analysis, namely, k 8 , k 9 , k 10 and k 11 , are shown in Table 4. Monocytes flow from the lumen into the arterial wall and through the endothelium, which also depends on the haemodynamical stimulus that we consider. For TAWSS, it can be modelled with the Kedem-Katchalsky equation as (Malek and Alper, 1999;Gijsen et al., 2008): where m r is monocyte recruitment from the lumen to the endothelium. TAWSS was modelled as a sigmoid function with maximal and minimal values equal to 2 and 0 Pa, respectively, to allow LDL flux across the endothelium. To completely define the sigmoid, an average value called TAWSS 0 of 1 Pa is necessary.
For the case of using OSI as the mechanical stimulus, another equation was developed in terms of the maximal and minimal fluxes of monocytes obtained with the TAWSS equation: J s,m (OSI) = m r ·(8.503·OSI 2 −3.741·OSI +0.7449)·C LDL,ox,w C m,l (37) The same procedure was done for the combination of TAWSS and OSI: J s,m (TAWSS, OSI) = m r · (0.8588 · e −0.6301·TAWSS + 0.1295 · e 3.963·OSI ) · C LDL,ox,w C m,l (38)

Plaque Initiation and Growth
The mass balance for open systems can be written as: where ρ o is the total density of the tissue in the reference configuration (Garikipati et al., 2004), inter-conversion of species, implying that the total density in the reference configuration, ρ o , changes with time.
As mass transport alters the reference density, ρ i o , assuming that these volume changes are isotropic, it leads to the following growth kinematicsḞ i g =ρ i o ρ i orig I where ρ i orig means the original concentration of a specie in the undeformed configuration (Garikipati et al., 2004) and I is the second-order unit tensor. For a small strain hypothesis and isotropic growth, we can write: where v is the velocity of the material points. Finally, knowing all substance distributions in the arterial wall, we can compute the growth of plaques. The arterial wall change of volume is due to the contribution of all the cells and substances that are present in the inflammatory process, but the influence of most of them is negligible, so we considered that only larger cells and collagen contribute to plaque formation. Therefore, only FCs, SSMCs and collagen fibres contribute to plaque volume in our model.
In addition, we considered isotropic growth of plaques, so atheroma plaque volume change can be written as: where ∂C i,w is the concentration variation with respect to the initial concentration of the considered substance. Vol FC and Vol SSMC are volumes of an FC and an SSMC, respectively, which can be calculated by knowing their radius. Finally, ρ G is collagen density. Foam cells were assumed to have a spherical geometry, whereas synthetic smooth muscle cells were modelled as ellipsoids, so their volumes can be calculated with Equations (42), (43).
where R FC and R SSMC are the FC and SSMC radii, respectively, and l SSMC its length. The parameters for the growth of plaques process are shown in Table 5.
To validate the results, the stenosis ratio (SR) in the areas with maximum plaque was computed, defining the stenosis ratio as the percent area stenosis in a section. It relates the area of the healthy lumen without the presence of plaque with the area of the lumen with plaque, and can be calculated as: Lumen area with plaque Lumen area without plaque · 100 (44) As can be seen, although the geometries are patient-specific, the parameters are based in literature due to the impossibility of determine their value for each patient. Therefore, there is some variability in them, which was checked to see how it can affect to the model. The parameters that are related to LDL have more influence in the plaque growth given that LDL is the substance that initiates all the process. These parameters were calibrated with the patient "A" and used later for the rest of geometries. In other cases, such as the parameters referred to the cell size, an average value of the parameters given in literature was taken. Regarding measurable parameters for each patient, the most important parameters whose variation would suppose a different behaviour of the model are LDL and monocytes concentration in blood, as well as the arterial pressure of the patient since there are studies that correlate changes in the endothelial permeability as a function of the arterial pressure (Tedgui and Lever, 1984), and other factors, e.g., if the patients are taking medication or not.
In the case of different vascular regions with the same arterial pressure, the parameters that could vary are the hydraulic conductivity of the normal junctions, Lp, nj (Tedgui and Lever, 1984), the monocytes recruitment, mr (Steinberg et al., 1997), as well as the thickness of the arterial wall (Olgac et al., 2008;Sommer et al., 2010).

RESULTS
In Figure 1, we can see the seven real geometries with their corresponding plaques indicated by arrows and the healthy one without plaque. All the geometries, with the exception of E, have large atheroma plaques at both the CCA close to the bifurcation and at the ICA. Geometry E presents a unique plaque in the ICA.
First, we analysed the haemodynamical stimuli effect and the growth on the healthy artery. The results of the simulation for SI and growth of plaques computed with TAWSS, OSI and the new proposed variable are presented in Figure 6. As we expected, small areas with high SI are usually accepted as atheroprones -presented in the healthy geometry -and the corresponding growth is reduced. In particular, for OSI stimulus, negligible growth is presented.
The SI obtained with TAWSS, OSI and the new variable is represented for all the pathological carotid bifurcations in Figure 7. As seen, OSI predicts these areas with high values near the bifurcation but lower than TAWSS, which also predicts areas of high SI in the CCA close to the bifurcation as well as in some areas of ICA and ECA. Finally, the new index also predicts high SI in these areas but in a more localized way than TAWSS. Note that for E geometry, none of the stimuli predicts the location of the plaque.
In Figure 8, growth of plaques after 30 years of the inflammatory process is represented for all the geometries considering the different haemodynamical stimuli. Generally, the location of plaques was better estimated using the new proposed stimulus. For all cases, OSI underestimates the area of the plaques and shows the worst prediction of the location. TAWSS predicts non-physiological growth on the CCA due to the very low values of WSS in this area. For patients A, C and D, the location of the plaque matches the clinical evidence using only the new stimulus and the prediction fails for TAWSS and OSI. For patient B, the new index predicts the location of the plaque of the ICA but fails on the length of the disease. For patient C, the new stimulus matches the location of the plaques on ICA and ECA; however, it predicts plaque on CCA that it is not presented on the clinical images, likely due to an excessive influence of TAWSS in this area. Plaques on the ICA where not predicted by any of the stimuli for patients E and G.
Finally, Table 6 shows the computational stenosis ratio obtained using our mechanobiological model for all the carotid bifurcations and compared with the real ones taken from the clinical images. It can be observed that the stenosis ratio after 30 years was better predicted using the new variable than the other two for the seven studied patients. For example, the stenosis ratio predicted with the new variable for patient C were 52.65 and 22.71 % for the CCA and ICA, respectively, and the real values are 61.14 and 32.91 % for the CCA and ICA, respectively, showing differences of 8.49 and 10.2 % of stenosis, respectively. In general, for all geometries, the new variable has values of maximal and minimal errors of estimated stenosis of 32.89 and 2.77 %, respectively. TAWSS overestimates the stenosis ratio with a maximal and minimal error of 45.74 and 3.58 %, respectively. In contrast, OSI underpredicts the stenosis ratio with an error higher than 9.53 % for all cases.

DISCUSSION
This work extends the mechanobiological model developed by Cilla et al. (2014) from an axisymmetric to a 3D model and is essential as a prior step to its application to patient-specific images. It is worth highlighting that atheroma plaques are usually eccentric, and this feature cannot be captured with 2D axisymmetric models. This works also updates some terms of the equations to enhance its convergence (apoptosis of macrophages into foam cells, differentiation of CSMCs into SSMCs and proliferation of SSMCs, and adding a new term in SSMCs due to their apoptosis) and is solved in 3D instead an axisymmetric formulation. However, the computational cost is higher, and some simplifications are necessary. Another improvement of this model is the computation of blood flow in the transient mode.
This transient mode for the blood flow allows us to analyse and compare some of the most common mechanical stimuli that are normally used for predicting atheroma plaque locations, FIGURE 6 | Shape index distribution (first row) and detail of the growth of plaque in the bifurcation area (second row) in the healthy geometry with the different mechanical stimuli studied. The first column depicts the use of TAWSS, the second column the OSI and the last column the proposed combination of TAWSS and OSI.
TAWSS and OSI. We also propose a new mechanical stimulus as a combination of TAWSS and OSI to better predict the location of plaques. The model was used to check which mechanical stimulus is more appropriate to predict the location of atheroma plaques on carotid geometries.
The model was computed in a healthy geometry without significant plaque in which there is no relevant haemodynamical stimuli to develop atheroma plaques to verify that the model is stable and plaques only grow when there is an atheroprone stimulus. Additionally, atheroma plaque growth was computed in different geometries obtained from pathological patientspecific images in order to obtain the equivalent healthy geometries to validate the results of the model.
As we can see in Figure 8, plaques predicted with different haemodynamical stimulus growth in different locations and present distinct stenosis ratio; thus, the choice of the haemodynamical stimulus to predict the location of plaque results crucial. TAWSS predicts excessively large plaques in the CCA branch of the carotids, far from the bifurcation, that do not match with the real ones. On the other hand, TAWSS adequately predicts the size of plaques appearing in its own bifurcation location compared with the real geometries. In contrast, OSI locates plaques with a better precision than TAWSS, but the growth ratio is very low in comparison to the actual ones. Finally, the new variable proposed in this study combines the results achieved with TAWSS and OSI, predicting in a more adequate way the location of atheroma plaques as in OSI -limiting the growth in areas of the CCA far from the bifurcation-and its growth is similar to the real one observed in the clinical images, similar to TAWSS. The growth model fails to predict the location and size when haemodynamics cannot predict high SI on the location of the plaques. Although most of plaques predicted with the model correspond to the ones clinically observed, there are a few plaques that cannot be explained with this model. For patients E and G, the stenosis appears in a zone where no haemodynamical disturbance is observed; therefore in these patients, haemodynamics are not the main trigger of the disease and could be due to other systemic or genetic conditions of the patient in this location, e.g., external lesion or pathological weakness of the intimal layer. Even though the model correctly predicts the plaques appearing in the specific bifurcation and the ECA, some plaques presented at the ICA are not very well predicted by any mechanical stimulus studied in this work. We hypothesized that it could be due to the uncoupled form between haemodynamics and growth, where fluid is computed at the beginning of the process and it is not updated with the plaque growth. First, plaques appear at the CCA branch, and then they cause a change of the blood flow downstream and the variable stimuli -TAWSS and OSI-could be affected at an area behind them and a new plaque is likely to appear in the ICA. To validate this hypothesis, we computed the growth of B geometry until 15 years, and then we actualized the geometry and computed it again. In Figure 9, we can see the different images that represent growth for a total of 30 years, growth for only half of the time (15 years), reconstruction of an updated geometry of the bifurcation including the stenosis for 15 years, growth for 30 years with the updated geometry and finally the actual clinical image. It can be observed that the plaque in the ICA better matches the real plaque. Due to the high computational cost of the 3D version of our mechanobiological model, full FSI for 3D real carotid images is not available.
As seen, the location of plaques was better predicted than the stenosis ratio, which was predicted with errors between 2.77 and 32.89 %, depending on the analysed geometry. This could be because the location of plaques only depends on haemodynamics, while the stenosis ratio also depends on specific parameters of the patients that are unknown, such as real concentration of LDL in the blood and the time that they have had high LDL levels. Moreover, it is also dependent on specific parameters of the inflammatory model. It could obviously rely on other factors not taken into account in the model that could affect plaque growth, e.g., age, gender, blood pressure level or genetic conditions. However, in general, our model using the new haemodynamic stimulus should predict the location and the growth of the plaques.
Apart from the LDL molecules, many studies have focused on the governing mechanics interaction of the different biological species that play a role in the atheroma plaque development from a computational (see e.g., Ougrinovskaia et al., 2010;Di Tomaso et al., 2011, among others) point of view. Furthermore, there are greatly varying degrees of complexity in these computational studies depending on the number of species considered and the development of the equations proposed. (Zohdi et al., 2004) modelled the adhesion of monocytes to the endothelial surface, which is controlled by the intensity of the blood flow and the adhesion molecules stimulated by the excess of LDL, the penetration of the monocytes into the intima and subsequent inflammation of the tissue, and the rupture of the plaque accompanied by some degree of thrombus formation or even subsequent occlusive thrombosis. Their modelling approach predicts a priori the time to rupture as a function of arterial geometry, diameter of the monocyte, adhesion stress, bulk modulus of the ruptured wall material, blood viscosity, flow rate and mass density of the monocytes. Di Tomaso et al. (2011) considered the interaction between just two species, LDL and monocytes, but the monocyte behaviour was modelled in a very simple way. Fok (2012) proposed a mathematical model of intimal thickening, posed as a free boundary problem. Intimal thickening was driven by damage to the endothelium, resulting in the release of cytokines and migration of SMCs. More complex studies were presented by Siogkas et al. (2011), who included in their model oxidized LDL, macrophages and cytokines, considering that all the LDL molecules and the monocytes were oxidized and differentiated, respectively, at the instant in which these agents pass through the endothelium. A similar study was presented by Calvez et al. (2009) from a mathematical point of view, but their study also included the foam cell. Ougrinovskaia et al. (2010) explored the uptake of cholesterol by different scavenger receptors of macrophages during earlystage atherosclerosis using an ordinary differential equation (ODE) model. It was found that macrophage proliferation rather than an increased influx of LDL particles drives lesion instability. Finally, Bulelzai and Dubbeldam (2012) presented a qualitative mathematical model consisting of a number of ordinary differential equations for the concentrations of the most relevant constituents of the atherosclerotic plaque: macrophages, monocytes, foam cell and oxidized LDL. More complex 2D mechanobiological models have been presented by other authors. For example, Filipovic et al. (2011) used axisymmetric models and thus they could not obtain eccentric plaques. In other studies, such as Filipovic et al. (2013), they used three-dimensional carotid geometries with Kedem-Katchalsky equations and convection-diffusion-reaction equations, but they only considered three substances on the arterial wall and focused the study in the location of plaques and not in their growth. Alimohammadi et al. (2017) and Díaz-Zuccarini et al. (2014) used the aortic bifurcation and the left femoral artery, respectively, in their studies and focused more on the location of plaques than on their growth. Only Alimohammadi et al. (2017) analysed the typical haemodynamical stimuli TAWSS or OSI and proposed a combined index, termed HOLMES, to emphasize regions of highly oscillatory and low-magnitude WSS; however, they focused only on the location of calcifications.
Moreover, it is important to highlight the advantages and disadvantages of our model versus agent-based models (Bhui and Hayenga, 2017;Corti et al., 2020). The main advantage of our model against versus agent-based models is that using a continuum model allows us to simulate the plaque growth in a real complex 3D geometry with an accessible computational cost. Moreover, average values of the parameters and sensitivity analysis are easily implemented. However, the disadvantage is that with the continuum models we cannot take into account the random behaviour of the cells and relation with micro-constituents, which is the one of the main strengths of agentbased models.
The findings of this study should be interpreted within the context of its limitations. For example, our model would be improved by implementing a fluid-structure interaction to better approximate the real pathology. Concerning the mathematical model, only the main processes were included in the model, while other important processes in the development of the atheroma plaque, such as the degradation of collagen with age, were considered. The haemodynamics were considered the main trigger of atherosclerosis initiation. Thus, the cyclic stretch effects of vessel compliance or curvature were disregarded. Another limitation of our model is that we did not have the geometries before plaques were developed, so the real geometries without plaques were unknown. It is a limitation because the real healthy geometries can have some differences with the reconstructed ones and it can have an influence in the obtained results. We did not have real plaque growth monitoring to see plaque geometry evolution, and, in the same way, we also did not have real data from each patient, e.g., blood pressure and flow, LDL levels and number of years with high LDL levels and other pathologies that they may have, so we did not use patient-specific parameters to solve the problem. A fixed plaque growth of 30 years was considered for our analysis based in literature (Insull, 2009). In addition, there are no experimental data about OSI and TAWSS and OSI combination influences on SI, and we used a simplified model in which we only consider the more important substances for plaque growth and do not consider different kinds of cytokines (IL-4, IL-10, IL-13, or TFG beta) or T-cells or free radicals that oxidize LDL. Adjusted parameters come from different species and vessels in which the studies were developed, and they were not all human carotids. We also did not consider other processes that may have an important role in atherosclerosis, such as mechanotaxis. This model assumes that the substances can move from the lumen to the arterial wall, but not in the reverse direction, and the transport properties are set as constant, but in fact, they are very likely to change during plaque formation.
In conclusion, despite the limitations, our model can predict the location and the growth of plaques in the main cases. The results show that prediction of plaques is, in most cases, better using the new mechanical stimulus proposed in this study than using TAWSS or OSI, with a maximal error of 32.89 % on the stenosis ratio computed at the areas of higher occlusion of the lumen due to the plaques. Based on the results, it can be concluded that the functional regulation of the endothelium by local haemodynamic shear stress provides a model for FIGURE 9 | Results for plaque growth in patient B with 30 and 15 years of continuous process (A,B, respectively), the updated geometry after 15 years of the growth process (C), the growth after 30 years in the updated geometry (D) and the real plaque of the patient (E).
understanding the focal propensity of atherosclerosis in the setting of systemic factors and this may help to guide future therapeutic strategies.
In the future, the model could be used to predict if a patient is susceptible to develop atheroma plaques and if so, to determine the places where plaques are likely to appear. In this way, it would be possible to take the necessary treatment to prevent atheroma plaque development and all the consequences derived from atherosclerosis.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation. FUNDING Support was obtained from the Spanish Ministry of Science and Technology through research projects DPI2016-76630-C2-1-R and PID2019-107517RB-I00 and financial support to PH-L from the grant BES-2017-080239, and the regional Government of Aragón support for the funding of the research project T24-20R.