Illuminating Atherogenesis Through Mathematical Modeling

Mathematical Medicine is a relatively new and expanding area of Applied Mathematics research with a growing number of mathematicians, experimentalist, biomedical engineers, and research physicians involved in collaborative efforts on a global scale. Mathematical models are playing an increasing role in our understanding of such complex biological processes as the onset, progression, and mitigation of various diseases. The cardiovascular system is particularly intricate, and the formulation and analysis of mathematical models presents a myriad of challenges to the investigator. (See (Quarteroni, 2001) for a survey on the subject.) Mathematical studies of the cardiovascular system have included continuum mechanical models of vascular soft tissue (Holzapfel et al., 2000; Humphery & Rajagopal, 2002; Taylor & Humphrey, 2009), fluid dynamical models of the interaction between blood flow and vessel walls (Baek et al., 2007; Quarteroni, 2001; Veneziani, 1998), and mathematical models, such as that of the present work, of biochemical characteristics of the vasculature (Ibragimov et al., 2005; Neumann et al., 1990; Saidel et al., 1987). The disease of atherosclerosis, and its initiation atherogenesis, involves a complex interplay between mechanical, genetic, pathogenic, and biochemical processes. A comprehensive view of atherosclerosis will ultimately require integration of these various modeling perspectives. Herein, we focus on the inflammatory component of atherogenesis, in particular the role of immune cells–primarily macrophages–in the presence of oxidatively modified low density lipoproteins (LDL cholesterol) within the intimal layer of large muscular arteries. We present a mathematical model of the key inflammatory spiral that characterizes the initiation of atherosclerosis, and perform some analyses of this model. It is well accepted that atherosclerosis is marked by chronic inflammation (Creager & Braunwald , eds.; Fan & Watanabe, 2003; Ross, 1995; 1999; Wilson , ed.). Changes in the permeability of the endothelial layer and subsequent deposition of lipids in the intima cause an up-regulation of chemoattractants such as monocyte chemotactic protein 1, interleukin-8 and macrophage colony-stimulating factor that are secreted by the endothelial and other cells. In addition, LDL molecules become trapped in the subendothelial intima where they 3


Introduction
Mathematical Medicine is a relatively new and expanding area of Applied Mathematics research with a growing number of mathematicians, experimentalist, biomedical engineers, and research physicians involved in collaborative efforts on a global scale. Mathematical models are playing an increasing role in our understanding of such complex biological processes as the onset, progression, and mitigation of various diseases. The cardiovascular system is particularly intricate, and the formulation and analysis of mathematical models presents a myriad of challenges to the investigator. (See (Quarteroni, 2001) for a survey on the subject.) Mathematical studies of the cardiovascular system have included continuum mechanical models of vascular soft tissue (Holzapfel et al., 2000;Humphery & Rajagopal, 2002;Taylor & Humphrey, 2009), fluid dynamical models of the interaction between blood flow and vessel walls (Baek et al., 2007;Quarteroni, 2001;Veneziani, 1998), and mathematical models, such as that of the present work, of biochemical characteristics of the vasculature (Ibragimov et al., 2005;Neumann et al., 1990;Saidel et al., 1987). The disease of atherosclerosis, and its initiation atherogenesis, involves a complex interplay between mechanical, genetic, pathogenic, and biochemical processes. A comprehensive view of atherosclerosis will ultimately require integration of these various modeling perspectives. Herein, we focus on the inflammatory component of atherogenesis, in particular the role of immune cells-primarily macrophages-in the presence of oxidatively modified low density lipoproteins (LDL cholesterol) within the intimal layer of large muscular arteries. We present a mathematical model of the key inflammatory spiral that characterizes the initiation of atherosclerosis, and perform some analyses of this model. It is well accepted that atherosclerosis is marked by chronic inflammation (Creager & Braunwald , eds.;Fan & Watanabe, 2003;Ross, 1995;Wilson , ed.). Changes in the permeability of the endothelial layer and subsequent deposition of lipids in the intima cause an up-regulation of chemoattractants such as monocyte chemotactic protein 1, interleukin-8 and macrophage colony-stimulating factor that are secreted by the endothelial and other cells. In addition, LDL molecules become trapped in the subendothelial intima where they 3 www.intechopen.com are subject to oxidative modification by reactive oxygen species. Macrophages begin to accumulate in the region where they assume a pro-inflammatory phenotype. The stimulation of macrophages may be due to the presence of inflammatory extra-cellular matrix fragments. In addition, oxidized LDL is recognized by a macrophage scavenger receptor with some degree of interindividual variation (Boullier et al., 2001;Martín-Fuentes et al., 2007;Mosser & Edwards, 2008;Podrez et al., 2002). Macrophages attempting to internalize these particles may become engorged with cholesterol and transform into foam cells. In this state, these immune cells are incapable of performing the customary immune function and become part of a developing atherosclerotic lesion. The immune response is mediated by those chemical signals emitted by endothelial cell, immune cells, and immune cell derived foam cells. The corruption 1 of the immune process caused by ingestion of oxidized LDL can trigger an inflammatory response which results in increased immune cell migration to the site, possible further corruption, and ultimately accumulation of debris (necrotic, apoptotic, and lipid laden cells) characterizing plaque onset. This inflammatory spiral facilitated by chemotaxis, the process modeled herein, is a hallmark of atherogenesis. It will become evident that our model incorporates many parameters characterizing such things as the rate at which macrophages move within the intimal tissue (independent of and in response to chemokines), rates of phenotype changes for macrophages, rates of phagocytosis and uptake of lipids by immune cells, degratation rates of various chemicals, chemical reaction rates and so forth. Some of these have value ranges that are known in vitro or in vivo, but many are unknown. The analytical techniques employed at present are linear stability studies. This allows us to obtain criteria based on the relative values of parameters and to interpret these criteria in terms of the propensity for a lesion to initiate-or not. These criteria will take the form of various inequalities in section 4. In the next section, we lay out the disease paradigm and the assumptions upon which the mathematical model is constructed. This is followed by a presentation of the general model in the form of a system of nonlinear, primarily parabolic partial differential equations with mixed third type boundary conditions. In section 4, we perform stability analyses of the model under two different assumptions regarding the source of inflammatory components. Two stability theorems are given along with a bio-medical interpretation of the criteria derived. Also included is a discussion of the existence of unstable equilibria with a focus on the role of an antioxidant presence and the competing processes of macrophage motility (unrelated to chemotaxis) and chemotaxis. The chapter closes with a brief conclusion.

