Stability analysis of a model of atherosclerotic plaque growth.

Atherosclerosis, the formation of life-threatening plaques in blood vessels, is a form of cardiovascular disease. In this paper, we analyze a simplified model of plaque growth to derive physically meaningful results about the growth of plaques. In particular, the main results of this paper are two conditions, which express that the immune response increases as LDL cholesterol levels increase and that diffusion prevails over inflammation in a healthy artery.

1. Introduction

Atherosclerosis is a type of cardiovascular disease characterized by the accumulation of fatty deposits within the arterial intima, which can eventually expand into the lumen and obstruct blood flow. Stiff plaques are susceptible to rupture due to blood shear stress, leading to the formation of blood clots which can block the lumen entirely and lead to myocardial infarction (heart attack) [1]. According to the Center of Disease Control, cardiovascular disease is the leading cause of death in America, leading to 1 in every 4 deaths. Furthermore, 49% of all Americans possess at least one of the three major risk factors for cardiovascular disease: smoking, high low density lipoprotein (LDL) levels, and high blood pressure [2]. Thus, examining which factors in the plaque formation process are most important in causing plaque growth is of utmost importance in being able to treat patients with cardiovascular disease.

Atherosclerosis can be modeled as a multistage disease, starting with inflammation, leading to plaque growth and finally to plaque rupture and cardiac infarction [3]. Herein, we focus on inflammation and deem that plaques are stable if the inflammation cascade culminates in a return to the equilibrium state.

In this work, we study two types of lipids: high density lipoproteins (HDL) and low density lipoproteins (LDL). Both are composed of dense lipids covered by a surface protein, with molecules such as vitamin E attached for defensive purposes (3–15 for LDL molecules, 0-1 for HDL molecules) [4]. When free radicals come in contact with HDL or LDL molecules, they preferentially oxidize the vitamin E molecules attached to the protein coat. When these outer molecules have already all been oxidized, the core is oxidized. This forms oxidized LDL (ox-LDL), which acts as an irritant to the body, as well as oxidized HDL (ox-HDL), which does not seem to have any adverse effects [3]. However, HDL molecules are six times more prevalent than LDL molecules in the bloodstream and can thus act as a buffer against the buildup of ox-LDL [4].

The presence of ox-LDL causes endothelial cells to release monocyte chemoattractant protein (MPC-1), which attracts monocytes into the intima. These monocytes mature into macrophages, which phagocytize ox-LDL molecules [5]. These macrophages lose motility after having ingested the offensive molecules and turn into so-called foam cells [6]. Foam cells then cause positive feedback, releasing chemokines that lead to the recruitment of more macrophages [5]. Prompted by growth factors emitted by macrophages, smooth muscle cells from the media migrate into the intima [3].

Foam and smooth muscle cells within the intima eventually necrotize, breaking apart into detritus that coagulates into the fatty core of a plaque. This core is covered by a brittle cap of extracellular matrix (ECM) excreted by smooth muscle cells. As the plaque expands into the artery over time, blood stress on these types of plaques, called thin-capped fibroatheromas, along with the degrading effects of matrix metalloproteinases produced by macrophages, causes the plaque to rupture. The rupture of plaque releases blood clotting factors, with blood clots potentially blocking all blood flow and causing heart attack and death. Rupture, rather than complete blockage of the artery, is estimated to cause a vast majority of heart attack cases [6].

From a mathematical standpoint, much work has been done in creating realistic simulations of plaque formation via systems of partial differential equations, with novelties including different types of grids [6], inclusion of many different types of molecules [5], and the inclusion of fluid structure interaction [7].

However, comparatively little has been done in the way of mathematical stability analysis. El Khatib et al. studied a system of two partial differential equations in 2 dimensions and showed that small deviations from equilibria can lead to travelling wave solutions leading to runaway inflammation [8]. Ibragimov et al. linearized their system of differential equations and developed an ansatz based on the eigenfunctions of the von Neumann problem to find specific conditions for stability [9]. They found that the equilibrium point of the system is unstable when disease effects dominate diffusion. The primary drawback of this approach is that the ansatz relies on certain initial conditions and thus does not fully capture the dynamics of plaque growth.

