MATHEMATICAL MODELLING OF THE EFFECTS OF STATINS ON THE GROWTH OF NECROTIC CORE IN ATHEROSCLEROTIC PLAQUE

. A large necrotic core increases the risk of atherosclerotic plaque instability. Statins can delay the growth of necrotic core in plaques, but the kinetic mechanism of statins in slowing down the necrotic core has not yet been addressed in detail. In this paper, a mathematical model is governed by a system of advection-diﬀusion-reaction equations coupling of the porous nature of vessel wall is established and applied to illustrate the plaque growth with lipid-rich necrotic core (LRNC) with and without statins using ﬁnite element method. We study the inﬂuence of LRNC plaque growth for diﬀerent drug concentrations at diﬀerent time intervals. The results showed that the drug use at diﬀerent time points has a signiﬁcant impact on the treatment eﬃcacy. Compared with short-term, low-dose treatment, early statin treatment with high dose showed more pronounced eﬀects on reducing the low-density lipoprotein (LDL) cholesterol, decreasing the volume of necrotic core, changing the characteristics of plaques, and improving the plaque stability. The model is validated by comparing with the clinical data, and may be used to predict the progression of LRNC plaque and the eﬀects of statin therapy.


Introduction
Atherosclerosis (AS) in the carotid artery can lead to cardiovascular diseases and its formation mechanism is complicated [1,38,59,62]. It is generally accepted that the development of AS plaque originates from the endothelial permeability to lipid components (such as HDL, LDL, etc.) under blood flow [47]. Low-density lipoprotein (LDL) begins to accumulate in the vessel wall and gradually evolves into oxidized low-density lipoprotein through oxidation reaction [47]. At the same time, foam cells and macrophages are immoderately accumulated in the lumen wall [24,47], and AS plaque is formed consequently. Human anatomo-pathological studies revealed that plaques with a thin fibrous cap and a large necrotic core produced by apoptosis and necrosis of vascular smooth muscle cell and macrophage foam cells [20,39], were more susceptible to rupture [4]. The formation of LRNC is the essential mechanism in the development of the rupture-prone plaque [20,39], and will increase the possibility of plaque rupture when the necrotic core volume accounts for 30-40% of the A number of clinical studies have confirmed that statins have the potential to reduce LDL level, enhance coronary blood flow, and reduce the risk of cardiovascular disease [3,46,55]. Different statins have different pharmacokinetics, mainly considering the safety and clinical application of different drugs [7,10,30], but the most commonly used types of statin drugs are lipophilic [30] (except Pravastatin and Rosuvastatin). Fatos et al. established a pharmacokinetic/ pharmacodynamic (PK/PD) model to describe the effect of statins (Such as simvastatin, fluvastatin, and atorvastatin) on LDL in patients [22]. Bailey et al. reported the therapeutic effect of statin drugs on atherosclerotic plaque in rabbits was similar to those of human inflammatory diseases [41]. Wang et al. studied the effect of rosuvastatin on macrophages in atherosclerotic plaques in mice [28]. Research shows, rosuvastatin drugs can promote the differentiation of M2 macrophages to inhibit the formation of atherosclerotic plaque, and can have the direct effects on the necrotic core development in AS plaques [26]. Hwang et al. found that 54 patients with acute coronary syndrome were treated with different statins after 6 months of continuous treatment, the necrotic core volumes of vulnerable plaques were significantly reduced, and the plaques were more stable after treatment [31]. Masanori et al. investigated the tissue characteristics of coronary plaques using imaging technology after 6 months of continuous treatment with statins, and revealed that the lipid composition of patients with stable angina pectoris was significantly reduced [33]. Banach et al. also found high-dose statin can reduce plaque volume and necrotic core volume while low-dose statin had a relative low effect [6].
In addition to clinical trials and experimental measurements on the pathological process of AS, mathematical modeling of the AS formation and development has gained popularity in recent years because of some advantages of using numerical simulations. Many mathematical models in one-, two-, or three-dimensional spaces have been used to simulate the initiation and formation of atherosclerotic plaques due to chronic inflammatory reaction induced by the diffusion of LDL in the intima [1,9,12,16,40,50,54]. These models included some key substances: LDL and oxidized LDL, and cells: macrophages, foam cells, endothelial cells and smooth muscle cells. Biochemical reactions within blood wall are expressed by a coupled set of advection-diffusion-reaction equations. For example, Calvez et al. [9] used the reaction-diffusion equation to establish a growth model of early atherosclerotic plaques, and considered the diffusion of macrophages, foam cells, oxidized LDL, etc. Silva et al., developed a mathematical model to investigate the early stage of atherosclerosis based on the dynamics of endothelial permeability and wall-shear stress [54]. Younes investigated the inflammatory process upon LDL accumulation and endothelial permeability, especially depicted the effect of anti-inflammatory process. The results indicated that the regulation of atherosclerosis progression was subject to anti-inflammatory responses, which was conducive to plaque regression [1]. Cilla presented a mathematical model to reproduce atheroma plaque growth upon wall shear stress and qualitatively predicted the atheroma plaque development [16]. Base on the positive effect of HDL in AS plaque regression, Chalmers et al. demonstrated the effect of foam cell accumulation on the plaque regression by assumption an increasing HDL influx through the plaque [12]. Sugiyama et al. and Rostam-Alilou et al. investigated hemodynamic characteristics of atherosclerotic lesions in intracranial aneurysms (IAs) using CFD simulation. They reported the influence of the ACA wall on the rupture risk of IAs, the results provide a good reference for clinical treatment [20,39]. In addition, Chambers et al. developed a lipid-structured model of atherosclerotic plaque to demonstrate the progress of macrophage apoptosis, emigration and proliferation using dimensionless modulating functions [13]. However, these computational models mainly focused on the progression of atherosclerotic pathology as mentioned above. But theoretical study or explanation of the highly coupled model interaction of fluid-wall with porosity wall has not yet been attempted. Especially, there is a lack of quantitative understanding about the dynamics of pharmacological drug-target interactions on AS treatment from theoretical models.
The delivery of statin drug involves multiple processes. On the one hand, there are a series of biochemical reactions of statins and molecular substances (such as HDL, LDL, foam cells, and macrophages, etc.) in AS development. On the other hand, penetration of macromolecular substances and macromolecular drugs through vascular endothelial cells serve as determinants in these processes. However, sensible understanding of the mechanisms of action of drugs in the AS treatment based on experimental measurements are not enough. Mathematical models may play an increasingly important role in gaining mechanistic understanding of the effects of statins on the growth of AS plaque and necrotic core in the AS development. Although significant progress has been made over the last decade in mathematical modelling of AS formation and growth in blood flow [9,45,50,56], very few attempts have been made to combine blood flow and drug transport in the AS development, and also accounting for the effect of statin use on AS in two-dimensional space. This article aims to explore the effects of statins on LDL conditions and establish a time-dependent coupling system of partial differential equations to demonstrate the interactions between reacting and diffusing of chemical species under blood flow, and understand the formation and evolution of AS plaque and necrotic core.

