## Abstract

Atherosclerosis is an inflammatory disease. The atherosclerosis process starts when low-density lipoproteins (LDLs) enter the intima of the blood vessel, where they are oxidized (ox-LDLs). The anti-inflammatory response triggers the recruitment of monocytes. Once in the intima, the monocytes are transformed into macrophages and foam cells, leading to the production of inflammatory cytokines and further recruitment of monocytes. This auto-amplified process leads to the formation of an atherosclerotic plaque and, possibly, to its rupture. In this paper we develop two mathematical models based on reaction–diffusion equations in order to explain the inflammatory process. The first model is one-dimensional: it does not consider the intima’s thickness and shows that low ox-LDL concentrations in the intima do not lead to a chronic inflammatory reaction. Intermediate ox-LDL concentrations correspond to a bistable system, which can lead to a travelling wave that can be initiated by certain conditions, such as infection or injury. High ox-LDL concentrations correspond to a monostable system, and even a small perturbation of the non-inflammatory case leads to travelling-wave propagation, which corresponds to a chronic inflammatory response. The second model we suggest is two-dimensional: it represents a reaction–diffusion system in a strip with nonlinear boundary conditions to describe the recruitment of monocytes as a function of the cytokines’ concentration. We prove the existence of travelling waves and confirm our previous results, which show that atherosclerosis develops as a reaction–diffusion wave. The results of the two models are confirmed by numerical simulations. The latter show that the two-dimensional model converges to the one-dimensional one if the thickness of the intima tends to zero.

## 1. Atherogenesis process

Several theories have been developed to explain the pathogenesis of atherosclerosis, but none of these theories can explain the whole process because of the large number of risk factors involved (Fan & Watanabe 2003). Despite these multiple theories, the concept that is now accepted is that atherosclerosis is an inflammatory disease (Ross 1999).

The initiation of atherosclerosis begins with the entry of low-density lipoproteins (LDLs) into the intima of the blood vessel (Ross 1999; Fan & Watanabe 2003; Østerud & Bjørklid 2003), where they are oxidized (ox-LDLs). The latter is considered as a dangerous agent, hence an anti-inflammatory reaction is launched: monocytes adhere to the endothelium, then they penetrate into the intima, where they are transformed into macrophages.

The macrophages phagocytose the ox-LDLs but this eventually transforms them into foam cells, which in turn have to be removed by the immune system. Hence, they set up a chronic inflammatory reaction (auto-amplification phenomenon) by secreting pro-inflammatory cytokines (tumour necrosis factor TNF-a, interleukin IL-1) that promote the recruitment of new monocytes and the production of new pro-inflammatory cytokines.

This auto-amplification phenomenon is compensated by two anti-inflammatory reactions: a ‘biochemical’ one mediated by the anti-inflammatory cytokines (IL-10), which inhibit the production of pro-inflammatory cytokines; and a ‘mechanical’ one when the smooth muscle cells migrate and proliferate to create a fibrous cap over the lipid deposit, which isolates this deposit centre from the blood flow.

The atherosclerotic plaque formed by the fibrous cap and the lipid deposit changes the geometry of the vessel and interacts with the blood flow. This interaction may lead to a thrombus, or to the degradation and rupture of the plaque (Li *et al.*2006*a*,*b*).

In this paper we consider only the setting up of the inflammatory reaction with its biochemical and mechanical inhibitions. The first model is one-dimensional. It is a simplified model and captures some essential features of atherosclerosis development, but it does not take into account the finite width of the blood vessel wall. This approximation implies that the vessel wall is very narrow and the concentrations across it are practically constant.

The second model is two-dimensional. It will allow us to describe in more detail the recruitment of monocytes from the blood flow. The flux of monocytes depends on the concentration of cytokines at the surface of endothelial cells that separate the blood flow and the intima. This should be described by nonlinear boundary conditions that change the mathematical nature of the problem. We will study it in this paper.

## 2. One-dimensional model

