**2. Geometry**

Due to the high computational requirements of the model and because the objective of this study was to analyse the influences of the different parameters on plaque growth, a 2D axisymmetric approximation was considered enough for this study.

The geometry was developed based on that of Olgac et al. [11] as it reproduces the mechanical stimuli to which real patients are subjected, leading to the appearance of plaque (*TAWSS* and *OSI*). The geometry was adapted to carotid arteries by modifying the diameter and thickness of the vessel to the corresponding values for carotid arteries: 3.63 mm and 0.7 mm, respectively [22]. Figure 1 shows the geometry, which is composed of the lumen and the arterial wall, considered a monolayer. It also has an obstacle plaque that disturbs blood flow, causing low *TAWSS* and high *OSI* downstream and favouring the appearance of a new plaque in that area. It is a phenomenon that was observed in real patients [23,24]. This second plaque is the one that is going to be analysed.

The model was discretised using the finite element method. A sensitivity analysis of the length of the geometry and the mesh was developed, obtaining that the minimal length of the artery to ensure a fully developed blood flow is 90 mm. The domains of the model were meshed using quadrilateral elements, with three boundary layers in the lumen near the endothelium. A sensitivity analysis of the mesh was performed in both lumen and arterial wall meshes, in order to determine the most suitable mesh to solve the model with minimal computational costs, but without affecting the results. The total number of elements is equal to 36,731 for the lumen and 36,036 for the arterial wall, resulting in a total number of elements in the geometry of 72,767. A discretisation with P1-P1 elements was made for the calculation of blood flow. Therefore, linear interpolation was considered for both the velocity and pressure of blood flow in the lumen of the artery. To perform

the sensitivity analysis of the mesh, a coarse mesh was selected in the first place and progressively refined, until the results obtained did not vary more than 5%. The same procedure was done with the number of boundary layers of the lumen mesh. In addition, linear and quadratic serendipity elements were used for the inflammatory process in the arterial wall and for plasma flow through the endothelium and the growth of the plaque, respectively. Two details of the selected mesh are shown in Figure 1. The solution to the transient problems was obtained using the Backwards differentiation formula (BDF) method, which is implicit, following Newton's method for nonlinear problems.

**Figure 1.** Two-dimensional axisymmetric geometry composed of lumen (in red), the arterial wall, and obstacle plaque (both in white color). Details (**A**,**B**) show the mesh near the endothelium at both sides (lumen and arterial wall) for the area near the obstacle plaque (**A**) and the area of development of the new plaque (**B**).

The total number of degrees of freedom of the lumen can be determined by multiplying the components of the velocity of blood and the number of elements in the lumen mesh. For the case of the arterial wall, the number of degrees of freedom would be equal to the multiplication of the components of the plasma flow velocity, the number of substance concentrations that we calculate and the number of elements of the arterial wall mesh.

### **3. Mathematical Model**

In this section, the mathematical model is described, showing the equations in different subsections, referring to blood flow in the lumen, plasma flow across the endothelium, the inflammatory process in the arterial wall, and plaque growth.

### *3.1. Blood Flow in the Lumen*

Blood flow along the arterial lumen is considered laminar due to the Reynolds number in the arteries under physiological conditions (*Re* ≈ 950 for the mean diameter of the considered artery) [25]. Furthermore, blood is modelled as a Newtonian incompressible fluid, since the diameter of the lumen in the arteries is greater than 0.5 mm [26,27]. Moreover, it is modelled as a homogeneous fluid due to the small size of the particles with respect to the lumen diameter of the arteries [25].

The Navier–Stokes and continuity equations govern blood flow along the lumen:

$$
\rho\_b(\boldsymbol{u}\_l \cdot \nabla)\boldsymbol{u}\_l = \nabla \cdot \left[ -P\_l I + \mu\_b (\nabla \boldsymbol{u}\_l + (\nabla \boldsymbol{u}\_l)^T) \right] + F\_l \tag{1}
$$

$$
\rho\_\flat \nabla \cdot \mathbf{u}\_I = 0,\tag{2}
$$

where subscripts *b* and *l* refer to blood and lumen, respectively. The parameters *ρb* and *μb* are, respectively, the density and dynamic viscosity of blood, while *ul* and *Pl* are the velocity and pressure of blood flow. *Fl* is the internal force of blood, which is negligible compared to the friction of blood flow with the arterial wall. All parameters referring to blood flow along the lumen are in Table 1.

**Table 1.** Parameters to calculate blood flow.


Blood flow is modelled in a transient step. An analysis of the number of cardiac cycles needed to model blood flow was performed, determining that three cardiac cycles are enough to develop blood flow in the lumen.

Transient mass flow and pressure are imposed, respectively, at the inlet and outlet of the geometry as boundary conditions [25], as can be seen in Figure 2. These boundary conditions influence the location and size of the generated plaque as they directly affect the recirculation that the obstacle plaque produces downstream. Therefore, they directly affect the mechanical stimuli that initiate plaque growth. For example, when increasing the inlet mass flow, the recirculation is bigger and, therefore, the areas with the mechanical stimuli that initiate plaque growth move to the end of the geometry. Finally, a no-slip condition is imposed at the endothelium.