Geometrical model
Assume that the model is a section of a single-layer artery wall with a local lesion, which describes the transport of LDL, HDL, and statin drugs through the porous wall. Figure 1 shows the schematic diagram of an artery with a length of 50 mm and a diameter of 6 mm. In addition, the length of the lesion is 10 mm. In the arterial lumen and walls, a primary process involves biochemical reactions like cell migration, proliferation, differentiation and apoptosis, secretion and degradation of substances, and mutual transformation between species.

Fluid dynamics
A two-dimensional form of the Navier-Stokes equation is used to describe the blood flow in the vessel as illustrated in Figure 1. We also assume that the blood is a laminar, incompressible, Newtonian flow with negligible gravitational body forces and a constant viscosity. Suppose that the blood flow at the inlet is directed in x-axis with the velocity profile assumed parabolic and zero velocity in y-axis, the velocity U is given as: where U max = 0.2 m/s is the stable average velocity at the entrance of blood flow, R l is the radius of the blood vessel.

Transport of statins in the vessel wall
Based on the species transport in homogeneous porous medium [2], the convection-diffusion equation of statins in the vessel wall is: Statins are input from the outside in a specific way, and the concentration changes according to equation (2.2), where C Sx,w is the concentration of statins in the blood vessel wall, D Sx,w (= 1.5×10 −14 m 2 /s) is the diffusion coefficient of statins in the blood vessel, ε p represents the porosity of the porous medium of vessel walls, which has a constant porosity of 0.45-0.49 for LDL if the walls are viewed as a multilayer model and normal tissues [17]. In our simulations, ε p is assumed to have a constant porosity of 0.983 [15] when considering a single layer model and tissue lesion, which increases the vascular permeability in inflammation [53]. u w is the convection coefficient of the conserved flux in the walls, and its value is 1e −15 m/s [49]. On the right of the equation (2.2), k in is the supply rate of statins, k out means the rate of degradation of statins, k dec represents the natural decay rate of statins. This article assumes that the entrance and exit of the arterial wall are insulated boundary conditions.