Energy estimate analyses are more applicable in proving stability for large systems of differential equations, where relying on an ansatz can become unwieldy. Ibragimov et al. performed an energy stability analyses in [10, 11], respectively, and derived conditions sufficient for linearized stability. However, this approach entails the imposition of a large number of conditions, many of which can not be readily interpreted physically. The main novelty of this work is to consider a somewhat simplified differential equation system that still captures the key aspects of the plaque formation process and to derive 2 simple and sufficient criteria for linear stability.

The presentation of this work is as follows. Governing differential equations for a model of plaque growth along with the formulation of the computational domain are outlined in Section 2. Section 3 presents a rigorous stability analysis for the linearization of the differential equations derived as well as sufficient conditions for the system equilibrium to be asymptotically stable. Section 4 then numerically validates the criteria found in Section 3.

2. Mathematical Modeling

2.1. Domain of Interest

Here, we study the evolution of species concentrations in the intima, as pictured in Figure 1. We assume that the intima is of fixed volume. This is a reasonable assumption as we are studying perturbations from equilibrium, more likely when total plaque growth is limited to relatively small fibrous atheromas, which do not extend extremely far into the intima. As in [9], we take the intima to be approximately shaped like a thin annulus, with some intimal thickening. We assume the outer boundary of the intima (separating it from the media) is impermeable to all chemical species.

Figure 1

Computational domain.

2.2. Modeling Assumptions

We make a number of simplifying assumptions to obtain a differential equation system capturing the key aspects of the plaque formation process. The most important aspect is that we only consider the species LDL/ox-LDL (L), macrophages (M), and necrotic lipids (N). We assume that all LDL molecules entering the intima are quickly oxidized, becoming ox-LDL. We also assume that monocytes entering the intimal region quickly mature into macrophages. We assume that macrophages are attracted directly to the presence of ox-LDL molecules. We assume that the time span between the creation and eventual necrosis of foam cells is relatively short in comparison to the length of the disease progression. Finally, we assume that neutrophils and other immune cells clear lipid accumulation at a rate proportional to the concentration of lipids present.

2.3. Governing Differential Equations

We model plaque growth with a system of reaction-diffusion equations of the form

We assume that the concentration flux for species f is of the form (where f can represent LDL/ox-LDL (L), macrophages (M), and necrotic lipids (N)) with terms addressing the effects of Keller-Segel chemotaxis (towards species g) and Fickian diffusion, respectively. D is the coefficient of diffusion, and μ is the chemotactic sensitivity coefficient of f towards g [10]; both are assumed to remain constant for any fixed species across the domain due to relative intimal homogeneity. Macrophages are the only species considered that responds to chemotactic influences, so we assume that μ = μ = 0. Necrotic lipids are assumed to be immobile due to their dense structure, and so we take D = 0.

The principal contributor to the reaction terms is the reaction L + M → N. We make the simplifying assumption that this reaction is 1st order in each reactant (i.e., that the rate of formation of necrotic lipids ∝LM), with reaction rates of k and k . We also assume that the rate of formation of lipids is the same as the rate of destruction of macrophages (i.e., that LDL compose an insignificant portion of the lipid mass) and that neutrophils are already dispersed throughout the intima and will destroy plaque at a rate proportional to the amount present with reaction rate k .

Finally, for mathematical simplicity, we impose no-flux conditions on the boundary. As in [10], we instead add source terms to our system to represent LDL and macrophage inputs into the system. We model LDL molecule input to be constant at a rate of L 0 and consider a variable influx of macrophages depending on the LDL concentration: f(L), where f is sigmoidal to account for saturation as the LDL concentration increases without bound.

With the above considerations, we have the following governing differential equation system:

3. Stability Analysis

In order to study the asymptotic stability of system equilibria, we linearize our system, introducing the equilibrium values of all reactants in the process, L , M , and N . We letwhere the equilibrium values, in terms of system parameters, are