The disease paradigm and model basis
The large muscular arteries most vulnerable to atherosclerotic lesions can be considered as thick walled tubes consisting of three distinct layers. The outermost layer, called the adventitia, provides structural integrity through a strong collagen network. The middle layer, the media, provides flexibility and adaptability through layers of smooth muscle cells enmeshed in an elastin and collagen network. And the thin, innermost layer, called the intima, is where the atherosclerotic lesions begin to develop. A monolayer of endothelial cells forms an interface between the intima and the lumen through which blood flows. These endothelial cells are highly active in the circulatory process providing a smooth surface for fluid flow, secreting anticoagulant, procoagulant, and inflammatory factors and regulating the exchange of cells and molecules between the blood and arterial wall. Insult to the endothelium and so called endothelial dysfunction is a precursor to atherosclerosis. A number of pathological conditions, genetic factors, and behavioral practices may result in endothelial dysfunction (Davignon & Ganz, 2004). This process appears to be characterized by a change in the permeability of the endothelial layer that allows lipids to migrate into the subendothelial layer followed by an influx of the cells that comprise the immune response. Following endothelial dysfunction and migration of lipoproteins and immune cells, we identify two significant steps to atherogenesis. These are oxidative modification of LDL, and the initiation of an inflammatory spiral.

Lipoprotein oxidation
Lipoproteins are micellar particles which contain regulatory proteins that direct the blood trafficking of cholesterol and other lipids to various cells in the body.. There are four major classes of lipoprotein-chylomicrons, very low-density lipoproteins (VLDL), low-density lipoproteins, and high-density lipoproteins (HDL)-but the bulk of cholesterol is contained in the latter two. Low density lipoproteins consist of a lipid core, a surface protein and a number of antioxidant defenses. LDLs deposit cholesterol in the tissues for cell metabolism. High density lipoproteins contain most of the remaining cholesterol in the body. These particles take excess lipids from tissues and return them to the liver for processing-the process referred to as reverse transport. Elevated plasma levels of LDL indicate a high risk of disease primarily because of their susceptibility to becoming trapped within the intima and subsequently attacked by radical oxygen species. The inflammation of atherosclerotic lesions occurs in areas of intimal thickening enriched by deposits of oxidized LDL. The modification of LDL is a complicated process that has been the subject of several studies, and the reader is directed to the articles (Parthasarathy et al., 1992;Steinberg, 1997) and the review (Young & McEneny, 2001) and the references therein for a more complete and detailed description. Cobbold, Sherratt and Maxwell provided a mathematical model of the in vitro cascade of oxidation of LDL cholesterol in 2002 formulated according to a linear chemical reaction process (Cobbold et al., 2002). This model is adapted and included in the present model of atherogenesis. In brief, the mechanics of the process can be described as follows: In the tissue, where the concentration of reactive oxygen species (ROS) may be relatively high and external antioxidant defenses low, each interaction of an ROS and an LDL molecule will result in oxidizing one of the vitamin E molecules on the lipid surface. It is also possible that an oxidized vitamin E molecule (α-tocopherol radical) may be reduced back to a vitamin E molecule by an antioxidant present (Niki et al., 1984;Watanabe et al., 1999). If, through a finite sequence of oxidation of vitamin E molecules, an LDL molecule losses all of its innate defense against free radical attack, it is susceptible to peroxidation of its lipid core. Once fully modified, the oxidized LDL is both attractive to macrophages and unable to leave the intima (unlike oxidized HDL particles (Tall, 1998)).