Transport of LDL in the vessel wall
The transport of LDL can be described by the following equation: where C LDL,w and D LDL,w are the LDL concentration and the LDL diffusion coefficient in the vessel wall, respectively, and the value of D LDL,w is 8×10 −13 m 2 /s [16]. Apart from the diffusion and convection terms of equation (2.3), d LDL C LDL,w corresponds to the LDL oxidation formation, indicating the concentration of LDL oxidized per second. Parameter d out is the degradation rate because of drug absorption.

Transport of HDL in the vessel wall
The transport equation of high-density lipoprotein in the vessel wall is [12]: where D HDL,w and C HDL,w are the diffusion coefficient and concentration of high-density lipoprotein in the blood vessel wall, respectively, and the value of D HDL,w is 8×10 −13 m 2 /s. As for the right side of equation (2.4), r h1 represents the rate of high-density lipoprotein absorbed from foam cells due to reverse cholesterol transport, r h2 is the degradation rate of HDL in the blood vessel wall.

Transport of LDLox in the vessel wall
For oxidized low-density lipoprotein, the convection-diffusion equation is expressed as: where C LDLox,w and D LDLox,w are the concentration and diffusion coefficient of oxidized LDL, respectively. The value of D LDLox,w is 8×10 −13 m 2 /s. The D LDLox,w value has been considered as similar to that used for the LDL substances by Prosi et al. [48], and the d LDL C LDL,w models the oxidized LDL process from the LDL. Finally, LDL oxr C LDLox,w C M,w represents the amount of LDL that macrophage consume per second.

Transport of macrophages in the vessel wall
The transport equation of macrophages in the vessel wall is [12,40] Here C M,w is the concentration of macrophage, D M,w is the diffusion coefficient of macrophage in the blood vessel wall. D M,w is assumed to be 1×10 −14 m 2 /s. The first term on the right side of equation (2.6) represents the rate of macrophages in the formation of macrophages by ingesting LDLox, the second corresponds to the rate of foam cells for reverse cholesterol transport (RCT) macrophage pools, and the third is the death and differentiation of macrophage, r m5 C LDLox,w is the increase rate of macrophage.

Transport of foam cell in the vessel wall
The evolution equation of foam cell transport in the vessel wall is given as: where C F,w and D F,w are the concentration and diffusion coefficient of foam cells. Foam cells are usually less diffusive than macrophages [25,37]. To simplify the calculation in the paper, we assume the diffusion coefficient D F,w is assumed to be the same as that of macrophage. The first term at the right end of Equation (2.7) corresponds to the rate at which macrophage transform into foam cells, and then represents the rate at which foam cells undergo RCT return to inflammatory macrophage pools [12]. The parameter combination r F 3 C LDLox,w is the increase of LDLox oxidation, and the last term represents the apoptosis of macrophage [16].

Transport of monocyte chemical attractant in the vessel wall
In this paper, the concentration equation of monocyte chemical attractant agent is given as: where D P,w and C P,w refer to the diffusion coefficient and the concentration of monocyte chemical attractant in blood vessel wall, respectively. We use the parameter D P,w = 1×10 −20 m 2 /s since its diffusion coefficient is very small [16]. On the right side of equation (2.8), r P 1 represents the number of monocyte chemical attractant produced by macrophages consuming oxidized low-density lipoprotein per second, the second term is the number of apoptosis per second for monocyte chemical attractant, and the last term represents the monocyte chemotaxis.