for perturbations a, b, and c. We now substitute these into (2), effecting a multivariate taylor series expansion of the different terms and preserving solely linear terms (it can be seen that the constant terms cancel due to the stipulations on the values of L , M , and N ). We then obtain the system where we have the parameters

Next, we will perform a complete energy estimate analysis on our linearized system.

Our first step is to multiply the equation for all species f by f and integrate by parts over the domain using Green's identities with no-flux boundary conditions. We obtain

Now we prove a lemma deriving a bound on an energy functional.

Lemma 1 .

Given equations (9), one has that

Proof

In order to separate products of terms, we apply Cauchy's inequality, which states that for any a, b ∈ ℝ and ϵ ∈ ℝ+ Thus, for the weights ϵ   (i ∈ {1,2, 3,4, 5,6}), we have the system of inequalities: Now, we add together (12), (13), and (14), weighted by x, y, z ∈ ℝ+, respectively, to obtain the desired result.

Stability requires all perturbations to go to zero at all points in space as time increases. We first note that it is a necessary (but not sufficient) condition for lim∫f = 0 for all species f involved.

Theorem 2 .

The ordinary differential equation system associated with (5)–(7) is stable if and only if k 1 k 4 > k 2 k 3.

Integrating over the entire region and using the homogenous Neumann boundary conditions, we have the ordinary differential equation systemwhere f = ∫f for any chemical species f. In order for our ODE system (15) to be stable, the eigenvalues of the above matrix must all be negative [12]. The characteristic polynomial of the above matrix, A, isThe zeroes of this equation are λ 1 = −k 5, which is always negative, as well asClearly both of these zeros are real, with taking the negative sign above leading to a negative value for λ 2. In order for the third zero to also be negative, we require

We now build up a series of lemmas in order to prove the existence of an everywhere positive energy functional that decays to zero as time increases without bound. Herein, we assume thatwhich follows from Theorem 2, as well asand the reasons for the sufficiency of this condition will become more apparent via the following lemmas.

Lemma 3 .

There exist x, y, and z positive such that one simultaneously satisfieswhere one sets ϵ 1 = k 6/k 5, ϵ 2 = k 4/k 5, and ϵ 4 = ϵ 5 = k 1/2k 2.

We may reformulate the stability condition derived in Theorem 2 asAs all of the k 's are positive, we may then writeRearranging, we may then say thatAs we have set ϵ 4 = ϵ 5 = k 1/2k 2, this is equivalent to

As this inequality holds for all parameter values in our parameter space, we may choose x and y to be positive constants such that the ratio y/x lies in between these two bounds. Then,

Now, ϵ 5 = k 1/2k 2 implies that k 1 − k 2 ϵ 5 > 0. Similarly ϵ 4 = k 1/2k 2 along with the stability condition derived in Theorem 2 implies that 4k 4 − k 3/ϵ 4 > 0. This means we may multiply out terms to obtain the system of inequalities:

Then, we have thatWe can then choose z to be some positive constant less than the right hand sides of the above equation. Rearranging, we then have giving the desired result.

Lemma 4 .

There exists ϵ 3 positive such that one can satisfy

From (20), we may assert the bound:Choosing ϵ 3 such that it lies between these two bounds, we have

Rearranging this inequality using the bound on y/x derived earlier gives the desired result.

The next theorem describes the conditions for asymptotic stability.

Theorem 5 .

The equilibrium solution (L , M , N ) is asymptotically stable given that

f′(L ) > 0,

χ 2 < 4D D k 3/k 2.

We define an energy functional, ℱ(a, b, c) = ϕ(t), equal to the left hand side of  (10). From lemmas provided earlier, we know that there exist constants ϵ and x, y, and z such that each of the coefficients in the definition of the functional is positive. As each of the coefficients on the right hand side is negative given the appropriate choice of constants, we may conclude for some positive M (the absolute value of the maximum of the coefficients of the RHS) thatThis implies that the functional decays exponentially and that it is thus asymptotically stable.