**Figure 2.** Mass flow (black line) and pressure (blue line) imposed as boundary conditions at the inlet and outlet of the geometry, respectively [25].

### *3.2. Plasma Flow across the Endothelium*

As a result of the porous nature of the arterial wall, plasma can enter the bloodstream, causing a plasma flow through the endothelium. Plasma flow across the arterial wall is modelled with Darcy's law and continuity equations [11,29–33]. Darcy's law relates the

velocity of the plasma flow with the pressure drop in the arterial wall and considers some parameters related to the porous nature of the arterial wall (permeability and porosity).

$$
\mu\_w = -\frac{k\_w}{\mu\_p} \nabla p\_{w\prime} \tag{3}
$$

$$\frac{\partial(\rho\_p \epsilon\_w)}{\partial t} + \rho\_p \nabla \cdot u\_w = f\_v \tag{4}$$

Parameters related to the arterial wall are *uw*, *kw*, *w*, and ∇*pw*, and represent the velocity of the plasma flow, the permeability and porosity of the arterial wall, and the pressure gradient, respectively. *ρp* and *μp* are the dynamic viscosity and density of the plasma. Finally, *Jv* is the plasma flow across the endothelium, which is calculated with the Kedem–Katchalsky equations [34] and the three-pore model [11]:

$$J\_v = Jv\_{vj} + Jv\_{lj} + j\mathfrak{p}\zeta\_p \tag{5}$$

As can be seen in Equation (5), the three-pore model considers three different types of plasma flow through the endothelium: plasma flow through normal junctions (*Jvnj*), leaky junctions (*Jvlj*), and vesicular pathways (*Jvvp*). Nevertheless, vesicular transport is intended for molecular transport, so plasma flow is negligible, relative to that across normal and leaky junctions.

Plasma flow through normal and vesicular pathways can be defined as:

$$Jv\_{n\rangle} = Lp\_{n\rangle} \cdot \Delta P\_{\to n d} \cdot (1 - \Phi\_{lj}) \tag{6}$$

$$Jv\_{lj} = Lp\_{lj} \cdot \Delta P\_{End} \tag{7}$$

*Lpnj* and *Lplj* depend on the considered artery and are the hydraulic conductivity of normal and leaky junctions [35], while Δ*P* is the pressure drop across the endothelium, which depends on the intraluminal pressure [35]. Finally, <sup>Φ</sup>*lj* is the endothelial fraction occupied by leaky cells [36–38].

The hydraulic conductivity of leaky junctions can be calculated as [38,39]:

$$L p\_{lj} = \frac{Ap}{S} \cdot L p\_{slj\prime} \tag{8}$$

with *ApS* being the fraction of the total area occupied by leaky junctions, and *Lpslj* being the hydraulic conductivity of a single leaky junction. To determine these parameters, we considered that leaky cells are enclosed by leaky junctions, which have ring shapes and are randomly distributed with a distance equal to *lj*, being <sup>2</sup>*lj* of the permeability of a leaky junction. Moreover, the spaces between endothelial cells are cylindrical [38,39]. Therefore, *Ap S* can be calculated as:

$$\frac{A\_p}{S} = \frac{A\_{slj}}{\pi \cdot \mathbf{e}\_{lj}^{2}},\tag{9}$$

where *π* · *lj*<sup>2</sup> is the total area occupied by leaky junctions, and *Aslj* is the area of a single leaky junction that can be calculated, assuming reduced thickness, as:

$$A\_{slj} = 2\pi \cdot R\_{celll} \cdot 2w\_{l\prime} \tag{10}$$

where *Rcell* is the radius of endothelial cells and *wl* is the half-width of a leaky junction [36,38,39]. On the other hand, <sup>Φ</sup>*lj* can be defined as [36,38]:

$$
\Phi\_{lj} = \frac{{R\_{croll}}^2}{{}^2\_{lj}} \tag{11}
$$

Combining Equations (9)–(11), the fraction of area occupied by leaky junctions can be rewritten as (12):

$$\frac{A\_p}{S} = \frac{4w\_l}{R\_{cell}} \cdot \Phi\_{l\bar{j}} \tag{12}$$

The ratio <sup>Φ</sup>*lj* is known to be a function of some mechanical stimuli, such as the wall shear stress (*WSS*), the time-averaged wall shear stress (*TAWSS*), and the oscillation shear index (*OSI*) [40]. As blood flow is modelled in transient mode, *WSS* was not considered due to its changes according to the cardiac cycle. Instead of *WSS*, *TAWSS* was defined. This variable takes into account the tangential stresses that blood produces in the endothelium during a cardiac cycle and that causes changes in the shape index of the endothelial cells. Therefore, in areas with low *TAWSS*, endothelial cells will be rounder, allowing plasma and substances to flow between their pores [16,40–42]:

$$TAWSS = \frac{1}{T} \int\_{o}^{T} |\mathbf{r}(t)| \cdot dt,\tag{13}$$

where *T* is the period of a cardiac cycle and |*τ*(*t*)| the magnitude of *WSS* dependent on time. On the other hand, high values of *OSI* imply that blood is highly oscillatory there. Thus, the tangential stresses of blood flow in the endothelium are also oscillatory and, due to their variation, they cause a change in the shape index of the endothelial cells, making them rounder and allowing the flow of plasma and substances through the endothelium, being areas of high *OSI* considered as atheroprones [8,40]. *OSI* is defined as:

$$OSI = 0.5 \cdot \left(1 - \frac{|\frac{1}{T} \int\_{o}^{T} \tau(t) \cdot dt|}{TAWSS} \right) \tag{14}$$

Moreover, the ratio <sup>Φ</sup>*lj*, *TAWSS*, and *OSI* are related by experimental correlations [6–8,11,40]. The first correlation is a parameter named shape index (*SI*), which can be calculated as:

$$SI = \frac{4\pi \cdot Area}{Perimeter^2} \tag{15}$$

where *Area* and *Perimeter* correspond to a single endothelial cell.

To determine the behaviour of endothelial cells with *TAWSS* and *OSI*, we proposed in a previous work [8] a numerical correlation based on the experimental results of Levesque et al. [40].

$$SI = 0.0264 \cdot e^{5.647 \cdot OSI} + 0.5513 \cdot e^{-0.1815 \cdot (TAWSS)^2} \tag{16}$$

The total number of mitotic cells (*MC*) in the endothelium depends on the *SI* [43], so the next experimental correlation was developed, with a unit area of 0.64 mm2:

$$\text{MC} = 0.003797 \cdot e^{(14.75 \cdot SI)} \tag{17}$$

Finally, the number of leaky cells (*LC*) and mitotic cells (*MC*) is correlated by the following empirical equation [11,44]:

$$L\text{C} = 0.307 + \frac{0.805 \cdot \text{MC}}{0.453} \tag{18}$$

On the other hand, <sup>Φ</sup>*lj* can be defined as the ratio between the area of leaky cells and the total area of cells. It can be calculated as:

$$\Phi\_{lj} = \frac{L\mathbb{C} \cdot \pi \cdot R\_{cell}^2}{A\_{unit}} ,\tag{19}$$

considering *Aunit* equal to the unit area in all the experimental correlations before (0.64 mm2). By calculating blood flow in the model, it is now possible to determine <sup>Φ</sup>*lj* and therefore *Ap S*.

Finally, the hydraulic conductivity of a single leaky junction, *Lpslj*, is defined according to Olgac et al. [11]:

$$L p\_{slj} = \frac{w\_l^2}{3 \cdot \mu\_p \cdot l\_{lj}},\tag{20}$$

where *wl* and *llj* are, respectively, the width and the length of a leaky junction.

So, the plasma flow across the endothelium is completely defined, and the parameters needed to calculate it are shown in Table 2.

The boundary conditions of the plasma flow across the endothelium are the normal velocity of the plasma flow (*Jv*), which has already been defined, and the pressure in the adventitia, which [11] must be defined. Both an increase in the normal velocity of the plasma flow and a decrease in pressure at the adventitia would produce an increase in plasma flow, which will also cause an increase in substances flow, with the consequent increase in plaque size [8].

**Table 2.** Parameters to model plasma flow across the endothelium.


#### *3.3. Inflammatory Process in the Arterial Wall*

There are many substances involved in the inflammatory process of the arterial wall. In the model, we consider molecules, such as low-density lipoproteins (LDL) and oxidised LDL (oxLDL). We also consider cells, such as monocytes (m), macrophages (M), foam cells (FC), and contractile and synthetic smooth muscle cells (CSMC and SSMC). We also model cytokines (C) and collagen fibre (Cg).

All of these substances can have convection and diffusion and interact with the others, so the inflammatory process of the arterial wall can be described with convection–diffusion– reaction equations.

In addition, the flow of substances across the arterial wall can be defined as:

$$N = -D\_{X\_i} \nabla X\_i + \mu\_w X\_{i\prime} \tag{21}$$

where *DXi* , *Xi*, and *uw* are the diffusion coefficient of the considered substance, its concentration and its convection velocity in the arterial wall. Due to the structure of the arterial wall, the diffusion is anisotropic, as the longitudinal diffusion is 3 times higher than radial diffusion [39].

At the beginning of the process, in the healthy artery, all muscle cells in the arterial wall are contractile. Thus, the only substance that has an initial concentration in the arterial wall is CSMC.

A high level of LDL concentration in blood is considered, according to a pathological model (2.7 *mg ml* [47], which corresponds to 6.98 *mol m*<sup>3</sup> ). Moreover, we impose a physiological monocyte inlet concentration in the lumen of 550 · 10−<sup>9</sup> *cells m*<sup>3</sup> [48]. An increase in both inlet concentrations would cause an increase in the growth of the resultant plaque.