We consider a reaction–diffusion system of equations on an interval representing the intima:
2.1
for *x*∈[0,*L*]. Here *M* is the density of immune cells (monocytes, macrophages) and *A* is the density of the cytokines secreted by the immune cells. The function *f*_{1}(*A*)=(*α*_{1}+*β*_{1}*A*)/(1+*A*/*τ*_{1}) models the recruitment of the immune cells from the blood flow, promoted by the inflammatory cytokines, with *α*_{1}=*f*_{1}(0) corresponding to the beginning of the inflammation, that is, the recruitment of monocytes due to the presence of ox-LDLs, promoted by endothelial adhesion molecules, chemoattractants and growth factors. The factor *β*_{1} represents the auto-amplification of the recruitment of monocytes due to the inflammatory cytokines secreted by the monocytes themselves. The factor 1+*A*/*τ*_{1} represents the mechanical saturation of the recruitment of *M*, with *τ*_{1} being the characteristic time for the fibrous cap formation. The term *f*_{2}(*A*)*M* models the cytokine production rate, with *f*_{2}(*A*)=*α*_{2}*A*/(1+*A*/*τ*_{2}), where *α*_{2}*A* represents the secretion of pro-inflammatory cytokines promoted by the pro-inflammatory cytokines themselves and 1+*A*/*τ*_{2} represents the inhibition of the pro-inflammatory cytokines’ secretion mediated by the anti-inflammatory cytokines, with *τ*_{2} being the necessary time for this inhibition to act. The terms −*λ*_{1}*M* and −*λ*_{2}*A* represent the degradation of the immune cells *M* and the cytokines *A*, respectively, and *d*_{1}(∂^{2}*M*/∂*x*^{2}) and *d*_{2}(∂^{2}*A*/∂*x*^{2}) describe their diffusion (or cell displacement) in the intima.

All the parameters of the model, *α*_{1}, *β*_{1}, *τ*_{1}, *α*_{2}, *τ*_{2}, *λ*_{1}, *λ*_{2}, *d*_{1} and *d*_{2}, are assumed to be non-negative. Besides, for *f*_{1} to be an increasing function of *A*, we impose the condition that
2.2

### (a) Kinetic system

In order to determine the conditions for the setting up of the inflammatory reaction, we first study only the reaction part of system (2.1), that is, we consider the kinetic system of equations
2.3
The point *E*_{0}=(0,*α*_{1}/*λ*_{1}) is an equilibrium point for any value of the parameter. The existence of other equilibrium points is determined by the equation
2.4
Three cases are possible:

There is only one positive solution of equation (2.4), denoted by

*A*_{r}, with*E*_{r}=(*A*_{r},*M*_{r}) its corresponding equilibrium.There exist two positive solutions of equation (2.4), denoted by

*A*_{l}and*A*_{r}, where*A*_{l}<*A*_{r}. The corresponding equilibrium points are denoted by*E*_{l}=(*A*_{l},*M*_{l}) and*E*_{r}=(*A*_{r},*M*_{r}).There is no positive solution to equation (2.4).

Denoting *E*_{1}=(0,*λ*_{2}/*α*_{2}), we observe that, if *E*_{0} is situated below *E*_{1}, then *E*_{0} is a stable equilibrium, otherwise it is unstable, *E*_{l} is unstable and *E*_{r} is a stable equilibrium.

### (b) Biological interpretation

The equilibrium points of the kinetic system admit the following biological interpretation: *E*_{0} (no cytokines and low concentration of immune cells) corresponds to the non-inflammatory state, whereas *E*_{r} (large concentrations of cytokines and immune cells) corresponds to the inflammatory state; *E*_{0} can be stable or unstable, and *E*_{r} is always stable when it exists.

The intermediate equilibrium *E*_{l} is always unstable when it exists. It represents a threshold that the system has to overcome in order to move from the non-inflammatory state *E*_{0} to the inflammatory state *E*_{r}.

Hence, the biological interpretation can be given in terms of the parameter *α*_{1} (ox-LDL concentration):

If

*α*_{1}is small (low ox-LDL concentration),*E*_{0}is the only equilibrium and it is stable. No chronic inflammatory reaction can be set up.If

*α*_{1}is intermediate, there are three equilibrium points:*E*_{0}and*E*_{r}are stable and*E*_{l}is unstable. This is called the*bistable case*. The system will reach*E*_{r}if the initial conditions are large enough and*E*_{0}otherwise. Hence, a chronic inflammatory reaction may be set up, but for that it has to overcome a threshold.If

*α*_{1}is large, there are two equilibrium points:*E*_{0}is unstable and*E*_{r}is stable. This is called the*monostable case*. Even a small perturbation of*E*_{0}will lead to*E*_{r}. Hence, even a small perturbation of the non-inflammatory state leads to the setting up of a chronic inflammatory reaction.

### (c) Existence of travelling waves

In this subsection we study the conditions for the *propagation* of the chronic inflammatory reaction. To this aim we turn back to system (2.1). For the theoretical study the space domain will be the entire real line and for the numerical simulations it will be the segment [0,*L*]. A travelling-wave solution of system (2.1) is a particular solution representing a front with constant velocity connecting two equilibria of system (2.3), that is, a solution of the form *W*(*x*−*ct*), where the constant *c* is the speed of the wave and , with *W*_{+} and *W*_{−} being two equilibria of system (2.3).

### Theorem 2.1.

*In the bistable case, there exists a unique travelling-wave solution connecting the non-inflammatory state E*_{0} *and the inflammatory state E*_{r}*, i.e. a constant c and a vector-valued function W(x−ct), the solution of system (2.1) on the real line* *and satisfying W*_{−}*=E*_{r} *and W*_{+}*=E*_{0}*.*

*In the monostable case, there exists a constant c** *such that, for all* *, there exists a travelling-wave solution of velocity c connecting the non-inflammatory state E*_{0} *and the inflammatory state E*_{r}*, i.e. a vector-valued function W(x−ct), the solution of system (2.1) on the real line* *and satisfying W*_{−}*=E*_{r} *and W*_{+}*=E*_{0}*.*

The results on the existence of waves are confirmed by numerical simulations (El Khatib *et al.* 2007).

## 3. Two-dimensional model

### (a) Mathematical model

We consider the system of equations 3.1 and 3.2 in the two-dimensional strip, , with the boundary conditions 3.3 and the initial conditions 3.4

The functions *f* and *g* are sufficiently smooth and satisfy the following conditions:
and *g*′(*A*)>0. We put *A*_{0}=*b*/*γ*. This is a constant level of cytokines in the intima such that the corresponding concentration of the monocytes is zero, and they are not recruited through the boundary. The constant vector (*M*,*A*)=(*M*_{0},*A*_{0}), where *M*_{0}=0 is a stationary solution of problem (3.1)–(3.3).

The function *f*(*A*) represents the flux of monocytes into the intima through the endothelium of the blood vessel; it depends on the concentration of cytokines.

### (b) Stationary solutions in the interval

Consider the problem in the perpendicular section of the strip:
3.5
3.6
3.7
Here a prime denotes the derivative with respect to *y*. It has a constant stationary solution
We linearize (3.5)–(3.7) about this solution and consider the corresponding eigenvalue problem:
3.8
3.9
3.10
We consider the case *λ*=0. From (3.8) and (3.9) we get
where , and

From the boundary conditions at *y*=0 and the boundary at *y*=*h* we get
and
3.11
where
Solutions of equation (3.11) give zero eigenvalues of problem (3.8)–(3.10).

### Proposition 3.1.

*Suppose that μ*_{i}≠0, *σ*_{i}≠0, *i*=1,2, *and σ*_{1}≠*σ*_{2}*. For all h sufficiently small, the principal eigenvalue of problem (3.8)–(3.10) is in the right half plane. If f(A*_{0}*) or g′(A*_{0}*) are sufficiently small and h sufficiently large, then the principal eigenvalue is in the left half plane.*

### Proposition 3.2.

*If the principal eigenvalue of problem (3.8)–(3.10) crosses the origin from negative to positive values, then the stationary solution M*=0, *A=A*_{0} *of problem (3.5)–(3.7) becomes unstable and two other stable stationary solutions bifurcate from it. For one of these solutions, M*_{s}*(y),A*_{s}*(y), the inequality*
3.12
*holds.*

The existence and stability of a bifurcating solution follows from the standard arguments related to the topological degree (Mawhin 1979). Inequality (3.12) follows from the positiveness of the eigenfunction corresponding to the zero eigenvalue. The proofs of these propositions and the assertions below will be presented elsewhere.

### (c) Existence of waves in the monostable case

In this section we consider problem (3.1)–(3.3) assuming that the stationary solution (*M*,*A*)=(*M*_{0},*A*_{0})=(0,*A*_{0}) is unstable and that there exists a stable stationary solution *M*_{s}(*y*),*A*_{s}(*y*) in the section of the cylinder such that
Here, we will study the existence of waves with the limits (*M*_{0},*A*_{0}) at and (*M*_{s},*A*_{s}) at . We assume that there are no other stationary solutions *M*(*y*), *A*(*y*) to system (3.8) and (3.9), such that
3.13

Consider the problem
3.14
3.15
3.16
Here *c* is the wave velocity. We will look for a solution (*M*(*x*,*y*,*t*),*A*(*x*,*y*,*t*)) such that, for all *y*∈[0,*h*] and all *t*≥0:
3.17

Let *μ*(*x*,*y*) and *α*(*x*,*y*) be some functions continuous together with their second derivatives and such that
3.18
3.19
Denote

### Proposition 3.3.

*Let functions μ(x,y) and α(x,y) satisfy conditions (3.18) and (3.19). If*
3.20
*then there exists a solution of problem (3.14)–(3.17).*

The main result of this section is given by the following theorem.

### Theorem 3.4.

*Problem (3.14)–(3.17) has a solution if and only if c satisfies the inequality*
*where the infimum is taken with respect to all functions satisfying conditions (3.18) and (3.19). These solutions are strictly monotone with respect to x.*

## 4. Numerical simulations

In this section we present numerical simulations of problem (3.1)–(3.3) in the bounded domain Ω_{s}=(*x*,*y*), 0≤*x*≤1, 0≤*y*≤1, with the additional boundary conditions at the sides of the rectangle:
The functions *f*(*A*) and *g*(*A*) are taken in the form

We carry out the simulations using the software COMSOL Multiphysics.

In the approximation of a thin domain, for such functions *f* and *g* we obtain the one-dimensional system (2.1) in the monostable case. Therefore, we can expect the monostable behaviour in the two-dimensional case. This means the absence of the threshold where even small perturbations of the disease-free solution lead to disease development. In this case, concentrations *A* and *M* grow and spread in the form of travelling waves. Figure 1*a* presents the propagation of the travelling wave in both the one-dimensional and two-dimensional models. The comparison shows a good agreement between these two cases when the strip thickness is small. Figure 1*b* demonstrates how the speed of propagation in the two-dimensional case depends on the strip thickness. The speed of the two-dimensional wave converges, numerically, to the speed of the one-dimensional wave as the width of the domain goes to zero. Figure 2 shows the wave propagation for wide domains. In this case the wave is essentially two-dimensional. To watch videos of the simulations, please see the electronic supplementary material.

## 5. Discussion

Atherosclerosis and other inflammatory diseases develop as a self-accelerating process that can be described using the reaction–diffusion equations. In the first part of this paper, we have developed a one-dimensional model for the early stage of atherosclerosis. The model is applicable for the case of a small thickness of the intima (blood vessel wall), which corresponds to biological reality. We prove the existence of a travelling-wave solution of the reaction–diffusion system and explain the chronic inflammatory reaction as the propagation of a travelling wave.

In the second part, we take into account the thickness of the intima by developing a two-dimensional model in which the recruitment of monocytes is described by a nonlinear boundary condition. The latter is a function of the cytokines’ concentration in the intima. We proved the existence of a travelling-wave solution in the monostable case and explain, as in the one-dimensional system, the inflammatory reaction as a wave propagation. Then we verified these theoretical results with numerical simulations. We also verified numerically the limiting passage from the two-dimensional model to the one-dimensional one.

Further development of atherosclerosis results in remodelling of the vessel. This means that the lumen (the channel where the blood flow takes place) can retract and the vessel wall takes a specific bell shape. This can essentially modify the characteristics of the flow, and mechanical interaction of the flow with the vessel walls become crucial because it can result in the rupture of the plaque. There are numerous studies of these phenomena (e.g. Li *et al.*2006*a*,*b*). The blood flow influences the development of the plaque: the shear stress activates the receptors of the endothelial cells and accelerates the recruitment of monocytes.

Another important question is related to risk factors like hypercholesterolaemia, diabetes or hypertension. They determine some parameters of the mathematical model. A more complete description would consist in supposing that this influence increases slowly during life. The parameters of the model would then evolve slowly, and the system would pass from the disease-free case to the bistable state to finally reach the monostable state. In each state, the ignition itself would be due to an accidental disturbance, such as an injury that can initiate infection.

Another approach to modelling atherosclerosis is based on cellular automata (Poston & Poston 2007). The authors investigate ‘the hypothesis that plaque is the result of self-perpetuating propagating process driven by macrophages’. The macrophage recruitment rate is considered as a steeply rising function of the number of macrophages locally present in the intima. The main result of Poston & Poston (2007) confirms the conclusion of this work that atherosclerosis development can be viewed as a wave propagation.

## Footnotes

One contribution of 17 to a Theme Issue ‘From biological and clinical experiments to mathematical models’.

- © 2009 The Royal Society