The first criterion that f′(L ) > 0 is necessary for the equilibrium of the linearized system of partial differential equations to be stable. Physically, as f is the macrophage response to an LDL presence in the intima, this condition implies that an increased number of macrophages must permeate the intima if LDL/ox-LDL levels increase. This is characteristic of a healthy immune response, as the presence of LDL causes endothelial cells to release monocyte chemoattractant protein, attracting monocytes that later mature into macrophages [5].

The second criterion that χ 2 < 4D D k 3/k 2 represents loosely that the force of chemotaxis (represented by χ ) should be countered by diffusive forces (represented by the product D D ). As chemotaxis results in an aggregation of chemical species in one place, a large value of χ could result in a runaway immune cascade leading to plaque growth and eventual rupture. A strong propensity towards diffusion indicated by large values of the diffusion coefficients would temper this aggregation by opposing the concentration of species in one region.

4. Numerical Results

In this section, we evaluate the robustness of the stability criteria:

4.1. Bipolar Coordinates

Similar to [6], we implement our original differential equation system (2) in a 2-dimensional bipolar coordinate system. The conversion between bipolar and rectangular coordinates is given by where curves of constant τ describe nonconcentric circles (with greater τ corresponding to smaller radius) and where the angle σ can vary between −π and π.

The computational domain, τ < τ < τ int and −π ≤ σ ≤ π, is thus an off-center annulus and so meets the requirements for our modeling assumptions to hold. We implement our equations in MATLAB using the method of lines, discretizing in the space variables σ and τ and solving a system of ordinary differential equations, with each equation corresponding to a grid point in our discretization.

Specifically, our grid points arewhere we take σ = σ for consistency.

Denoting the value of L at point (τ , σ ) by L ( (and similarly for M and N), we have the system of equations: where we discretize the Laplacian operator asand similarly for M ( for (i = 2,3,…, m − 1) and (j = 1,2,…, n).

At the boundaries τ int and τ , we employ the one-sided approximation:where F and G are the τ and σ components of the flux. Due to the no-flux boundary conditions, we may assert thatA similar approximation scheme is used for M at the boundary.

4.2. Numerical Experiments

Setting parameter values that satisfy the second stability criterion, we first attempt to investigate the robustness of the strict first criterion. Letting

our previous stability results should show that I → C as t → ∞, where C is some constant (as the perturbations from equilibrium should approach 0). Approximating this integral as a Riemann sum over grid points over time, we have numerical confirmation that this criterion holds strongly. For instance, the first two images are for f′(L ) > 0 while the latter two are for f′(L ) < 0, starting with different initial conditions.

Mathematically, f′(L ) < 0 implies the existence of an eigenvalue with positive real part of the coefficient matrix A in the proof of Theorem 2. It is well known that this causes the space-integral of at least one of the reactants to diverge over time in the context of the linear system. f′(L ) > 0 is thus a necessary condition for stability of the linearized system and seems to also be necessary for the stability of the initial system, as illustrated by Figures 2–5.

Figure 2

Spatial integrals of squared perturbations over time for f′(L ) > 0.

Figure 5

Spatial integrals of squared perturbations over time for f′(L ) < 0.

On the other hand, setting parameter values that satisfy the first stability criterion and investigating the second stability criterion results in no such confirmation: the criterion does not seem to be necessary, only sufficient.

5. Conclusion

The model that we have developed encompasses the salient features of atherosclerotic plaque development, including lipid oxidation, foam cell formation, and necrosis, as well as plaque buildup. We provide sufficient conditions for a stable artery during the first inflammatory stages of atherosclerosis. Specifically, our first stability criterion reproduces the intuitive result that the influx of macrophages to phagocytize oxidized molecules should increase when the LDL concentration is increased from equilibrium in a healthy artery. The second criterion loosely says that diffusive forces must dominate inflammatory forces of chemotaxis in a healthy artery, in accordance with the results of [10]. We then numerically verify that these criteria are sufficient. A forthcoming study will include parameter estimation to increase the reliability of the model.