Transport of ES cytokine attractant in the vessel wall
In the vessel wall, the concentration change of ES cytokine attractant reacts according to the following equation: where D Q,w and C Q,w refer to the diffusion coefficient and concentration of ES cytokine attractant in the blood vessel wall, respectively. The diffusion coefficient D Q,w is also very small [16], and is the same as that of monocyte chemical attractant mentioned above. r Q1 represents the rate of ES cytokine attractant produced by macrophages, and the last term r Q2 indicates the loss rate of ES cytokine attractant in the vessel wall. Figure 2 depicts the reactions between various types of substances in equations (2.2)-(2.9), and all reactions take place in the vessel wall. It should be noted that the red boundary represents the damaged boundary (see Fig. 1), and the response is transformed by the formula marked on the arrow. The presentation within this section mainly follows [20,39].

Initial conditions
In the vessel wall, the initial conditions of the above equation at the initial moment are listed below:

Boundary conditions
We establish a mutually perpendicular x and y coordinate system in a two-dimensional coordinate system, Ω = {(x, y) , 0 ≤ x ≤ L, 0 ≤ y ≤ 2R0} where R0 is the radius of the vessel lumen. The boundary conditions In this paper, we assume that the materials (such as LDL and HDL) enter from the damage boundary, and the other boundaries are zero flux. Whereθ is the boundary of the injury at the intima,n is the outer surface of the blood vessel wall, Γ1 and Γ2 represent the boundary of the vessel wall and the intima, respectively. J is the boundary flux, and the values of L 0 and h 0 are listed in the table. The equation (2.14) is the boundary condition of equation (2.2) and we assume that the value of K S is 1×10 −14 m/s. Wall permeability of macromolecules such as LDL/HDL is correlated to the evolutionary process of atherosclerotic lesions [29]. In addition to the blood pressure during the cardiac cycle, shear stress may have a significant influence on wall permeability [29]. In order to simplify the calculation, we assume that HDL and LDL have the same wall permeability [54], and the relationship between shear stress and wall permeability is given as: without the effects of statins: where WSS is the shear stress in the lumen, the value of W SS 0 is 1 N/m 2 [54].

Formation processes of the plaque and necrotic core
In this paper, finite element method (FEM) is used to solve the multi-physical layer coupling problem of convection diffusion reaction equation, and matter transport equation in porous media by using COMSOL Multiphysics [14]. In Figure 1, we use the unit size free triangle mesh to divide the model: the maximum grid size is 1.25 mm, the blood is composed of 1622 mesh vertices and 3018 triangular units, and blood vessel wall is composed of 808 mesh vertices and 1204 triangular units. The boundary conditions of the blood and the intimal boundary (including the lesion boundary) have been given above in the article. The growth area of the patch adopts the ALE moving grid pattern.

Growth model of atherosclerotic plaque
In order to simulate the growth and evolution process of AS plaque in the area modified by macrophages and foam cells, we use the "moving grid" model, that is to say, the grid will move locally due to change with the concentration of foam cells. The mesh gets the smooth deformation in the calculation region by means of a diffusion mode, in which the boundary moves and spreads to the inner area. Meanwhile, we use a Laplace smoothing mesh deformation scheme to represent additional degrees of freedom, which satisfies the following equations [34].
This paper uses the function x = x (X, Y,t) and y = y(X, Y, t) to explore the spatial position of the deformation node (x, y), where t is the time, (x, y) represents the coordinates of the fixed space, (X, Y) is the initial coordinates, and the radial component is given: (2.20) In order to connect the modified high-density LDL to the arteriosclerosis phenomenon of the vessel wall, free parameters, kr and a in equation (2.20), must be set properly because there is a lack of reliable clinical data to study the nature of the internal and external deformation of the vessel wall. We adjust the parameters kr = 0.03 and a = 1/0.014 to get a better description of the phenomenon.