The inflammatory response
Accompanying the permeability changes to the endothelium and the influx of lipids is the immune response. Various white blood cells (monocytes, T-lymphocytes, neutrophils) migrate into tissues in response to chemical signals. Once in the subendothelial intima, monocytes differentiate into macrophages. Under normal healthy conditions, these immune cells aid in the degradation of apoptotic cells as well as the removal of foreign agents such 51 Illuminating Atherogenesis Through Mathematical Modeling www.intechopen.com as bacteria or viruses through phagocytosis. As stated, macrophages have a high affinity for oxidized LDL. However, attempts to take up the modified lipids by the process of phagocytosis are unsuccessful, and the lipid laden macrophages transform into foam cells (Goldstein et al., 1977;Podrez et al., 2002;Steinberg, 1997). Unable to perform their normal immune function, these lipid-laden cells signal other immune cells to the site precipitating accumulation of fatty tissue and the progression toward plaque growth. Additional chemical signals secreted by the foam cells and endothelial cells summon more immune cells to the site. Additional macrophages migrate to the localized site of inflammation. The chemical mediators of inflammation can increase binding of oxidized LDL to cells in the arterial wall (Hajjar & Haberland, 1997). Hence, the new macrophages become engorged with oxidized LDL and the cycle of chemical signaling continues. The role of macrophages in initiation of an atherosclerotic lesion is complicated and far from singular. 2 In addition to foam cells, apoptotic macrophages are regularly found in lesions. Apoptosis of cells (macrophage and others) within a plaque is found to have both stabilizing as well as destabilizing affects (Cui et al., 2007;Tabas, 2004). Phagocytosis of apoptotic cells (not necessarily macrophages) may induce resistance to foam cell formation among macrophages. This occurs when during phagocytosis, the macrophage takes in high levels of membrane-derived cholesterol as opposed to lipoprotein-derived cholesterol. In Cui et al. , the authors report that ingestion of apoptotic cells induced a survival response in the macrophages in their experiments (Cui et al., 2007). It is also known that macrophages appear in different phenotypes that are non-static in the sense that they may change types-a process that is reversible (Kadl et al., 2010;Mosser & Edwards, 2008;Stout et al., 2005), and that the different types serve opposing functions (e.g. inflammatory versus anti-inflammatory). Moreover, the sources of additional immune cells include transport across the endothelium as well as migration via the vasa vasorum that provides blood to the artery wall. The mathematical model is constructed to allow for the diverse functions of the immune cells. (For a mathematical study similar to that presented here that focuses primarily on the competing role of inflammatory and anti-inflammatory macrophages, we direct you to the article (Ibragimov et al., 2008).)

The mathematical model
We begin by identifying the key chemical and cellular species involved in atherogenesis. For each species, an evolution equation is derived through the classical approach of imposing a mass balance in an arbitrary control volume and subsequent reduction to a pointwise statement. We do not consider here the volume of a lesion but rather the concentration of each species at any point. Our model consists of five classes of generalized species-two cellular and three chemical-that have critical roles in the initiation of an atherosclerotic lesion. These classes are labeled and denoted as follows: I Immune cells: These are primarily monocyte derived macrophages but may include other white blood cells (T-cells and perhaps neutrophils).
D Debris: This is the bulk of a forming lesion consisting of apoptotic cells, macrophage derived foam cells, and potentially necrotic tissue . Our use of the term debris is unconventional in the sense that we do not intend to suggest that these are inert cells or simply a byproduct of some process and merely occupy space. As will be seen in the mathematics to follow, this species type plays a pivotal role in the inflammatory feedback.
C Chemoattractant: This chemical species represents any of a number of cytokines and chemotactic molecules including macrophage colony stimulating factor, monocyte chemotactic protein-1, and various interleukin proteins. Any chemical that is used in the regulation of the immune response primarily through inducing chemotaxis is included in this species type.
L & L ox Low density lipoproteins: The LDL species consists of two major sub-types, those in a native (un-oxidized) state, and those molecules that have undergone full peroxidation of the lipid core.
R Reactive oxygen species: These are free radical molecules that induce oxidative damage to the lipoproteins present. This species is a byproduct of various metabolic processes within the arterial wall.
Also included in the model are several input parameters. Of particular interest are the parameters A ox , that is a level of antioxidants such as vitamins C, E, and beta-carotene, and L B representing the serum concentration of LDL. Each of the representative variables here is a vector with each component representing a specific member of the class. For example, ..,N I may be a different specific white blood cell, a different phenotype, or may represent cells in different roles. We allow for I to have N I components, D to have N D , C to have N C , L to have N L + 1 3 , and each of L ox and R are scalar valued. If we isolate any representative variable u from this list, we construct an equation of the form that equates the evolution of the concentration of species u to a spatial flux field J u and any net source Q u due to cellular interactions, chemical secretion or uptake, chemical reactions, and the like. The flux fields and source terms are outline below for each variable.

