STABILITY ANALYSIS USING AN ENERGY ESTIMATE APPROACH OF A REACTION-DIFFUSION MODEL OF ATHEROGENESIS

Abstract. This paper considers modeling the initiation of atherosclerosis, as an inflammatory instability. Motivated by the disease paradigm articulated by Russell Ross, atherogenesis is viewed as an inflammatory spiral with positive feedback loop involving key cellular and chemical species interacting and reacting within the intimal layer of muscular arteries. The inflammation is modeled through a system of nonlinear reaction/diffusion/convection partial differential equations. The inflammatory spiral is initiated as an instability from a healthy state which is defined to be an equilibrium state devoid of certain key inflammatory markers. Disease initiation is studied through a linear, asymptotic stability analysis of a healthy equilibrium state. Conditions on system parameters guaranteeing stability of the health state and conditions on system parameters leading to instability are given. Among the questions addressed in the analysis is the possible mitigating effect of anti-oxidants upon transition to the inflammatory spiral.

1. Introduction.Mathematical models have a significant role to play in understanding the structure, functioning, evolution and diseases of the cardiovascular system.Moreover, formulating, simulating and analyzing such models offer a vast array of challenges.(See [12] for an interesting survey on the subject.)This article is a continuation of a program to develop, analyze and simulate mathematical models of atherosclerosis initiated by the authors in [7].
Atherosclerosis is a very complex chronic disease of the arterial system with many manifestations and many routes to initiation and progression [3,11,14,13,15].Biochemical, genetic, mechanical and pathogenic factors conspire to initiate and promote the disease.The focus of [9] and the present contribution is the role played by inflammation in atherogenesis [14,4].This is not to suggest that genetic, mechanical and pathogenic factors are unimportant or are subordinate to the inflammatory processes considered herein and in [7] and [9].Account is taken of them through parameters in the equations studied below that model a particular inflammatory cycle thought to play a fundamental role in atherogenesis.In [9] and the present study, attention is focused is upon the inner most layer of the arterial wall, the tunica intima since the very beginning stages of atherosclerosis are largely confined to this layer.
The first extension to the analysis in [7] given by the authors appears in [9] where we used an energy estimate to analyze the stability of a model of atherogenesis that focused on only four species involved in the inflammatory process, and which considered the interplay between viable and apoptotic immune cells.Herein, we consider an extended model that includes the role of low density lipoproteins (LDL) in both a native and chemically modified state (oxLDL), as well as reactive oxygen species (referred to throughout as "free-radicals") present in the subintimal layer.Anti-oxidant effects are also introduced through a parameter.While these species were considered as part of the original model proposed by the authors in [7], they were ignored in the numerical and analytical studies appearing in [7], [8], and [9].
The inflammatory process modeled in [7] involved the following ingredients: two cellular species (smooth muscle cells and macrophages), lesion debris (necrotic cells, lipid core of foam cells, smooth muscle cells) and three molecular species (low density lipoproteins, chemically modified LDL, a chemical signally species).Each of the cellular and molecular species are to be viewed as representative of large classes of cells or molecules exhibiting the functional response attributed to the respective representative.For example, while a number of immune system cells play a role in the inflammatory processes occurring during atherogenesis, the monocyte derived macrophages are probably the dominant players in the creation of the lipid laden foam cells that collect in the lipid rich core of atherosclerotic plaques.Hence to simplify the model, macrophages are the only immune systems cells included in the modeling.Similarly, the LDL species should be viewed as a generic representative of a large class of lipid molecules and oxLDL as a generic representative from the corresponding class of lipids that have been oxidized (chemically modified) by free radicals.
The point of view articulated in [14] and motivating the model adopted in [7] is that atherosclerotic plaques form as a consequence of chronic inflammation sustained through a positive inflammatory feedback loop.The heart of this disease paradigm consists of the following process elements.Through unspecified means, a portion of the endothelial layer of a muscular arterial wall develops a "leaky" spot permitting accelerated transport of LDL (and other macromolecular species) through the endothelial barrier into the intima where they tend to concentrate due to the difficulty of further passage through the inner elastic lamina into the media.Simultaneously, monocytes also enter the intima in response to chemical signaling from an initiating inflammatory reaction (possibly due to viral or bacterial insult, for example).The LDL is eventually chemically modified by free-radicals produced through natural metabolic processes occurring in various cellular species within the arterial wall (e.g.smooth muscle cells, endothelial cells, fibroblasts, etc).Macrophages have an affinity for the oxLDL resulting from this chemical modification process (Indeed, there is strong experimental evidence that macrophages exhibit positive chemotactic sensitivity to these oxLDL species.),eventually becoming foam cells (i.e., macrophages engorged with oxLDL particles).These engorged macrophage derived foam cells are no longer capable of doing their customary job of removing the debris produced by the inflammatory processes; in fact they become components of the growing lesion debris.The growing lesion debris produces various chemical signaling species that attract additional macrophages to the lesion site which then get "corrupted" by the oxLDL species resulting in a chronic inflammatory spiral.
A number of important issues were not addressed in [7] including how to model plaque growth with significant luminal occlusion and how to determine under what conditions the runaway inflammatory/plaque growth spiral occurs and conversely under what conditions the natural defense mechanisms of the body prevent it.The latter question is the subject of [9] and the present paper.
The perspective taken in [9] and extended herein on the latter question is that it is one of stability of the nonlinear reaction-diffusion-chemotaxis system used to model the inflammatory processes initiating atherosclerotic plaque growth.More specifically, the question investigated is whether certain equilibrium states of the governing system of nonlinear partial differential equations, referred to as "healthy states", are linearly, asymptotically stable.These healthy states are characterized by the absence of inflammatory markers, which in the context of the model described above, correspond to equilibrium states in which the macrophage, debris and chemical signal species are at some base-line level in the intimal layer that is commensurate with normal immune function.As stated, the results presented here differ from those obtained in [9] as we account here for LDL, oxLDL, and free-radical interaction and reactions.In addition, we consider herein a more realistic system allowing for boundary transport of some species.The mathematical methods employed in [9] are adapted to account for the increased mathematical complexity introduced.
2. Mathematical Model.The model for atherogenesis of interest here tracks the evolution of six generic "species" which are major contributors to the initial stages of atherosclerosis.These species are generic in that they are representative of classes of factors contributing to the inflammatory processes leading to disease initiation.In this spirit, these representative species are given the labels: immune cells (principally macrophages), debris (developing lesion), chemo-attractant, native LDL, oxidized LDL and free radicals, and denoted Î , D , Ĉ , L, Lox , and R, respectively.
The governing equations for this simplified model are of the form: Here, div and ∇ denote the usual divergence and gradient operators.The various terms appearing on the right side of these equations require some discussion.The term −χ( Î , Ĉ )∇ Ĉ in equation ( 1) is the portion of immune cell flux due to chemotaxis, and the coefficient χ( Î , Ĉ ) is the chemotactic sensitivity.In the classical Keller-Segal model for dictyostelium discoideum, for example, χ( Î , Ĉ ) has the form χ( Î , Ĉ ) = Î / Ĉ [10].At present, there is no need to specify a particular form for χ( Î , Ĉ ).The term d 11 Î represents natural turnover of immune cells.The two terms a 15 Î Lox and a 12 Î D in (1) give the rate at which the macrophage population is diminished through foam cell formation ( through binding with oxidized LDL), and through normal immune function.The latter, for example, could be accounted for by viable macrophages binding with debris for eventual processing in the liver. 1inally, in the stability analyses that follow, we will be considering a perturbation off of a constant level of macrophages.In essence, we are looking at a small time window.The term Mφ 0 in (1) represents a base line level of immune cells present.In general, M φ 0 could depend on the level of chemo-attractant, especially at the boundary where transport across the endothelial layer can occur.We can assume that over the time scales of interest, the value is constant.Mass transport through the endothelial layer will be considered.
The term b 15 Î Lox appearing in equation ( 5) represents conversion of oxidized LDL into foam cells.The balance of mass is captured by c 15 Î Lox which appears in equation ( 2); thus we have c 15 = a 15 + b 15 .The term a 21 Î D is the rate at which debris is removed by uncorrupted macrophages while d 22 D is a natural turnover rate for debris.In (3), p 32 D is the rate at which chemo-attractant is produced by the lesion debris, while a 31 Ĉ Î is the rate by which the chemo-attractant concentration is diminished by binding with macrophages.The term d 33 Ĉ is a natural chemical degradation rate for the chemo-attractant.
In (4), (5), and (6), a 46 LR and b 46 LR are the rates at which the native LDL and free radical concentrations are diminished by free radical oxidation of the native LDL (and their sum c 46 = a 46 + b 46 added to the Lox concentration), while A ox r 4 Lox is the rate at which the anti-oxidant concentration, A ox , is able to reverse the oxidative damage done to LDL by the free radicals.The coefficient b 4 (with 0 < b 4 < 1) is an efficiency parameter representing the fraction of the products of the A ox -L ox reverse reaction feeding back into the native LDL population2 .Finally, p R in ( 6) is the rate of free radical production-Free radical production is considered here to be a metabolic byproduct which for present purposes will be assumed to be constantand A ox b 6 R is the rate at which the anti-oxidant concentration is able to reduce the free radical concentration through direct reaction.
We will consider the equations ( 1)-( 6) to hold in a domain Ω with inner and outer boundaries Γ 1 and Γ 2 , respectively.While we will not specify the geometry exactly, Ω can be taken as a deformed annulus in two dimensions, or an annular (deformed) cylinder in three dimensions.The inner boundary Γ 1 represents the endothelial layer which is the interface between the intima and the blood stream.
The outer boundary is the interface between the intima and the muscular media layer.In our analysis, we will assume that all species are subject to homogeneous Neuman boundary conditions on the outer boundary.On Γ 1 , however, we have The parameters α 1 and α 3 are positive constants as is Ĉ * .The value Ĉ * is a baseline level of chemo-attractant in the blood.If the level of chemo-attractant at the endothelial layer is greater than the baseline level, chemo-attractant enters the blood stream while immune cells enter into the subendothelial intima.
3. Linear stability analysis.We begin by assuming that there is a constant equilibrium state ( Îe , De , Ĉe , Le , Loxe , Re ), and introduce the perturbation variables u, v, w, z, y, s which are defined by Substituting the assumed form for Î -R into (1)-( 6) and keeping only terms that are linear in the perturbation variables results in the system of equations ∂z ∂t = div (µ 4 ∇z) − P 1 z + P 2 y − P 3 s (10) For v, z, y, and s ∂v ∂n For u and w we note that Similarly For ease of notation, we have introduced a number of parameters.The new parameters are defined to be: A ox , and χ = χ( Îe , Ĉe ).Each of these constants is assumed to be nonnegative.Note that due to balance of mass In our analysis, we will make the simplifying assumption that the mobility and diffusions coefficients µ i , i = 1, . . ., 6 are constant.
For the stability analysis that follows, we consider classical solutions of the system ( 7)-( 12) subject to the boundary conditions ( 13)- (15).That is, each of u-s is in C 2,1 (x,t) (Ω × (0, ∞)).Note that in the special case w = 0 and Ku = Lv, such a solution can be obtained as an eigenfunction expansion.However, in the more general case, we are interested in the qualitative nature of solutions as perturbations off of the stated equilibrium state.We also note here that the nonlinear convection due to chemotaxis is well known to exhibit blow-up in finite time in domains of dimension 2 or 3 (see for example the reviews [6,5] and the references therein).This was observed numerically in a simplified version of the nonlinear model presented here [7].The standard linearization requires an assumption that perturbations are small enough so that only first order terms are appreciable.As a result the convection appearing in ( 1) is not present in (7).Our interest is in identifying sufficient criteria under which the constant state of the nonlinear system is stable-as defined below-to small disturbances.
Let U = (u, v, w, z, y, s).Before proceeding, we define stability in the following way: Definition 3.1.The equilibrium state is called monotonically stable if every solution of the linearized initial boundary value problem ( 7)-( 15) for the perturbation variables decays in the sense that there exists a positive functional 4. Construction of an appropriate norm.We will construct an inequality that allows us to conclude sufficient conditions under which the equilibrium state is stable.We will make use of the well known Cauchy inequality.Moreover, to account for the nonhomogeneous boundary terms, the following inequalities will be utilized: The coefficients C p , C 1 , C 2 , and C 3 depend on the geometry 3 of the domain Ω, and Γ is the boundary of the domain.
For ease of notation, we set Then we begin our construction by multiplying (7) by u, (8) by v, and so forth and integrating by parts to obtain (all integration that follows is over Ω except where 3 In the case of an annulus of inner and outer radii r I and r O , respectively, Cp is related to the first nonzero eigenvalue of the Laplacian when considering an We assume that the effect of foam cell formation on the concentration of oxLDL is negligible as compared to the competing oxidizing and anti-oxidant reactions.Thus, for simplicity, we set Q 1 = 0 and likewise Q 4 = 0.This of course requires that c 15 = a 15 .In addition, the following condition will be imposed: have negative real part We note these conditions ensure that z, y, and s vanish as t → ∞. Next, we can apply the Cauchy and Sobolev inequalities to (16) to arrive at the inequality 1 2 Let us note that large µ 3 should enhance the stability of the system in general, since this would indicate strong diffusive effects.Similarly, small α 3 and small χ would be associated with stability since this corresponds to weak cumulative (in the domain and on the boundary) chemotactic effects.If µ 3 is larger than χα 3 , we would expect this to be stabilizing.We impose the condition This condition coupled with (22) implies If we sum (18) and ( 23), we find that For ease of notation, we will introduce the parameter ᾱ Similarly, we introduce the parameter μ3 and impose the condition ).Then C(ᾱ, μ3 ) will increase if both ᾱ and μ3 increase, and is at least nondecreasing in ᾱ and μ3 independently.Letting C = C 2 , where C 2 is the other geometrically dependent constant appearing in the Generalized Frederich inequality, we have (after applying said inequality)