Growth model of necrotic core
Based on AS pathogenesis [20,39], we assume that the necrotic core is mainly composed of dead macrophages and foam cells, and its concentration is expressed as follows. This paper aims to study the development of necrotic core during plaque growth, we assume that when C dead,w reaches a concentration threshold, the necrotic core begins to appear and accumulate.
where C dead,w represents the total concentration of the dead macrophages and foam cells in the plaque, D dead,w is the diffusion coefficient and equals 0. The first term on the right side of equation (2.21) represents the rate of macrophage death, and the second term represents the rate of foam cell death. It is noted that the rate of cell death varies greatly depending on the distance to the endomembrane boundary. Parameter r is the distance to the intimal boundary of the wall. The last term indicates the rate of loss in the blood vessel wall, and the value of d de is 7 × 10 −7 1/s [8].
In order to describe the necrotic core growth in plaque, the necrotic core volume is calculated according to the following equation [8]: where l p is the length of the damage boundary, R (t) represents the height of the necrotic core at time t (the height of red area in Fig. 9).

Drug application and model validation
According to equation (2.2), we assume that the statin concentration in the bloodstream for a 1-day period satisfies the following curve evolution: Figure 4 shows the simulated and experimental results for statins with 80 ng/mL dose every 1rd day when using parameters in Table 1 and drug actions as shown in Figure 3. It found that the statins reduce LDL levels by 68.6% and led to a reduction in plaque volume of 3.6% predicted by our model. A 50.7% reduction in LDL levels and a 5.8% reduction in plaque volume can be observed in experimental results [20,39]. By comparing, the plaque volume from simulation is in accord with the experimental results, but there is a very large difference in LDL levels between them. One possible explanation is that the contribution of different kinds of statins to LDL is complex and is associated with different stages of atherosclerosis in patients. In fact, clinic evidence indicates, the average starting dose of all statins can reduce LDL by approximately 25% to 40%, while higher doses of some stronger statins can reduce LDL by 50% to 60% for long-term drug treatment [20,39].
As another example, Figure 5 shows the comparisons of the AS and necrotic core volume between the simulated and observed results for 40 ng/mL dose every 1rd day. After treatment with moderate-intensity statins for 6 months, our simulation results show that statins reduce the AS plaque volume by 3.3% and the necrotic core volume by 8.6%, respectively. Although the simulation results are slightly lower, they are basically consistent with the experimental values [31], which verifies the rationality of our model.

Model prediction and result analysis
The index values of LDL level, plaque volume and necrotic core volume are viewed as the important indicators of effects assessment on the development of atherosclerosis during statin treatment. Figure 6 shows the LDL concentration-time response curves. For the comparison with no-statin therapy, there is a reduction in LDL levels of 33.57% at 6-month following statin treatment with low doses of 10 ng/ml after 48 months of plaque growth. Interestingly, the LDL concentration change is similar to the concentration-time curves of stains in Figure 3. during statin treatment. The result illustrates that the higher the concentration of reactants, the faster the rate of a reaction will be. Clinical trials indicate that most patients are usually started with a moderate-intensity statin with the potential to lower LDL-C levels by about 30% [20,39]. On the other hand, lowering LDL Levels leads to a corresponding reduction of LDLox as demonstrated in Figure 7. This can be explained well in the LDL-LDLox relation of equation (2.5).
Decreasing LDL levels due to statin therapy has an additional effect on slowing down the transport of LDL through the arterial wall. Next, LDL reduction affects the production of foam cells and macrophages, leading to a slower plaque growth rate. Figure 8 shows the growth behavior of foam cells and macrophages for a "patient" without treatment (red line) versus the effect of the same "patient" with statins (the blue line). Changes in the growing slope of foam cells and macrophages are previously compared with the plaque growing trends from Pichardo-Almarza et al.'s results [20,39].   Figure 9 shows the growth process of the necrotic core in AS plaque in a two-dimensional space. The red part represents the necrotic core area. It can be seen that the necrotic core area gradually expands outward as the plaque grows over time. The area of the necrotic core is 0.52 mm 2 at 1170 days, and increased by 5.49 mm 2 at 1440 days and 13.87 mm 2 for 1800 days. The plaque growth at lesion has an unexpectedly direct influence on the development of AS plaque in the healthy wall regions, where the necrotic core area of the healthy wall is 5.85 mm 2 at 1800 days. This phenomenon has never been illustrated in previous numerical simulations, but confirmed at autopsy [51] and imaging research [44,51,58].