Governing equations
The equations governing these species and based upon the disease paradigm outlined in section 2 are The model of LDL oxidation presented by Cobbold, Sherratt and Maxwell includes LDL molecules in a fully native state containing N L vitamin E molecules. Studies show the number of such antioxidant defenses is on average 6 per LDL molecule but may vary from 3 to 15 (Esterbauer et al., 1992;Stocker, 1999). We then consider L i to contain i vitamin E molecules where L 0 represents LDL molecules completely depleted of native vitamin E that has yet to undergo full oxidation of its core.
The various parameters appearing in (3.1)-(3.8) require explanation; a succint description of each is given in table 1. Each species is subject to diffusion, or diffusive motility in the case of immune cells, and this is reflected in the flux terms μ u ∇ 2 u (u represents any of the various state variables I-R) with the coefficient μ with a subscript a measure of the motility or diffusive capability of the respective species.
χ ik chemotactic sensitivity of immune species i to chemical stimulant k a ik ,â ik binding of immune cells to the lesion for removal b ik measure of subspecies interaction for immune cells c i , f i rates of foam cell formation d ni cell turn over or chemical degradation rate p ik rate of chemical attractant production due to the lesion presence e ik uptake of chemoattractant during chemotaxis k R , k R 0 , k A rate of oxidation, peroxidation, and reverse (anti-oxidation), respectively d ni cell turn over or chemical degradation rate p ik rate of chemical attractant production due to the lesion presence e ik uptake of chemoattractant during chemotaxis τ i , h efficiency factors p R production of free-radicals due to normal metabolism Table 1. Bio-physiological Interpretation of Parameters The terms χ ik (C k , I i )∇C k are the contribution to the flux field for macrophages due to chemotaxis. The coefficient χ ik (C k , I i ) is the chemo-tactic sensitivity of immune cell i to chemoattractant k. This is the classic Keller-Segal model of chemotaxis (Keller & Segel, 1971). The dependence of χ ik on the immune cells is generally taken to be linear, however there is no present need to specify a particular form for these functions. Each of the immune cells, debris, chemoattractants, and native LDL species may undergo natural turnover or chemical degradation represented by the last terms in equations (3.1)-(3.6). The immune cell equations contain three significant cross interaction terms. The terms a ik D k I i capture binding of macrophages with debris-in particular, these and the analogous termŝ a ik I k D i in (3.2), account for phagocytosis of debris by healthy macrophages and removal for future processing in the liver. We also allow for inter-species interactions via the terms b ik I k I i in (3.1). This accounts, for example, for potential change of phenotype of macrophages during the inflammatory process. Recent studies have demonstrated that such changes may occur (reversibly) in vitro and in animal models (Kadl et al., 2010;Stout et al., 2005). Finally, the formation of foam cells through binding with oxidized LDL appears in equations (3.1) and (3.7) in the removal terms c i I i L ox and f k I k L ox . This foam cell formation appears as a source term in equation ( The parameter τ i allows us to catagorize different contributions to the lesion-different types of debris. The equation for the chemoattractants includes (in addition to those terms already mentioned) a source term reflecting production of these chemicals in response to the presence of debris p ik D k C i . The removal terms e ik I k C i represent the reduction of the chemoattractant concentration by binding with macrophages during chemotaxis. As stated, the equations governing the lipid oxidation reactions (3.4)-(3.8) are a modification of the model of lipoprotein oxidation presented by Cobbold, Sherratt andMaxwell in 2002 (Cobbold et al., 2002). The chemical kinetics are assumed to be a linear reaction model in which an LDL molecule containing i vitamin E particles reacts with a reactive oxygen species, with reaction rate k R to produce an LDL molecule with i − 1 vitamin E molecules. This model also allows for the reverse oxidation reaction in that an LDL molecule with i < N L vitamin E molecules may react with the antioxidant species A ox , with reaction rate k A ,t o produce an LDL molecule with i + 1 vitamin E defenses. Any LDL molecule that has been completely depleted of its native antioxidant defenses contributes to the concentration L 0 . A subsequent reaction of an L 0 molecule with an ROS (with reaction rate k R 0 ) results in peroxidation of the lipid core and a fully modified LDL particle. The ROS is depleted through these reactions and through direct reaction with the anti-oxidant species-the latter occurring with the rate of reaction h appearing in equation (3.8). The primary source of ROS is as a byproduct of metabolic processes within the intima. The term p R represents this source. The reader is encouraged to see (Cobbold et al., 2002) for a detailed construction of the model. The modifications of Cobbold-Sherratt-Maxwell model presented here are two fold. First, we allow for spatial variation through a standard Fickian diffusion. More significant to the study of atherogenesis, we include the uptake of modified LDL by macrophages leading to foam cell formation and subsequent inflammation. The terms f k I k L ox represent removal of oxidized LDL through macrophage binding and the contribution to the forming lesion as seen in (3.2), (3.7) and (3.8).

Domain and boundary conditions
The system (3.1)-(3.8) can be considered in two or three spatial dimensions. In (Ibragimov et al., 2005), the current authors performed numerical simulations of a simplified model accounting only for immune cells, debris, chemoattractant, and later smooth muscle cells (a species not considered here as we are interested only with the earliest onset of cellular aggregation). Such simulations were performed in a two dimensional annular domain, and demonstrated the ability of the model to produce such features as localization of immune cells during inflammation and localized aggregation. The subsequent focus has been on illuminating the interplay of the various parameters by considering the initiation of inflammation as due to an instability in an equilibrium state. The general spatial regime considered is a deformed annulus (in two dimensions) or a deformed annular tube (in three spatial dimensions). In either case, the mathematical domain Ω is intended to represent the tunica intima, the innermost subendothelial layer of an arterial wall. The annulus, or annular tube, has an inner and outer boundary denoted by Γ I and Γ O , respectively. The inner boundary Γ I corresponds to the monolayer of endothelial cells that form the interface between the arterial wall and the lumen, while the boundary Γ O represents the inner elastic lamina that separates the intima from the media. In the following analyses, we will assume that there is no transport of any species across the boundary Γ O . While there may well be some transport across this elastic lamina-in particular of free radicals due to metabolic processes within the media-we will assume here that any such contribution is negligible relative to production, consumption, and inter-species reactions within the intima. Influx through the inner boundary Γ I is for some species a significant source in the model. In particular, the chemoattractant and native LDL are subject to a third type boundary condition on Γ I modeling transport in response to a chemical potential across the endothelial cells. This corresponds mathematically to the conditions Here, n is the outward unit normal to Γ I , C * is a baseline level of chemoattractant present at the endothelium, and L B is the serum level of LDL. The parameters α C,i are assumed to be non-negative. However, the sign of α L is not specified so that (3.10) may correspond with either forward transport of native LDL into the subendothelial intima or reverse transport of native LDL into the blood. We assume here that LDL in the blood stream is fully native (has undergone no free radical attack) so that only native LDL is capable of either forward or reverse transport. The immune cells are also subject to transport across the endothelium. The mechanism here is a chemo-tactic sensitivity regulated by the level of chemoattractant at the endothelium. The boundary condition is therefore a mixed third type condition with the flux of immune cells dependent on the chemoattractant species.
Each function α I,i (C), i = 1,...,N I is a nonnegative monotone function of the vector C of chemoattractants 4 . The remaining boundary conditions are This is the mathematical representation of the previous statement that no transport of any species across the inner elastic lamina separating the intima and the media is considered significant relative to the interactions within the intima, and that only fully native LDL, 4 We can state the boundary condition for the immune cells in the more general form ∇C k is the flux field for the i th immune cell species, andᾱ I,i is a corresponding reformulation of the right hand side of (3.11).
immune cells and chemoattractant enter into the system via the endothelial layer. We may further consider the completely homogeneous Neumann conditions under the conditions that α C,i = α L = α I,i = 0. This closed system requires a modification of (3.1) to include a source term. This may be interpreted as modeling the vasa vasorum as the sole source of immune cells contributing to the inflammatory process. In reality, supply both via the vasa vasorum and via transport across the endothelium occur simultaneously. Study of the two extreme cases considered here is done to illuminate both the biological and mathematical differences these two delivery mechanisms make in the modeling and analysis.

Mathematical analysis of the model
There are several approaches to analyzing a particular mathematical model including numerical simulations, asymptotic and perturbation methods, and stability analyses. As suggested, the last of these, stability analyses, is particularly applicable under the present circumstances since we do not have experimental data from which to glean relevant ranges for many of the parameters. A classical approach to mathematical models of biological phenomena-especially those characterized by pattern formation, morphogenesis, and aggregation (Keller & Segel, 1971;Murray, 2002;Turing, 1952), is to consider significant state changes as resulting from a mathematical instability. This will result in the criteria based on relative parameter ranges. The inequalities will depend not only on the relationships between parameter ranges, but also on the source of inflammatory factors, and on the size of the domain (intimal thickness). We present stability analyses of the system (3.1)-(3.8) under some specified conditions. The system considered throughout this section will be simplified to account for one of each of the species types I, D, C, one native LDL species (which may be considered an averaging over each of L i ), an oxidized LDL species, and free-radicals. The system of equations is The modification to (3.1) appearing in (4.1) includes the source term of macrophages via the vasa vasorum as previously indicated (which may be set to zero if appropriate.) Since we are considering only one native LDL species, we also modify the equations to allow for reverse oxidation of oxidized LDL and allow for an efficiency factor r for such reactions. Subscripts have been eliminated where they are no longer needed. For ease of notationĉ = c + f . Our analysis of (4.1)-(4.6) consists of a linear stability analysis using an energy estimate-i.e. Lyapunov functional--approach. That is, we consider certain equilibrium solutions of this system as characterizing a healthy state free from certain inflammatory markers. We then ask whether such equilibria are linearly, asymptotically stable.

Stability with zero transport across the endothelium
We consider a uniform, healthy equilibrium solution of (4.1)-(4.6) subject to the boundary conditions (3.9), (3.10), (3.11), and (3.12) in the special case that α C = α L = α I = 0. We label this equilibrium solution (I e , D e , C e , L e , L oxe , R e ), and introduce the perturbation variables u, v, w, z, y, s which are defined by I = I e + u, D = D e + v, C = C e + w, L = L e + z, L ox = L oxe + y, and R = R e + s.
Substituting the assumed form for I-R into (4.1)-(4.6) and keeping only terms that are linear in the perturbation variables results in the system of equations with the boundary conditions (4.13) The various parameters appearing here are the rates at equilibrium given by A = d 1 , B = cL oxe , C = aD e , D = aI e , E = cI e , F = c 15 L oxe , G =âD e , H =âI e , I = d 2 , J =ĉI e , K = eC e , L = p, M = eI e , N = d 3 , P 1 = k R R e + d 4 , P 2 = k A A ox r, and χ = χ(I e , C e ). Each of these constants is assumed to be nonnegative, and due to balance of mass F = B + Q 1 , J = E + Q 4 , Q 2 =(P 1 − d 4 )+R 1 , and Q 5 = P 3 + R 2 . Let U =(u, v, w, z, y, s). Before proceeding, we define stability in the following way: Definition 4.1. The equilibrium state is called asymptotically stable if every solution of the linearized initial boundary value problem (4.7)-(4.13) for the perturbation variables vanishes at infinity in the sense that there exists a positive functional Our study of (4.7)-(4.13) requires the construction of an appropriate functional F , and this construction gives rise to the inequalities involving the parameters including intimal thickness. In the interest of brevity, much of the computational details are omitted here. The main results are stated with a discussion. We begin by assuming that the product terms uv and uw are nonnegative within Ω. Physically, this can be interpreted as saying that an increase in debris (v > 0) and an increase in chemoattractant (w > 0) results in an increase in immune cells (u > 0). Likewise a decrease in debris and chemoattractant (v < 0, w < 0) is met with a decrease in immune cells (u < 0). This is a rather minor and biologically reasonable condition. However it can be dropped, and a weaker stability theorem obtained (Ibragimov et al., 2010a). The transition matrix characterizing the species interactions associated with the system (4.7)-(4.12) is We will assume that the eigenvalues of Λ have negative real part. (The implication of this and other imposed conditions will be discussed later.) In the following construction, this ensures that integrals of the form Ω U i → 0a st → ∞ for U i = u, v, w, z, y,o rs. This follows from Green's theorem and the homogeneous Neumann boundary conditions. This constraint does not guarantee stability of the system or even point-wise boundedness of each U i . We will also assume here that μ D = 0 which is consistent with the immobile nature of the lesion core. A sequence of inequalities is obtained by multiplying (4.7) by u (4.8) by v, and so forth and integrating over the domain Ω to secure bounds on the rate of change of the total energy of the perturbations. In so doing, we introduce consideration of the geometry and size of the domain through use of the Poincaré inequality Here, |Ω| is the volume of the domain, and the parameter C p is dependent on the geometry of the domain 5 . For ease of notation, we set And in addition to the condition imposed upon the matrix Λ, suppose that Following a systematic construction of integral inequalities from the equations (4.7)-(4.12) we arrive at the principal inequality essential to the present analysis. (4.14) The coefficients on the right hand side of the inequality (4.14) are (4.15) We are now able to state our first major result.
The proof requires a definition of the functional as the obvious modification of the left hand side of (4.14). Of interest are the physical interpretations of the sufficiency conditions stated here. The meaning of the conditions on the products uv and uw has already been given. It can also be noted that each of the coefficients appearing in the array (4.15) is written as a positive term minus a non-negative term to highlight the relationships necessary between the parameters to guarantee stability. The condition on the matrix Λ-that its eigenvalues have negative real part-has distinct bio-medical interpretation. Parameters E, J, and Q 1 are the rates of foam cell production by binding of macrophages to oxidized LDL. If these are large, then they are a source to the lesion. If each of these is small (conditions 4.1.1 and 4.1.6), then to leading order Λ is block diagonal. The parameters Q 2 and R 1 are the oxidation rates of LDL. If Q 2 << 1 and R 1 << 1-so that C z , C s > 0, then the eigenvalues of the lower 3 × 3 block has negative eigenvalues −P 1 , −(Q 3 + Q 4 ), and −(R 2 + R 3 ). A healthy system would be dominated by the antioxidant reactions which correspond to large values of P 1 , Q 3 , Q 4 , R 2 ,a n dR 3 . If in addition L << 1 (condition 4.1.4), then the production of chemoattractant due to the presence of the lesion is small, and the eigenvalues of the upper 3 × 3 block are to leading order Large M 1 indicates a fast degradation of chemoattractant and sufficient uptake of chemokines by macrophages to minimize immune cell migration (reduce inflammation.) The parameters D and G 1 are large (condition 4.1.5) when the uncorrupted, healthy immune function dominates through normal phagocytosis of lesion debris. Similarly, large values for A 1 (due to dependence on A and C) and large H 1 correspond to degradation of the lesion and clearing by macrophages. The sufficient condition for stability is the inequality Several of the requirements for stability rest on the interplay between chemotactic effects and diffusion/motility. This is typical of systems characterized by chemotaxis. Conditions 4.1.2 and 4.1.3 as well as the positivity of each parameter in the array (4.15) provide a minimal requirement of the diffusivity of the intimal layer and the motility of macrophages-motility unrelated to chemotaxis-to guarantee that a perturbation off of the healthy equilibrium state decays.

Stability with transport of macrophages and LDL across the endothelium
We again consider the simplified system (4.1)-(4.6) and perturb off of a healthy equilibrium solutions (I e ,...,R e ). However, we consider the boundary conditions (3.9), (3.10), (3.11), and (3.12) with α C > 0, α L = 0, and the form of α I appearing in (3.11) as where C * is a base line serum level of chemoattractant and α 0 I is a positive constant. If the level of chemotaxis inducing agents at the endothelial interface is greater than an average level in the blood stream, then macrophages (or monocytes which differentiate) will enter into the subendothelial intima. The perturbation variables, u,...,s are defined in the same manner as in 4.1, and the linearized system (4.7)-(4.12) is again studied. However, the boundary conditions on Γ I for the variables u, w, and z (corresponding to immune cells, chemoattractant, and native LDL, respectively) in the present analysis are nonhomogeneous and must be derived from (3.9), (3.10), and (3.11). It should be noted that the existence of a spatially uniform equilibrium requires C e = C * and L e = L B with C * and L B the serum levels of chemoattractant and native LDL introduced in section 3.2. From (3.11) The additional boundary conditions on both Γ I and on Γ O remain as in section 4.1.
The approach applied previously must be modified here to account for the effect of the boundary terms on the total energy of each perturbation variable. In addition to the Poincaré inequality, we require the well known Sobolev trace and generalized Friedrich's inequalities (Sobolev Trace) ∂Ω u 2 ds ≤ C 1 Ω u 2 + |∇u| 2 dx, and We note that these inequalities also depend on the geometry of the domain through the constants C 1 , C 2 , and C 3 . The consideration of best estimates for these constants has received much attention (Acosta & Durán, 2003;Mazya, 1985). For the tubular domain considered herein, the present authors provide estimates of the constants appearing in each of these inequalities (including the Poincaré inequality) in (Ibragimov et al., 2010b). The procedure is similar to that used in 4.1, and we obtain a set of inequalities relating the parameters of the system that provide sufficient conditions under which the perturbations will decay. Again, much of the computational details are omitted (the interested reader is referred to (Ibragimov et al., 2010a) and to (Ibragimov et al., 2008;2010b) for similar results). Instead, we highlight a number of inequalities in light of the bio-medical significance and state the primary result.
To facilitate the analysis, we assume that the decrease of oxidized LDL due to attempted phagocytosis by macrophages is negligible compared the increase and decreases resulting from the chemical reactions with free-radical and antioxidant species. (This is to say that uptake by macrophages is a minor effect on the oxidized LDL concentration, not that foam cell formation is negligible especially as it relates to debris growth or decay.) This is equivalent to the previous case where Q 1 is small and corresponds to f = 0 so that c =ĉ, Q 1 = Q 4 = 0.
Here, we no longer consider the transition matrix because we cannot impose any physically reasonable constraints to guarantee that the integrals Ω y 2 and Ω s 2 can be ignored throughout the construction. (The terms 1 |Ω| Ω y and 1 |Ω| Ω s are the average values of the total perturbations of oxidized LDL and free radicals, respectively, over the entire domain.) Instead, these are treated in the same manner as each of the perturbation variables. The competing effects of diffusion (cellular motility) and chemotaxis are prominent in the result in section 4.1, and the same is true when considering boundary transport. In the present case, however, the sufficient conditions require the diffusion to overcome both chemotaxis within the intima as well as that across the endothelial layer. In particular, stability will rest on the conditions The latter condition arose in the previous result, however the former relates the impact of chemotaxis at the boundary through the parameters α 0 I and α C on the net motility of macrophages within the intima. A direct comparison with C ∇u appearing in (4.15) reveals the additional requirement on this motility to overcome chemotactic effects when boundary transport is accounted for. The diffusion and degradation of chemoattractant are also required to be significantly increased in this case. Set C(ᾱ,μ C )=min(ᾱ C 3 ,μ C ) wherē (C 3 is the constant from the Friedrich inequality.) The function C(ᾱ,μ C ) is nondecreasing in α andμ C independently, and will increase if both of these increase. The condition C w > 0in theorem 1 will be replaced by Recall that M 1 is the rate at which the chemoattractant is reduced within the intima due to natural degradation and through uptake by macrophages during chemotaxis. The added source of chemoattractant from the endothelial boundary is reflected in the new requirement on the size of this parameter.
Of particular interest are the two cases of LDL transport-forward and reverse-that can be admitted by allowing α L to be either positive or negative. When α L < 0, LDL enters into the intima through the endothelial layer. Stability in this case will require The second in condition 4.2.3(-) gives a specific requirement on the removal rate of LDL (P 1 ), especially due to chemical degradation, relative to influx across the endothelial layer (α L ), the oxidation kinetics within the intima ( P 2 +P 3 +R 1 +Q 2 2 ), and the size of the intima (|Ω|). This particular inequality indicates that intimal thickening is destabilizing mathematically. The role of diffuse intimal thickening (DIT) as a precursor to, and in the early stages of, atherosclerosis has been the subject of a number of studies (Nakashima et al., 2008). Those arteries that are prone to atherosclerotic lesions such as the abdominal aorta, carotid, and coronary arteries are observed to express DIT whereas arteries known to be resistant to atherosclerosis do not (Nakashima et al., 2002). Accumulation of oxidized LDL relative to native LDL in the deep region of DIT in human coronary arteries has been observed (Fukuchi et al., 2002). The stability requirement on the degradation of LDL in the case of reverse transport is significantly weaker provided the rate of diffusion of LDL and the rate at which LDL leaves the intima through the endothelial boundary are sufficiently high. Let Then φ 0 is a measure of the total oxidation rate of LDL and depends on the thickness of the intima. If LDL is transported from the intima back to the blood stream, then stability will require If conditions 4.2.1, 4.2.2, and the appropriate version of 4.2.3 (-or +) hold, and we follow the techniques used in section 4.1, we obtain our primary inequality (4.19) The coefficients appearing on the right hand side are defined by (4.20) and we define the parameters The primary result of the current anaylsis is The proof involves the pair of functionals ( u, v, w, z, y, s, Ω y, Ω s) and M = min{C z , C s , C y , C s }. The hypotheses of theorem 4.2 ensure establishing asymptotic stability for this case.

Instability of the equilibrium solution
The theorems obtained in sections 4.1 and 4.2 establish sufficient conditions under which the uniform healthy state is guaranteed to be stable to small perturbations. It is not readily clear whether the inequalities derived are tight-in the sense that they are nearly necessary conditions. One can ask the degree to which these conditions must be violated to result in the existence of a perturbation that will blow up. The existence of perturbations that blow up, in particular spatially nonhomogeneous perturbations, is typically addressed through construction of an explicit example. The classical approach is adopt the ansatz for the perturbation variables Expressing the perturbation as U as in definition 4.1, this gives U(x, t)=e σtŪ (x), and the system (4.7)-(4.12) can be written in the vector/matrix formulation as σŪ = ∇·(M e ∇Ū)+ΛŪ.
( 4.22) The diffusion-chemotaxis coefficient matrix M e has the diffusion coefficients on the main diagonal, χ(I e , C e ) in the first row and third column, and zeroes everywhere else. When considering the case without boundary transport, using the fact that M e and Λ are constant matrices, the ansatz can be further refined to seek solutions of the form Here ξ is a constant vector and φ λ an eigenfunction of the Laplacian on Ω, −∇ 2 φ λ = λφ λ , subject to the completely homogeneous Neumann boundary conditions (4.13). The system (4.7)-(4.12) reduces to the algebraic equation in σ, λ, and ξ σ ξ =(Λ − λM e ) ξ. (4.23) Solutions to (4.23) for which the real part of σ is positive will grow; this is the classic Turing instability problem (Turing, 1952). For the Turing stability problem, one considers the case for which Λ has only negative eigenvalues and M s e , the symmetric part of M e , has only positive eigenvalues. (This latter condition will hold if and only if χ(I e , C e ) < 2 √ μ I μ C , and this inequality follows from the conditions C ∇u > 0 and C ∇w > 0 appearing in theorem 4.1.) For any domain Ω, the first eigenpair is λ 0 = 0 and φ 0 = constant, and it is well known that there is an enumerable set of positive eigenvalues 0 < λ 1 < λ 2 < ··· and corresponding eigenfunctions {φ λ n } that form an orthonormal basis for L 2 (Ω). In the case that Λ has only negative eigenvalues, an instability must come from one of the larger eigenvalues. Unfortunately, finding these eigenvalues explicitly for a general domain is not possible. For an annulus (R 2 ), or an annular cylinder (R 3 ), they can be found by separation of variables. The corresponding eigenmodes in these cases are nonaxisymmetric suggesting that lesion initiation should also be nonaxisymmetric-this is consistent with clinical observations. For the present case with no transport of any species across the endothelial layer, we can study the effect of the antioxidant level on the stability of the healthy state. If Mφ 0 in equation (4.1) is replaced with Mφ 0 C (to make the source explicitly dependent on the chemoattractant), or if Mφ 0 = 0, then the equilibrium solution is (I e , D e , C e , L e , L oxe , R e )=( 0, 0, 0, L e , L oxe , R e ).
We can, after some lengthy calculations, show that in the limiting cases as A ox → 0 + and A ox → ∞ (Ibragimov et al., 2010b) L oxe ∝ A −1 ox , and R e ∝ p R ,a sA ox → 0 + and L oxe ∝ A −2 ox , and R e ∝ A −1 ox ,a sA ox → ∞. The latter result demonstrates that the antioxidants strongly control free radical production and LDL oxidation. When studying the spectrum of Λ − λM e , this asymptotic result shows that for A ox sufficiently large, the lower 3 × 3 block corresponding to the lipid chemistry produces no eigenvalues with positive real part. The question of most interest is what conditions are required for stability (or produce an instability) for the full system in the case that the lipid chemistry alone is stable. If (L e , L oxe , R e ) is a stable equilibrium for the lipid equations in isolation, what effect does the antioxidant level have on the stability of the equilibrium (0, 0, 0, L e , L oxe , R e )? In the limit as the antioxidant level vanishes, the equilibrium will be unstable whenever (4.24) The critical and competing roles of diffusion and chemotaxis are prominent in this criterion providing an unstable equilibrium. For any positive eigenvalue λ and diffusive capacity of the chemoattractant μ C , if the chemotactic sensitivity χ(0, 0) is large enough, the perturbation will grow away from the healthy equilibrium to some other state. An analysis like the above can be employed with any variation of the system (3.1)-(3.8) provided the boundary conditions considered are completely homogeneous of Neumann type. For example, in (Ibragimov et al., 2008) the present authors consider a system characterized by two distinct macrophage phenotypes each subject to diffusion, chemotaxis, the potential to change phenotypes, but for which only one subspecies was subject to foam cell formation. We showed that the stability result provided therein-that analogous to theorem 4.1 here-was strongly dependent on the dominance of diffusion over chemotaxis. As is seen here, for any set of other parameter values, the chemotactic sensitivity coefficient can always be taken sufficiently large to produce an unstable equilibrium. The question of unstable equilibria for the case with boundary transport can likewise be considered. Not surprisingly this presents a far more complicated situation mathematically. Even if we only consider constant equilibria, the special approach based on the ansatz (4.21) and the spectral analysis above does not yield any instability examples. Moreover, the coupled boundary conditions provide a Laplacian that is not self-adjoint and does not allow us the option to expand all of the perturbation variables using any single family of eigenfunctions. Nevertheless, a careful construction within the appropriate mathematical framework will provide conditions for which an unstable equilibrium exists. The effect of antioxidant concentration on stability can be analyzed. The reader is encouraged to see (Ibragimov et al., 2010b) for a construction in this case.

Conclusion
The purpose of this work is twofold. We have formulated a mathematical model of the inflammatory process that characterizes atherogenesis. This model given in equations (3.1)-(3.8) is presented in general terms to provide a framework for the ongoing modeling process. With this in mind, adaptations are easily included as our understanding of this complex medical process increases. We believe that mathematical modeling provides a useful tool to meet the goals of medical research on atherogenesis-identifying vulnerability to disease, development of treatments, and promotion of preventative interventions. Computer simulation (in silico analysis) requires a model consistent with and able to capture the characteristics of disease as observed in vivo.
Here, we have also studied the model by performing stability analyses under two different assumptions regarding the supply of inflammatory components-macrophages, chemotactic chemokines, and LDL. Taking the vasa vasorum as the sole source of these species, we arrive at a distinct set of inequality conditions on the system parameters that will guarantee that perturbations off of the healthy equilibrium state will decay. Bio-medically, the perturbations are interpreted as the start of inflammation, and the starting equilibrium as a disease-free state. A stable equilibrium is then seen as representing a cellular configuration that is robust-where a lesion is unlikely to develop in the short term. An unstable one suggests that (bio-chemically) the location is vulnerable to atherogenesis and the potential for development of a fibrofatty lesion or latter fibrous plaque. In addition to the positive stability criteria obtained using the energy estimate, we offer a negative result in the form of construction of an instability example. This latter condition highlights the inflammation mitigating effects of antioxidant presence and the significant interplay between chemotaxis and diffusion when the antioxidant level becomes negligible. We also raise this same stability question under the assumption that the supply of inflammatory components is from influx from the blood flow via the endothelial interface. We again produce several inequalities that when satisfied by the system parameters ensure that the equilibrium solution is linearly asymptotically stable to small perturbations. Of particular interest in this latter case is the stabilizing effect of reverse transport of native LDL from the intima back to the blood stream. That reverse transport of LDL is stabilizing is not surprising given the corruptive nature of oxidized LDL on macrophage function. Our finding further supports the development of treatment modalities aimed at not only reducing serum LDL levels but at facilitating reverse transport of cholesterol (Superko, 2006). Although the conditions are numerous, clinical values of the various parameters can be easily compared in light of the various inequalities derived and presented in theorems 4.1 and 4.2. The availability of clinical values for several parameters is lacking. Moreover, the parameters appearing in table 1 need not be constant, and determination of appropriate functional forms is an important and difficult task. This will require a process of "fine tuning" through collaboration with clinicians and experimental scientists.