Application of the Cauchy and Poincaré inequalities to
Combining ( 24), ( 17), ( 25)-( 27) and applying the Poincaré inequality to the terms involving |∇z|, |∇y|, and |∇s| we arrive at the major inequality for the case under consideration.The details are included in the appendix; here we state the result.Define the functional U 2 i and each of the parameters Obviously, F 1 is positive and vanishes only if U = 0. Then let the functional F 2 be defined by The following conditions ensure stability:  Thus the solution to the linearized system ( 7)-( 12) subject to the boundary conditions ( 13)-( 15) is monotonically stable but not necessarily asymptotically stable.
In order to ensure that F 1 → 0, we would require some additional conditions on the coefficients of the original system.

5.
Discussion.A physical interpretation of the result is in order.Note that if Q 2 and R 1 are very small, then to leading order the eigenvalues of Λ 1 are −P 1 , −Q 3 , and −(R 2 + R 3 ) which are all negative.Now, Q 2 and R 1 are rates of oxidation of LDL, a destabilizing reaction, whereas Q 3 , and (R 2 +R 3 ) are rates of healthy restoration due to anti-oxidant reaction.So, these eigenvalues being negative indicates dominance of anti-oxidant reactions over oxidation of LDL.This is physically realistic as a stability-i.e.indication of healthiness-condition.Conditions 4 and 5 hold if A 1 , H 1 and M 1 are sufficiently large.Large M 1 indicates low levels of the chemo-attractant consistent with low inflammation, while large A 1 (due to A and C) and H 1 indicate healthy immune function since these are rates of decrease of immune cells and debris due to normal immune response.Conditions 4 through 7 also require the diffusion dominates over chemotaxis-this is well known to be stabilizing in systems involving chemotaxis reactions.Together, these conditions provide specific relationships between parameters that can be investigated experimentally as data becomes available.