Model discussion
The formation and growth of the necrotic core in AS plaque is closely associated with macrophage apoptosis and foam cell apoptosis, which is usually differentiated by macrophages [20,39]. The effectiveness of statin use in the inhibition of macrophage differentiation and the associated process of the necrotic core growth is influenced  by serval factors, including dose [44], duration of treatment [6], frequency of administration [18], biological properties of the microvasculature, etc. Figure 10 shows the curves of the necrotic core volume at the lower vessel wall after 6 months of statin treatment during different treatment cycles. For the comparison with nostatin therapy, applying statins for 1 days, 2 days and 3 days reduce the necrotic core volume by 6.06%, 4.61%, and 4.35%, respectively. After 1 year, the necrotic core volume decreases by 14.82%, 10.51% and 9.42%, respectively. The treatment cycle is negatively correlated with the necrotic core volume. Figure 11 shows the evolution of necrotic core volume in AS plaque with and without statin use. Compared with no drug intervention, after 6 months of treatment with statins at concentrations of 2 ng/mL, 10 ng/mL, and 20 ng/mL, the necrotic core volume is reduced by 4.51%, 6.06%, 6.86%, respectively. After 1 year, the necrotic core volume was reduced by 10.46%, 14.82%, and 16.03%, respectively. The results indicate early statin treatment has a great impact on the necrotic core growth.    The curves in Figure 12 illustrate the necrotic core volume at the lower vessel wall after continuous treatment with statins for 6 months at different points over time. Compared to the curve without statin intervention, the necrotic core volume of the lower vessel wall decreased by 34.83%, 28.61%, and 14.82% at 1800 days, when treatment started at 720 days, 1080 days, and 1440 days, respectively. The results show that treatment effect is Figure 11. The growth curves of the necrotic core volume in the inferior blood vessel wall for three different statin concentrations. The main goal of this paper is to develop a model that provides enough flexibility to the "user" for creating different case scenarios to test the pharmacological hypotheses in the study of atherosclerosis, and the effect of statins on the progression of atherosclerotic plaque to treatment. One of the main advantages of this approach is the possibility to expand or simplify the model according to different assumptions or requirements, made by the user. The model can either predict the plaque formation and development in coronary arterial segment, or predict the treatment status of the user with statin therapy. It should be noted that this article also has some limitations. For example, there is no inclusion of the mechanical properties of vessel walls in the vessel deformation. Multilayer vascular-specific structure (media and adventitia) is not included in the process of substance transportation, and the effect of hypoxia and angiogenesis, which have a significant role in necrosis and drug delivery is neglected. In addition, the model does not take into account the development of fibrous cap, which is associated with the high residual risk of atherosclerosis in later stages of the atherosclerotic process. We will improve the computation model and provide a more accurate prediction of complex pathological process in the atherosclerotic plaque with and without statin sue in a future study.

CONCLUSION
In this work, we established a mathematical model, which combining chemical reactions, fluid dynamics, and mass transfer in porous wall. The model took into account the pharmacodynamics and pharmacokinetics of statin drug, and was solved by using finite element method to describe the AS plaque growth and to evaluate the effect of statin therapy on the growth and evolution of the necrotic core in AS plaques.
It was found that the model supports previous experiments, in which statins slowed down the volumes of plaque and necrotic core, and explained the mechanism of the process. The therapeutic effect of statins depended on the dose of the drug, treatment time [19]. High-frequency and high-dose statins were highly effective in slowing down the lipid-core plaque growth. Especially, early statin treatment had a great impact on regression of AS plaque rupture, causing serious cardiovascular diseases such as a heart attack or stroke. In addition, the model could be used to predict the volume growth of AS plaque and necrotic core, and the volume changes for various drug diffusion problems in pharmaceutical studies. The study might have many applications in drug control, drug dosage and other related problems of pharmaceutical research.