N El Khatib

S Genieys

B Kazmierczak

V Volpert

Reaction-Diffusion Model of Atherosclerosis Development

Keywords: atherosclerosis, reaction-diffusion waves, nonlinear boundary conditions, existence, numerical simulations AMS subject classification: 35K57, 92C50

Atherosclerosis begins as an inflammation in blood vessels walls (intima). Antiinflammatory response of the organism leads to the recruitment of monocytes. Trapped in the intima, they differentiate into macrophages and foam cells leading to the production of inflammatory cytokines and further recruitment of white blood cells. This self-accelerating process, strongly influenced by low-density lipoproteins (cholesterol), results in a dramatic increase of the width of blood vessel walls, formation of an atherosclerotic plaque and, possibly, of its rupture. We suggest a 2D mathematical model of the initiation and development of atherosclerosis which takes into account the concentration of blood cells inside the intima and of pro-and anti-inflammatory cytokines. The model represents a reaction-diffusion system in a strip with nonlinear boundary conditions which describe the recruitment of monocytes as a function of the concentration of inflammatory cytokines. We prove the existence of travelling waves described by this system and confirm our previous results which show that atherosclerosis develops as a reaction-diffusion wave. The theoretical results are confirmed by the results of numerical simulations.

Introduction 1.Biological background

High plasma concentration of low density lipoprotein (LDL) cholesterol is one of the principal risk factors for atherosclerosis. Its mechanism can be sketched as follows [START_REF] Ross | Atherosclerosis -An inflammatory disease[END_REF][START_REF] Bjarne | Role of Monocytes in Atherogenesis[END_REF]: the process of atherosclerosis begins when LDLs penetrate into the intima of the arterial wall where they are oxidized. Oxidized LDL (ox-LDL) in the intima is considered by the immune system as a dangerous substance, hence an immune response is launched: chemoattractants, which mediate the adhesion of the monocytes to the endothelium and the penetration of the monocytes through the endothelium, are released and endothelial cells are activated. As a result, monocytes circulating in the blood adhere to the endothelium and then they penetrate to the arterial intima. Once in the intima, these monocytes are converted into macrophages.

The macrophages phagocytose the ox-LDL but this eventually transforms them into foam cells (lipid-ladden cells) which in turn have to be removed by the immune system. In the same time they set up a chronic inflammatory reaction (auto-amplification phenomenon): they secrete proinflammatory cytokines (e.g., TNF-α, IL-1) which increase endothelial cells activation, promote the recruitment of new monocytes and support the production of new pro-inflammatory cytokines.

This auto-amplification phenomenon is compensated by an anti-inflammatory phenomenon mediated by the anti-inflammatory cytokines (e.g., IL-10) which inhibit the production of proinflammatory cytokines (biochemical anti-inflammation). Next, the inflammation process involves the proliferation and the migration of smooth muscle cells to create a fibrous cap over the lipid deposit which isolates this deposit center from the blood flow (mechanical anti-inflammation).

This mechanical inhibition of the inflammation may become a part of the disease process. Indeed the fibrous cap changes the geometry of the vasculature and modifies the blood flow. The interaction between the flow and the cap may lead to a thrombus, or to the degradation and rupture of the plaque liberating dangerous solid parts in the flow [START_REF] Li | How Critical Is Fibrous Cap Thickness to Carotid Plaque Stability? A Flow Plaque Interaction Model[END_REF][START_REF] Li | Stress analysis of carotid plaque rupture based on in vivo high resolution MRI[END_REF].

In this study we do not address the fluid-structure interaction between the blood flow and the plaque, and only consider the set up of the chronic inflammatory reaction with its biochemical and mechanical inhibitions. In our previous work, we have developed a simplified model of the reactions arising in the arterial intima [START_REF] Khatib | Atherosclerosis initiation modeled as an inflammatory process[END_REF]. The model represents a reaction-diffusion system in one space dimension:

∂M ∂t = d M ∂ 2 M ∂x 2 + g(A) -βM, (1.1)
∂A ∂t = d A ∂ 2 A ∂x 2 + f (A)M -γA + b, (1.2)
where M is the concentration of monocytes, macrophages and foam cells in the intima, A is the concentration of cytokines. The function g(A) describes the recruitment of monocytes from the blood flow, f (A)M is the rate of production of the cytokines which depends on their concentration and on the concentration of the blood cells. The negative terms correspond to the natural death of the blood cells and of the chemical substances, while the last term in the right-hand side of equation (1.2) describes the ground level of the cytokines in the intima. This model allows us to give the following biological interpretation: at low LDL concentrations the auto-amplification phenomenon does not set up and no chronic inflammatory reaction occurs. At intermediate concentrations a perturbation of the non inflammatory state may lead to the chronic inflammation, but it has to overcome a threshold for that. Otherwise the system returns to the disease free state. At large LDL concentrations, even a small perturbation of the non inflammatory state leads to the chronic inflammatory reaction.

We show that inflammation propagates in the intima as a reaction-diffusion wave. In the case of intermediate LDL concentrations, where a threshold occurs, there are two stable equilibria. One of them is disease free, another one corresponds to the inflammatory state. The travelling waves connects these two states and corresponds to the transition from one to another. The second situation, where the concentration of LDL is high, corresponds to the monostable case where the disease free equilibrium is unstable.

Though this model captures some essential features of atherosclerosis development, it does not take into account a finite width of the blood vessel wall. This approximation signifies that the vessel wall is very narrow and the concentrations across it are practically constant. In a more realistic situation, we should consider a multi-dimensional model and take into account the recruitment of monocytes from the blood flow. The flux of monocytes depends on the concentration of cytokines at the surface of endothelial cells which separate the blood flow and the intima. This should be described by nonlinear boundary conditions which change the mathematical nature of the problem. We will study it in this work. We present the mathematical model in the next section. Section 2 is devoted to positiveness and comparison theorems which appear to be valid for the problem under consideration. We take into account here the particular form of the system and of the boundary conditions. In the general case these results are not applicable. In Section 3, we use them to study the existence of travelling waves in the monostable case. The results concerning wave existence are confirmed by the numerical simulations (Section 4).

Mathematical model

We consider the system of equations

∂M ∂t = d M ∆M -βM, (1.3)
∂A ∂t = d A ∆A + f (A)M -γA + b, (1.4) in the two-dimensional strip Ω ⊂ R 2 , Ω = {(x, y), -∞ < x < ∞, 0 ≤ y ≤ h}
with the boundary conditions

y = 0 : ∂M ∂y = 0, ∂A ∂y = 0, y = h : ∂M ∂y = g(A), ∂A ∂y = 0 (1.5)
and the initial conditions

M (x, y, 0) = M 0 (x, y), A(x, y, 0) = A 0 (x, y). (1.6)
Here M is the concentration of white blood cells inside the intima, A is the concentration of cytokines, d M , d A , β, γ, and b are positive constants, the constant b describes a constant source of the activator in the intima. It can be oxidized LDL coming from the blood. Assuming that its diffusion into the vessel wall is sufficiently fast we can describe it by means of the additional term in the equation and not as a flux through the boundary as in the case of monocytes. The functions f and g are sufficiently smooth and satisfy the following conditions:

f (A) > 0 for A > 0, f (0) = 0, f (A) → f + as A → ∞, g(A) > 0 for A > A 0 , g(A 0 ) = 0, g(A) → g + as A → ∞,
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 values A = A 0 , M = 0 is a stationary solution of problem (1.3)-(1.5). The matching condition for the initial and boundary conditions is satisfied, that is the functions M 0 (x, y) and A 0 (x, y) satisfy (1.5). These conditions provide the existence of a unique solution of problem (1.3)- (1.6) in the space C 2+α,1+α/2 ( Ω), 0 < α < 1 of Hölder continuous functions with respect to x and t (Section 2.1).

Positiveness and comparison of solutions 2.1. Existence of solutions

We begin with the result on global existence of solution of problem (1.3)- (1.6). We note that it is considered in an unbounded domain and has nonlinear boundary conditions. Therefore, we cannot directly apply the classical results for semi-linear parabolic problems. Theorem 2.1. Suppose that f (A) ∈ C 1+α (R), g(A) ∈ C 2+α (R) for some α, 0 < α < 1, the initial condition (M 0 (x, y), A 0 (x, y)) of problem (1.3)-(1.6) belongs to (C 2+α ( Ω)) 2 and satisfies boundary conditions (1.5). Then this problem has a unique global solution (M (x, y, t), A(x, y, t)) with the C 2+α,1+α/2 ( Ω × [0, T ])-norm bounded independently of T .

The proof of this theorem is given in the Appendix. We first prove the existence of solutions in bounded rectangles and then pass to the limit as the length of the rectangle increases. A priori estimates of solutions independent of the length of the rectangles allow us to conclude about the existence of solution in the unbounded domain.

Positiveness for linear problems

Consider the linear parabolic problem

∂u ∂t = d 1 ∂ 2 u ∂y 2 -βu, (2.1
)

∂v ∂t = d 2 ∂ 2 v ∂y 2 + a(y, t)u + b(y, t)v -γv, (2.2)
y = 0 : ∂u ∂y = 0, ∂v ∂y = 0, y = h : ∂u ∂y = c(y, t)v, ∂v ∂y = 0, (2.3)
in the interval 0 ≤ y ≤ h with the initial conditions u(0, y, v(0, y). The coefficients a(y, t), b(y, t) belong to C α,α/2 ( ΩT ), and c(y, t) to C 1+α,(1+α)/2 ( ST ), a(y, t) and c(y, t) are non-negative. We assume the matching conditions between the boundary and the initial conditions. This means that u(y, 0), v(y, 0) satisfies the boundary conditions (2.3). Then there exists a unique solution of this problem, and it is continuous for t ≥ 0, 0 ≤ x ≤ h.

Proposition 2.2. Let the initial condition of problem (2.1)-(2.3) be non-negative functions, u(y, 0) ≥ 0, v(y, 0) ≥ 0.

Then the solution u(y, t), v(y, t) of problem (2.1)-(2.3) is nonnegative for all y and t. If moreover u(y, 0) ≡ 0, v(y, 0) ≡ 0, then the solution is strictly positive.

Proof. Consider first the Dirichlet problem

∂u 1 ∂t = d 1 ∂ 2 u 1 ∂y 2 -βu 1 , (2.4
)

∂v 1 ∂t = d 2 ∂ 2 v 1 ∂y 2 + a(y, t)u 1 + b(y, t)v 1 -γv 1 , (2.5)
y = 0, h : u = v = 0. (2.6)
From the positiveness theorem for monotone systems [START_REF] Volpert | Traveling Wave Solutions of Parabolic Systems[END_REF], it follows that if the initial condition (u 0 1 (y), v 0 1 (y)) is non-negative, then the solution is non-negative. If moreover the initial condition is not identically zero, then the solution is strictly positive. In this case for t > 0 y = 0 :

∂u 1 ∂y > 0, ∂v 1 ∂y > 0, y = h : ∂u 1 ∂y < 0, ∂v 1 ∂y < 0. (2.7)
We compare next the solution (u, v) of problem (2.1)-(2.3) with the solution (u 1 , v 1 ) of problem (2.4)- (2.6). Denote by (u 0 (y), v 0 (y)) the initial condition of this problem. Let

u 0 (y) = u 0 1 (y) + ǫ, v 0 (y) = v 0 1 (y) + ǫ
for some small ǫ > 0. We will prove that the solution of problem (2.1)-(2.3) is greater than the solution of problem (2.4)-(2.6). After that, we can pass to the limit as ǫ → 0. Therefore, we will obtain that the solution of problem (2.1)-( 2.3) with non-negative initial conditions is positive if it is not identically zero. However, the initial condition (u 0 (y), v 0 (y)) introduced above may not satisfy the boundary conditions. In this case, we can introduce a modified initial condition (û 0 (y), v0 (y)) such that it satisfies the boundary condition and

max y |û 0 (y) -u 0 (y)| ≤ ǫ 2 , max y |v 0 (y) -v 0 (y)| ≤ ǫ 2 . Hence û0 (y) > u 0 1 (y), v0 (y) > v 0 1 (y), 0 ≤ y ≤ h. The solution of problem (2.1)-(2.
3) with the initial condition (û 0 (y), v0 (y)) exists and it is continuous for t ≥ 0, 0 ≤ y ≤ h. Therefore u(y, t) > u 1 (y, t), v(y, t) > v 1 (y, t), 0 ≤ y ≤ h at least for some small positive t. Suppose that this inequality holds for 0 < t < t 0 and that it is not valid for t = t 0 . Hence at least one of the following equalities hold:

u(y 1 , t 0 ) = u 1 (y 1 , t 0 ), v(y 1 , t 0 ) = v 1 (y 1 , t 0 ) (2.8)
for some y 0 and y 1 . If these points are inside the interval (0, h), that is

u(0, t 0 ) > 0, u(h, t 0 ) > 0, v(0, t 0 ) > 0, v(h, t 0 ) > 0,
then (2.8) cannot hold because of the comparison theorem. Consequently, one of the functions u(y, t 0 ), v(y, t 0 ) equals zero at the boundary of the interval. Suppose that u(0, t 0 ) = 0. Since ∂u(0, t 0 )/∂y = 0 and u(y, t 0 ) ≥ u 1 (y, t 0 ), then ∂u 1 (0, t 0 )/∂y = 0. This contradicts (2.7).

If u(h, t 0 ) = 0, then ∂u(h, t 0 )/∂y ≤ 0 since u(y, t 0 ) ≥ 0. On the other hand, from the boundary condition (2.3) it follows that ∂u(h, t 0 )/∂y ≥ 0. Hence ∂u(h, t 0 )/∂y = 0. Since u(y, t 0 ) ≥ u 1 (y, t 0 ) ≥ 0 and u 1 (h, t 0 ) = 0, then ∂u 1 (h, t 0 )/∂y = 0. This contradicts (2.7).

Other cases, v(0, t 0 ) = 0, v(h, t 0 ) = 0 can be considered similarly. The proposition is proved.

A similar proposition holds for the two-dimensional problem

∂u ∂t = d 1 ∆u -βu, (2.9
)

∂v ∂t = d 2 ∆v + a(y, t)u + b(y, t)v -γv, (2.10)
y = 0 : ∂u ∂y = 0, ∂v ∂y = 0, y = h : ∂u ∂y = c(y, t)v, ∂v ∂y = 0 (2.11)
in the domain Ω.

Proposition 2.2'. Let the initial condition of problem (2.9)-(2.11) be non-negative functions, u(x, y, 0) ≥ 0, v(x, y, 0) ≥ 0. Then the solution u(y, t), v(y, t) of problem (2.9)-(2.11) is nonnegative for all y and t. If moreover u(x, y, 0) ≡ 0, v(x, y, 0) ≡ 0, then the solution is strictly positive.

Comparison of solutions

Let (M 1 , A 1 ) and (M 2 , A 2 ) be two solutions of problem (1.3)-(1.5) from C 2+α,1+α/2 ( QT ). We assume that all conditions of Section 1.2 are satisfied. Moreover, we will impose some additional conditions in the propositions below. Put

u = M 1 -M 2 , v = A 1 -A 2 . Then ∂u ∂t = d M ∆u -βu, (2.12
)

∂v ∂t = d A ∆v + a(x, t)u + b(x, t)v -γv, (2.13)
y = 0 : ∂u ∂y = 0, ∂v ∂y = 0, y = h : ∂u ∂y = c(x, t)v, ∂v ∂y = 0, (2.14) where a(x, t) = f (A 1 (x, t)), b(x, t) = f (A 1 (x, t)) -f 2 (A(x, t)) A 1 (x, t) -A 2 (x, t) M 2 (x, t) , c(x, t) = g(A 1 (x, t)) -g 2 (A(x, t)) A 1 (x, t) -A 2 (x, t) .
From Proposition 2.2 we immediately obtain the result about comparison of solutions.

Proposition 2.3. In addition to conditions of Section 1.2, suppose that f ′ (A) and g ′′ (A) satisfy the Lipschitz condition. Let (M 1 , A 1 ) and (M 2 , A 2 ) be two solutions of problem (1.3)-(1.5). If

M 1 (x, 0) ≥ M 2 (x, 0), A 1 (x, 0) ≥ A 2 (x, 0), x ∈ Ω,
then the same inequalities are valid for the solutions. If moreover

M 1 (x, 0) ≡ M 2 (x, 0), A 1 (x, 0) ≡ A 2 (x, 0), x ∈ Ω,
then the inequalities for the solutions are strict for t > 0.

Proof. It is sufficient to verify that the coefficients of problem (2.12)-(2.14) satisfy required regularity conditions. We have

f (A 1 (x, t)) -f 2 (A(x, t)) A 1 (x, t) -A 2 (x, t) = 1 0 f ′ (sA 1 (x, t) + (1 -s)A 2 (x, t)) ds.
If the derivative f ′ (A) satisfies the Lipschitz condition, then b(x, t) ∈ C α,α/2 ( QT ). Similarly, if the derivative g ′′ (A) satisfies the Lipschitz condition, then c(x, t) ∈ C 1+α,(1+α)/2 ( ST ). Finally, if M i , A i satisfy the matching conditions, then u and v also satisfy them. The proposition is proved.

Proposition 2.4. Suppose that the initial condition of problem (1.3)-(1.5) is such that

d M ∆M -βM > 0, d A ∆A + f (A)M -γA + b > 0.
Then the solution is strictly increasing with respect to t for each x ∈ Ω.

Proof. Denote u = ∂M ∂t , v = ∂A ∂t .
Differentiating problem (1.3)-(1.5) with respect to t, we obtain

∂u ∂t = d M ∆u -βu, (2.15
)

∂v ∂t = d A ∆v + f (A)u + f ′ (A)M v -γv, (2.16
)

y = 0 : ∂u ∂y = 0, ∂v ∂y = 0, y = h : ∂u ∂y = g ′ (A)v, ∂v ∂y = 0.
(2.17)

Since M, A ∈ C 2+α,1+α/2 ( QT ) and the derivatives f ′ (A), g ′′ (A) satisfy the Lipschitz condition, then f (A), f ′ (A)M ∈ C α,α/2 ( QT ) and g ′ (A) ∈ C 1+α,(1+α)/2 ( QT ). Moreover, the derivatives ∂M/∂t and ∂A/∂t belong to C 2+α,α/2 ( QT ). Hence u, v satisfy the matching conditions and

u(x, y, 0) = d M ∆M -βM > 0, v(x, y, 0) = d A ∆A + f (A)M -γA + b > 0.
Since f (A) ≥ 0, g ′ (A) > 0, then u and v are positive for t > 0. The proposition is proved.

Existence of travelling waves

Stationary solutions in the interval

Consider the problem in the section of the strip:

∂M ∂t = d M M ′′ -βM, (3.1)
∂A ∂t = d A A ′′ + f (A)M -γA + b, (3.2)
y = 0 : M ′ = A ′ = 0, y = h : M ′ = g(A), A ′ = 0, (3.3)
where prime denoted the derivative with respect to y. It has a constant stationary solution

M = 0, A = A 0 .
We linearize (3.1)-( 3.3) about this solution and consider the corresponding eigenvalue problem:

d M M ′′ -βM = λM, (3.4)
d A A ′′ + f (A 0 )M -γA = λA, (3.5)
y = 0 : M ′ = A ′ = 0, y = h : M ′ = g ′ (A 0 )A, A ′ = 0. (3.6)
We consider the case λ = 0. From (3.4),

M (y) = c 1 e σ 1 y + c 2 e -σ 1 y ,
where σ 1 = β/d M . From (3.5),

A(y) = c 3 e σ 2 y + c 4 e -σ 2 y + kc 1 e σ 1 y + kc 2 e -σ 1 y ,

where σ 2 = γ/d A , k = - f (A 0 ) d A σ 2 1 -γ = f (A 0 ) d A (σ 2 2 -σ 2 1 )
.

From the boundary conditions at y = 0:

c 1 σ 1 -c 2 σ 1 = 0, c 3 σ 2 -c 4 σ 2 + k(c 1 σ 1 -c 2 σ 1 ) = 0. Therefore c 1 = c 2 , c 3 = c 4 .
From the boundary condition at y = h,

c 1 σ 1 (e σ 1 h -e -σ 1 h ) = g ′ (A 0 ) c 3 (e σ 2 h + e -σ 2 h ) + kc 1 (e σ 1 h + e -σ 1 h ) , c 3 σ 2 (e σ 2 h -e -σ 2 h ) + kc 1 σ 1 (e σ 1 h -e -σ 1 h ) = 0.
We express c 3 from the second equation and substitute into the first equation:

σ 1 sinh(σ 1 h) -kg ′ (A 0 ) cosh(σ 1 h) = - kσ 1 g ′ (A 0 ) σ 2 cosh(σ 2 h) sinh(σ 2 h) sinh(σ 1 h) or µ 1 coth(σ 1 h) = 1 + µ 2 coth(σ 2 h), (3.7)
where

µ i = kg ′ (A 0 ) σ i = f (A 0 )g ′ (A 0 ) d A σ i (σ 2 2 -σ 2 1 ) , i = 1, 2.
Solutions of equation (3.7) give zero eigenvalues of problem (3.4)- (3.6). Denote the left-hand side of equation (3.7) by s 1 (h). Then it is a decreasing function for h > 0,

s 1 (h) ∼ f (A 0 )g ′ (A 0 ) d A σ 2 1 (σ 2 2 -σ 2 1 ) 1 h , h → 0, s 1 (h) → f (A 0 )g ′ (A 0 ) d A σ 1 (σ 2 2 -σ 2 1 ) , h → ∞. (3.8)
The right-hand side s 2 (h):

s 2 (h) ∼ f (A 0 )g ′ (A 0 ) d A σ 2 2 (σ 2 2 -σ 2 1 ) 1 h , h → 0, s 1 (h) → 1 + f (A 0 )g ′ (A 0 ) d A σ 2 (σ 2 2 -σ 2 1 ) , h → ∞. (3.9) Proposition 3.1. Suppose that µ i = 0, σ i = 0, i = 1, 2, σ 1 = σ 2 .
For all h sufficiently small, the principal eigenvalue of problem (3.4)-(3.6) 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.

Proof. Denote by λ 0 the principal eigenvalue of problem (3.4)- (3.6), that is the eigenvalue with the maximal real part. Clearly, if we increase β and γ by the same value, then λ 0 is decreased by the same value. Therefore, for any h fixed and σ 1 , σ 2 sufficiently large, λ 0 becomes negative (or with the negative real part). On the other hand, by virtue of (3.8), (3.9), in this case, s 2 (h) > s 1 (h). Hence, if this inequality is satisfied, then λ 0 is in the left-half plane.

It can be easily verified that s 1 (h) > s 2 (h) for h small enough and σ 1 = σ 2 . Therefore, when we decrease h, the principal eigenvalue crosses the imaginary axes and passes in the right-half plane.

If f (A 0 ) or g ′ (A 0 ) are sufficiently small, then s 2 (h) > s 1 (h) for h large enough. The proposition is proved. Remark 3.2. From the Krein-Rutman theorem it follows that the principal eigenvalue is simple, real and the corresponding eigenfunction is positive. Contrary to the Dirichlet boundary conditions for which the principal eigenvalue growth with the length of the interval being negative for small h, in the problem under consideration it is positive for small h. It is related to the singular character of this problem as h → 0. 3) becomes unstable and two other stable stationary solutions bifurcate from it. For one of these solutions, M s (y), A s (y), the inequality

M s (y) > 0, A s (y) > A 0 , 0 < y < h (3.10) holds.
The existence and stability of a bifurcating solution follows from the standard arguments related to the topological degree. Inequality (3.10) follows from the positiveness of the eigenfunction corresponding to the zero eigenvalue.

At the end of this section, we will find explicitly stationary solutions of problem (3.1)-(3.3) in the particular case where F (A) = f 0 is a constant. From the first equation we have

M (y) = g(A 1 ) 2σ 1 sinh(σ 1 h) e σ 1 y + e -σ 1 y ,
where A h = A(h). Substituting this expression into the second equation, we find

A(y) = - f 0 g(A h ) 2d A (σ 2 1 -σ 2 2 ) × 1 σ 1 sinh(σ 1 h) e σ 1 y + e -σ 1 y - 1 σ 2 sinh(σ 2 h) e σ 2 y + e -σ 2 y + b d A σ 2 2 .
We obtain the following equation with respect to A h :

A h = - f 0 g(A h ) d A (σ 2 1 -σ 2 2 ) cosh(σ 1 y) σ 1 sinh(σ 1 h) - cosh(σ 2 y) σ 2 sinh(σ 2 h) + b d A σ 2 2 .
The number of its solutions is determined by the function g(A). For small σ 1 and σ 2 , using the approximation exp(±σ i h) ≈ 1 ± σ i h, we obtain approximate equation

A h = f 0 g(A h ) d A hσ 2 1 σ 2 2 + b d A σ 2 2 .
We obtain the same equation if we integrate the equalities

d M M ′′ -βM = 0, d A A ′′ + f 0 M -γA + b = 0
from 0 to h, use the boundary conditions (3.3) and suppose that M and A do not depend on y.

Existence of waves in the monostable case

We consider in this section problem (1.3)-(1.5) assuming that the stationary solution M = M 0 , A = 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 M 0 < M s (y), A 0 < A s (y), 0 ≤ y ≤ h.

We will study here the existence of waves with the limits (M 0 , A 0 ) at x = -∞ and (M s , A s ) at x = +∞. We assume that there are no other stationary solutions such that

M 0 ≤ M (y) ≤ M s (y), A 0 ≤ A(y) ≤ A s (y), 0 ≤ y ≤ h. (3.11)
Consider the problem

d M ∆M -c ∂M ∂x -βM = 0, (3.12
) Here c is the wave velocity. We will look for its solution (M, A) such that Proof. From inequality (3.18) it follows that

d A ∆A -c ∂A ∂x + f (A)M -γA + b = 0, (3.13
x = -∞ : M = M 0 , A = A 0 , x = +∞ : M = M s , A = A s . ( 3
d M ∆µ -c ∂µ ∂x -βµ < 0, (3.19)
d A ∆α -c ∂α ∂x + f (α)µ -γα + b < 0. (3.20) Denote Ω N = {(x, y) : x > -N, 0 ≤ y ≤ h}
and consider the initial-boundary value problem for the system We note that the boundary functions at the left boundary of the cylinder and the initials conditions are the same functions which depend only on the y variable. Their choice depends on N . We suppose that they satisfy the following conditions:

∂M ∂t = d M ∆M -c ∂M ∂x -βM, (3.21)
∂A ∂t = d A ∆A -c ∂A ∂x + f (A)M -γA + b (3.
M 0 ≤ M N (y) ≤ µ(-N, y), A 0 ≤ A N (y) ≤ α(-N, y), (3.26) d M M ′′ N -βM N ≥ 0, (3.27) d A A ′′ N + f (A N )M N -γA N + b ≥ 0, ( 3
∂v ∂t = d A ∆v -c ∂v ∂x + f (A)u + f ′ (A)M v -γv (3.31) y = 0 : ∂u ∂y = 0, ∂v ∂y = 0, y = h : ∂u ∂y = g ′ (A)v, ∂v ∂y = 0. (3.32)
For the solution of problem (3.10)-(3.13), the following estimate holds:

M N (x, y, t) ≥ M b (y), A N (x, y, t) ≥ A b (y), (x, y) ∈ Ω N , t > 0.
Therefore, u(-N, y, t) ≥ 0, v(-N, y, t) ≥ 0, 0 ≤ y ≤ h, t > 0.

(3.33)

If the boundary condition at x = -N was

u(-N, y, t) = v(-N, y, t) = 0, 0 ≤ y ≤ h, t ≥ 0,
then the solution of this problem wouldbe identically zero. Since we have inequalities (3.33) at the boundary, then the solution is non-negative. The lemma is proved.

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

c ≥ c 0 = inf µ,α max(S 1 (µ, α), S 2 (µ, α)),
where the infimum is taken with respect to all functions satisfying conditions (3.16), (3.17). These solutions are strictly monotone with respect to x.

Proof. Existence of a solution (M c , A c ) for c > c 0 follows from Proposition 3.4. Since these solutions are bounded in the C 2+δ ( Ω) norm for some δ ∈ (0, 1) independently of c for c close to c 0 , then we can pass to the limit as c → c 0 and obtain a solution (M c 0 , A c 0 ) for c = c 0 . If there exists a solution (M * , A * ) for some c * < c 0 , then

max(S 1 (M * , A * ), S 2 (M * , A * )) = c * < c 0 = inf µ,α max(S 1 (µ, α), S 2 (µ, α)).
This contradiction proves the theorem.

Numerical simulations

In this section we present numerical simulations of problem (1.3)-(1.5) in the bounded domain Ω s = (x, y), 0 ≤ x ≤ 1, 0 ≤ y ≤ 1 with the additional boundary conditions at the sides of the rectangle:

x = 0, 1 : A = M = 0.
The functions f (A) and g(A) are taken in the form

f (A) = A 1 + 43A/42 , g(A) = ǫ 2 + 8A 1 + A .
We carry out the simulations using the software Comsol c MultiPhysics.

In the approximation of thin domain, for such functions f and g we obtain the one-dimensional system (1.1), (1.2) in the monostable case [START_REF] Khatib | Atherosclerosis initiation modeled as an inflammatory process[END_REF]. Therefore we can expect the monostable behavior in the two-dimensional case. This means the absence of the threshold where even small perturbation of the disease free solution lead to the disease development. In this case, concentrations A and M grow and spread in the form of travelling waves. Figure 1 (left) shows the set up of the wave front. At the fist stage of the development of inflammation, monocytes spread across the domain, that is in the y-direction (Figure 1, left) and then along the intima, that is in the x-direction (Figure 1, right). The wave is essentially two-dimensional. When the domain width is sufficiently large, the wave propagation occurs near the surface where there is an excess of monocytes. Their concentration there becomes high leading to an essentially higher speed of propagation. Their concentration inside the intima remains low.

Figure 2 (left) presents propagation of the travelling wave in both 1D and 2D models. The comparison shows a good agreement between these two cases when the strip thickness is small. The right figure demonstrates how the speed of propagation in the 2D case depends on the strip thickness. The speed of the 2D wave converges to the speed of the 1D wave as the width of the domain goes to zero. We recall that the 2D and 1D models are different. The former takes into account the monocyte recruitment through nonlinear boundary conditions, the latter includes a nonlinear production term in the equation. In some cases, the limiting passage from 2D to 1D as the width goes to zero can be justified [START_REF] Khatib | Atherosclerosis initiation modeled as an inflammatory process[END_REF]. This is not proved for travelling waves. The results presented here show this convergence numerically.

Discussion

Atherosclerosis and other inflammatory diseases develop as a self-accelerating process which can be described with reaction-diffusion equations. In [START_REF] Khatib | Atherosclerosis initiation modeled as an inflammatory process[END_REF] 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 the biological reality. We prove the existence of travelling wave solution of the reaction-diffusion system and explain the chronic inflammatory reaction as propagation of a travelling wave.

During atherosclerosis development, the intima thickness grows and we need to take it into account. In this work we study the two-dimensional case where the second dimension corresponds to the direction across intima. Essential difference with the previous model is not only space dimension but also nonlinear boundary conditions which describe recruitment of monocytes through the epithelial layer of the intima. This is a new class of reaction-diffusion systems for which it appears to be possible to study the existence of travelling waves.

Numerical simulations confirm the analytical results. They show wave propagation and allow us to analyze its speed as a function of the parameters of the model. 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 the specific bell shape. This can essentially modify the characteristics of the flow, and mechanical interaction of the flow with the vessel walls becomes crucial because it can result in the plaque rupture. There are numerous studies of these phenomena (see, e.g. [START_REF] Li | How Critical Is Fibrous Cap Thickness to Carotid Plaque Stability? A Flow Plaque Interaction Model[END_REF][START_REF] Li | Stress analysis of carotid plaque rupture based on in vivo high resolution MRI[END_REF]). 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 hypercholesterolemia, 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 the lifetime. The parameters of the model would evolve then slowly, and the system would pass from the diseasefree case to the bistable state to reach finally the monostable state. In each state, the ignition itself would be due to an accidental disturbance, such as an injury that can initiate infection.

The action of these risk factors, which can be taken into account in the mathematical model, is as follows:

1. The influence of the hypertension: it changes the properties of the blood flow and creates a higher pressure on the vessel wall which can activate the receptors and accelerate the recruitment of monocytes. It can also provoke the plaque rupture, 2. The influence of diabetes II: the monocytes and the platelets can be already activated because of the hyperglycemia. The active state of monocytes increases their recruitment and the active state of platelets can cause spontaneous coagulation (thrombosis),

3. The influence of hypercholesterolemia: the cholesterol level in blood can slowly increase during the lifetime. The parameter α 1 of the model, which shows the level of bad cholesterol in blood vessel walls, becomes time dependent. It increases slowly, and so the system passes from the disease-free state to the bistable state and then to the monostable state. The other risk factors can modify the speed of these transitions.

Influence of the risk factors can be studied in relation to medical treatment. In particular, with statins, which are inhibitors of the low density lipoprotein cholesterol. Recent studies shows a reduction of 28% reduction in LDL-C and 5% increase in high-density lipoprotein cholesterol [START_REF] Larosa | Effect of Statins on Risk of Coronary Disease A Meta-analysis of Randomized Controlled Trials[END_REF]. The inhibition of the LDL by statins "appears to be directly proportional to the degree to which they lower lipids" [START_REF] Larosa | Effect of Statins on Risk of Coronary Disease A Meta-analysis of Randomized Controlled Trials[END_REF]. Its action can be taken into account through the parameters of the mathematical model.

Another approach to modelling atherosclerosis is based on cellular automata [START_REF] Poston | A Typical Atherosclerotic Plaque Morphology Produced in Silico by an Atherogenesis Model Based on Self-Perpetuating Propagating Macrophage Recruitment[END_REF]. 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. Smooth muscle cells dynamics also depend on the macrophage number. Macrophages can die with certain probability resulting in lipid accumulation. During the process, fatty streaks of macrophages set up at random sites, which may progress or regress. Some of them develop into progressive focal lesions, that is advanced pieces of plaque which are macrophage-rich and have a central fibrous cap-like region of smooth muscle cells. The main result of [START_REF] Poston | A Typical Atherosclerotic Plaque Morphology Produced in Silico by an Atherogenesis Model Based on Self-Perpetuating Propagating Macrophage Recruitment[END_REF] confirm the conclusion of this work and of the previous work [START_REF] Khatib | Atherosclerosis initiation modeled as an inflammatory process[END_REF] that atherosclerosis development can be viewed as a wave propagation.

Appendix. Existence of solutions of the evolution problem

We will prove here Theorem 2.1. The idea of the proof is quite standard. We first consider smoothed bounded rectangles and prove the existence of solutions in the bounded domains. We use here a priori estimates of solutions. Since they are independent on the length of the rectangle, we can construct a sequence of uniformly bounded solutions in the increasing domains and choose a convergent subsequence. The limiting function will be a solution of the problem in the unbounded strip.

A priori estimates

In order to obtain a priori estimates of solutions of problem (1.3)-(1.6) we will construct an appropriate supersolution. We consider the stationary equation (1.3) and look for its solution which depends only on the y variable. It satisfies the problem

dM ′′ -βM = 0, M (h) = m > 0, M ′ (h) = a , a > 0, (6.1)
where, for simplicity of notation, we replace d M by d. Its solution has the following form:

M (a; y) = e √ β √ d -a √ d + √ β m + e - √ β √ d + 2 √ β y √ d a √ d + √ β m 2 √ βe √ β y √ d .
We note that for

m = coth( √ β h √ d )a √ d/ β (6.2)
we obtain

M (a; y) = a √ d cosh( √ β √ d y) √ β sinh( √ β √ d h) = m cosh( √ β √ d y) cosh( √ β √ d h) (6.3)
and both of the relations M (y) ≥ 0 and M ′ (0) = 0 are satisfied. Let us take a * = g(∞) + ρ with ρ > 0 arbitrarily small. Let M * (y) := M (a * , y). Obviously, according to (6.3), M (a * , y)) > M (a, y) for a ∈ [0, g(∞)].

(6.4)

∂A ∂t = d A ∆A + f (A)M -γA + b (6.6) (6.5)
and construct the boundary conditions

∂M ∂ν = Ψ(z), ∂A ∂ν = 0, z = (x, y) ∈ ∂R l , (6.7)
in such a way that a) in the limit, as the length increases, we obtain the boundary condition (1.5), b) the function (M * (y), A * (y)) is a supersolution for the auxiliary problems in the bounded domains (A * (y) is defined below). The boundary conditions can be defined in the following way. First of all, we can take them symmetric with respect to x = 0, Ψ(x, y) = Ψ(-x, y). For x ≥ 0 the function Ψ(•) is defined as follows:

z ∈ B 1 \ B 1 * : Ψ(z) = 0, z ∈ B 3 \ B 3 * : Ψ(z) = g(A(z)), z ∈ B 3 * : Ψ(z) = (1 -s(x -(l -2δ)))g(A(z)), z ∈ G 2 ∪ B 2 \ B 2 * : Ψ(z) ≡ 0, z ∈ B 2 * : Ψ(z) = -qs (2yδ/(h/2 -δ)) , z ∈ G 1 : Ψ(z) = -q, z ∈ B 1 * : Ψ(z) = -qs(x -(l -2δ)).
Here s(τ ) is C ∞ function such that s(τ ) ≡ 0 for τ ≤ 0 and s(τ ) ≡ 1 for τ ≥ 1, q > 0 is a constant. Now, let A * denote the solution of the boundary value problem:

d A ∆A + f (A)M * (y) -γA + b = 0, in R l , ∂A ∂ν = 0 on ∂R l . (6.8)
It is obvious that, given M * (•), we can choose q > 0 sufficiently large, so that ∂M * (y)/∂ν ≥ Ψ(z) for z ∈ ∂R l . Hence the following lemma holds.

Lemma 1. Suppose that a classical solution to system (6.5)-(6.7) of class C 2+α,1+α/2 exists on Ω × (0, T ). Let M (x, y, 0) ∈ [0, M * (y)), A(x, y, 0) ∈ [0, A * (y)). Then, for t ∈ [0, T ), 0 ≤ M (x, y, t) < M * (y), 0 ≤ A(x, y, t) < A * (y).

Taking into account the non-negativity of solution (proven above) the proof can be carried out via the maximum principle. ✷ 6.2. Existence of solutions Given U = ( U 1 , U 2 ), let P ( U ) denote the solution of the system

(L 1 U 1 , L 2 U 2 ) = Φ( U )
with the above boundary conditions and initial conditions U (x, y, 0) = U 0 (x, y) = (M 0 (x, y), A 0 (x, y)).

Let us assume that 0 ≤ U 0 (x, y) ≤ (M * (y), A * (y)). The local in time existence of such solution is guaranteed by Theorem IV.5.3 in [START_REF] Ladyzhenskaya | Linear and quasi-linear equations of parabolic type[END_REF]. In the set Ω T = Ω × (0, T ) let us consider the mapping U = P ( U ). (6.9)

Let M T = C 1+α,(1+α)/2 (Ω T ) and B = {U ∈ M : U -U 0 M ≤ 1}. From the Schauder estimates (see Theorem IV.5.3 in [START_REF] Ladyzhenskaya | Linear and quasi-linear equations of parabolic type[END_REF]) it follows that for T sufficiently small, T ≤ T * with some T * , the mapping (6.9) acts from B into B and it is a contraction (see, e.g. [START_REF] Kazmierczak | Existence of global solutions to a model of chondrogenesis[END_REF]). Hence it has a unique fixed point U in B. The function U is in fact of the class C 2+α,1+α/2 (Ω T ) and it is a solution of system (1.3)-(1.6). Obviously 0 ≤ U (x, y) ≤ (M * (y), A * (y)). (6.10)

Knowing only L ∞ norm of the solution, we can obtain an a priori estimate of the solution in the C 1+α,(1+α)/2 (R l ) norm. First, from Theorem 6.49 of section VI in [START_REF] Lieberman | Second Order Parabolic Differential Equations[END_REF] we conclude that the following estimate for A holds (6.11) for some constant W , where

A C 1+β,(1+β)/2 x,t (R l ×(0,T )) ≤ W F L ∞ (R l ×(0,T )) + A 0 C 1+β x (R l )
F = f (A)M -γA + b.
According to (6.10), f ∈ L ∞ (R l ×(0, T )). This estimate has a local character. Thus W can depend on T , but does not depend on l. Having this estimate and using Theorem IV.5.3 in [START_REF] Ladyzhenskaya | Linear and quasi-linear equations of parabolic type[END_REF], we can also estimate the C 2+α,1+α/2 (R l × (0, T )) norm of the function M : M C 2+α,1+α/2 (R l ×(0,T )) ≤ c 3 (T ).

(6.12)

Finally, we can estimate the C 2+α,1+α/2 (R l × (0, T )) norm of A:

A C 2+α,1+α/2 (R l ×(0,T )) ≤ c 4 (T ). (6.13)

Global existence of solutions

According to a priori estimates (6.12) and (6.13), the vector function U (x, y, T * ) has its C 2+α (R l )norm bounded by a finite constant. Using U (x, y, T * ) as a new initial condition and repeating the procedure we obtain the solution of the considered system in Ω T * +T 0 with some T 0 > 0. Continuing in this way, we obtain a global in time solution in R l × (0, T ) for any T > 0.

As we mentioned above, a priori estimates necessary for the global existence of solutions do not depend on l. Hence, passing to the limit, we obtain a global in time solution to the problem (1.3)- (1.6). The proof of uniqueness of solution is standard and it is left to the reader.

Proposition 3 . 3 .

 33 If the principal eigenvalue of problem (3.4)-(3.6) crosses the origin from negative to positive values, then the stationary solution M = 0, A = A 0 of problem (3.1)-(3.

S 1 .Proposition 3 . 4 .

 134 .15) Let µ(x, y) and α(x, y) be some functions continuous together with their second derivatives and (µ, α) = sup(x,y)∈Ω d M ∆µβµ ∂µ ∂x , S 2 (µ, α) = sup (x,y)∈Ω d A ∆α + f (α)µγα + b ∂α ∂xLet functions µ(x, y), α(x, y) satisfy conditions (3.16), (3.17). If c > max(S 1 (µ, α), S 2 (µ, α)) ,(3.18) then there exists a solution of problem (3.12)-(3.15).

  22) in the domain Ω N , with the boundary conditions y x = -N : M = M N (y), A = A N (y) (3.24) and the initial conditions M (x, y, 0) = M N (y), A(x, y, 0) = A N (y). (3.25)

Theorem 3 . 6 .

 36 Problem (3.12)-(3.15) has a solution if and only if c satisfies the inequality

Figure 1 :

 1 Figure 1: Left: beginning of the wave development. Right: the wave propagation.

Figure 2 :

 2 Figure 2: Left: comparison between the 2D-model with a small thickness (upper part) and the 1D-model (below). Right: dependence of the wave speed on the thickness of the strip.

Figure 3 :

 3 Figure 3: Bounded domain and boundary conditions (see the explanation in the text).

6. 2 . 1 .

 21 Local existence of solutionsThe local existence of solutions follows from the application of the contraction mapping principle.Let L 1 = ∂/∂td M ∆, L 2 = ∂/∂td A ∆, U = (U 1 , U 2 ) = (M, A), Φ(U ) = (-βU 1 , f (U 2 )U 1 -γU 2 + b).

  .28) Proof. To prove the lemma, we will write the problem for the new unknown functions

	u(x, y, t) =	∂M N (x, y, t) ∂x	, v(x, y, t) =	∂A N (x, y, t) ∂x
	and will show that its solution is positive. Differentiating problem (3.21)-(3.24) with respect to x,
	we obtain	∂u ∂t	= d M ∆u -c	∂u ∂x	-βu,	(3.30)

Acknowledgments

The authors are grateful to John McGregor for discussions about the inflammatory aspect of atherosclerosis and clarifying some important details about the risk factors.

(3.29)

The existence of such functions follows from the instability of the solution M 0 , A 0 . Indeed, let (µ 0 (y), α 0 (y)) be the eigenfunction corresponding to the principal (positive) eigenvalue of problem (3.4)- (3.6). Then the functions

satisfy conditions (3.26)-(3.28) for τ N sufficiently small. We should note, however, that instead of the condition

. This is a minor technical detail which can be arranged by replacing the function g(A) in (3.29) by g ′ (A 0 )A for A 0 ≤ A ≤ A 0 + ǫ, and then considering a limiting procedure as ǫ → 0.

By virtue of conditions (3.26)-(3.29) and of Proposition 2.4 adapted for the problem under consideration, the solution of problem (3.21)-(3.25) increases in time for each (x, y) ∈ Ω N . On the other hand, from inequalities (3.19), (3.20) it follows that it is estimated from above:

Therefore, it converges to a stationary solution (u N , v N ) of problem (3.21)-(3.25). From Lemma 3.5 below it follows that the functions u N (x, y) and v N (x, y) are non-decreasing with respect to x. Therefore there exists their limits as x → +∞. Since the limiting functions satisfy the problem in the section of the cylinder, and there are no other solutions that satisfy inequality (3.11) except for (M 0 , A 0 ) and (M s , A s ), then

We consider the sequence of solutions (u N , v N ) as N → -∞ and choose a convergence subsequence in order to obtain a solution on the whole axis. For this we introduce the shifted functions

where k N is chosen in such a way that

Such values exists due the boundary conditions at x = -N and the limiting values of the solutions at +∞. These new functions are defined for -N -k N ≤ x < +∞. Since u N (x, y) ≤ µ(x, y), x ≥ -N , then -Nk N → -∞ as N → ∞. Thus, we can choose a subsequence of the sequence (u N , v N ), for which we keep the same notations, which converges locally to some limiting functions (u 0 , v 0 ). They are defined in the whole cylinder Ω and satisfy problem (3.12)-(3.14). Moreover, they are non-decreasing with respect to x and u 0 (0, h/2) = 1 2 (M 0 + M s (h/2)). Hence the solution has limits (3.15) for x = ±∞. The proposition is proved.


Alfredo Rodríguez

Iñigo Gómara

Gianni Bellocchi

Raphaël Martin

Adela Martínez-Fernández

Alfonso Caballal

Jordi Doltra

Del Agustín

Margarita Prado

Ruiz-Ramos

Gianni Bellocchi

Raphaël Martin

Giani Bellocchi

Raphael Martin

Agustín Del Prado

Margarita Ruiz-Ramos

The challenges of modelling mixed management grasslands in North Spain under climate change

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Permanent grasslands are a very relevant cropping system in the North of Spain and support the main dairy farms in the country. Adaptation to climate change will be required given the projected changes of regional precipitation. To support such adaptation, modelling of these systems to generate high quality projections of the system performance is required. In the region to be simulated, grasslands are managed with a mixture of cuts and grazing. Several issues hinder the modelling of this type of systems: 1) the available data of grazing intensity presents large uncertainties; 2) there are few grassland models that allows flexibility to define a variable combination of cuts and grazing; 3) soil heterogeneity. This exercise has the goal of exploring the potential of PASIM to simulate mixed management grasslands. The model has been successfully applied in a permanent grassland in the French Massif Central, for a management that only contemplated cuts without grazing (Gómara et al., 2020). The model was calibrated using data from Villaviciosa (Asturias, Spain, 5º 26' 27" W, 43°28' 50" N, 10 m a.s.l.), located at northern Spain with a temperate climate. This calibration was used to simulate several grasslands locations distributed along the Cantabrian Sea. The soil information was obtained from [START_REF] Trueba | CIEMAT[END_REF]. The model was configured for the optimum management for mowing and nitrogen fertilization. The 1976-2005 period and the 2030-2059 period were selected. For the future period two representative concentration pathways emission scenarios (RCP, van Vuuren et al., 2011) were selected (i.e. RCP4.5 and RCP8.5). An ensemble of climate models will be used from the Coordinated Regional Climate Downscaling Experiment (CORDEX, Giorgi and Gutowski, 2015) previously bias-adjusting them by using the European observational database EOBS (Haylock et al., 2008) with the empirical quantile mapping method included at the climate4R R package (Iturbide et al., 2019). Modelling was challenging due to a combination of complexity (many processes involved) and uncertainty (observed data are difficult to generate). The results of the simulation exercise allow for assessing PASIM skill to reproduce the performance of these complex systems, as well as to

1

  Universidad de Castilla-La Mancha, Department of Economic Analysis and Finances, Toledo, Spain (alfredo.rodriguez@uclm.es) 2 Departamento de Física de la Tierra y Astrofísica, Universidad Complutense de Madrid 3 Grassland Ecosystem Research, French National Institute for Agricultural Research, INRAe 4 Servicio Regional de Investigación y Desarrollo Agroalimentario (SERIDA)

Powered by TCPDF (www.tcpdf.org)


Michel Morvan

Sylvain Sené
email: sylvain.sene@ens-lyon.fr

A Distributed Trust Diffusion Protocol for Ad Hoc Networks

Keywords: trust diffusion, ad hoc network, interaction graph

In this paper, we propose and evaluate a distributed protocol to manage trust diffusion in ad hoc networks. In this protocol, each node i maintains a "trust value" about an other node j which is computed both as a result of the exchanges with node j itself and as a function of the opinion that other nodes have about j. These two aspects are respectively weighted by a trust index that measures the trust quality the node has in its own experiences and by a trust index representing the trust the node has in the opinions of the other nodes. Simulations have been realized to validate the robustness of this protocol against three kinds of attacks: simple coalitions, Trojan attacks and detonator attacks.

Introduction and preliminaries

Ad hoc networks can be defined as collections of mobile nodes, distributed and independent, being able to communicate by radio transmission and to self-organise. They constitute networks with unstable infrastructure. Since the neighbourhood of every node changes over time, it is important to develop protocols that will help each node to identify reliably other nodes with which it can interact safely.

Let us present two definitions of the concept of trust given by the Webster's classic 1913: "Firm reliance on the integrity, ability, or character of a person or thing" and "Certainty based on past experience". Trust is therefore a natural notion acquired and systematically used by anyone to decide if an exchange with somebody else is conceivable or not. Let us remark that the trust we can have in somebody usually depends on the personal knowledges we have about the person and also on his/her reputation according to other people we have already met. This is usually the way that trust diffuses in human groups. For example, a researcher can trust the scientific opinion of a colleague while not trusting the opinion of the same colleague about the work of other scientists. The idea of this paper is to transfer this human trust diffusion protocol to ad hoc networks. Literature about trust management in P2P systems [START_REF] Buchegger | A Robust Reputation System for P2P and Mobile Ad-hoc Networks. P2P and Mobile Ad-hoc Networks[END_REF][START_REF] Sepandar | The EigenTrust Algorithm for Reputation Management in P2P networks[END_REF] and Web Services [START_REF] Resnick | Trust among Strangers in Internet Transactions: Empirical Analysis of Ebay's Reputation System[END_REF] presents some solutions using both personal and external opinions. However, these protocols do not offer protection against attacks made by a set of nodes working in coalition, for example diffusing wrong positive opinions about a dishonest node (called Trojan attack further in the paper). That is why we propose to study in this paper a management policy that also measures the trust quality by introducing a new level of trust that adds weights both on a node own experiences and on external knowledges this node has received from the others.

Let us call efficient a protocol that leads to prevent bad interaction and reliable a protocol that aims at favouring good interactions. This work consists in creating, studying and validating by simulations a new efficient distributed trust diffusion protocol for ad hoc networks using the double level of trust described above.

Section 2 presents and explains the trust diffusion protocol. Some results obtained on this protocol are exposed in Section 3. In particular, this section shows how the protocol can resist to three different types of attacks. Section 4 gives perspectives of this work.

In this section, we first present the different attributes used by the protocol in order to have the requested features to adapt on ad hoc networks. Then, we expose the protocol dynamics.

Protocol attributes

In order to respect the constraints induced by ad hoc networks, all nodes have to store their own knowledges about the others. This subsection introduces the different attributes that compose the notion of knowledge used in the proposed protocol and that a node i has about a node j at time t.

• Trust mark, denoted by tm i,j (t) ∈ [0, 1].

Node i stores a trust mark valued in [0, 1] about node j with which it has already interacted. This trust mark is used to know if an interaction has a good chance to be benefic. It is updated as a function of two major information: the own experience of node i and the external knowledges node i has obtained from other nodes during past interactions. It is initialized to

1 2 . • Trust index in trust mark, denoted by itm i,j (t) ∈ [0, 1].
Node i stores a trust index about its trust mark about node j. This trust index gives an indication about the reliability of the trust mark it has about node j. It is initialized to 1 2 . • List of external trust knowledges, denoted by et i,j (t) with et i,j (t)[k] ∈ [0, 1] and k ∈ [0, n].

Node i stores about node j a list of external trust marks. Such a list contains the trust marks owned by j about other nodes. This list has been transmitted to node i during the last interaction between i and j (during this interaction, i has symmetrically transmitted its trust marks to j). This list is initialized to ∅. Moreover, when nodes i and j symmetrically transmit to each other their trust marks, they also transmit their trust indices in these trust marks. These lists of trust indices are called lists of external trust indices in trust knowledges. Such a list owned by node i about node j is denoted at time t by iel i,j (t) with iel i,j (t)[k] ∈ [0, 1] and k ∈ [0, n] and is initialized to ∅.

• Trust index in list of external trust knowledges, denoted by iet i,j (t) ∈ [0, 1].

Node i stores a trust index about each list of external knowledges it has received from other nodes. Such trust indices give an indication about the reliability of these lists and are initialized to 1 2 .

• List of past interactions marks, denoted by pil i,j (t) with pil i,j (t)[k] ∈ {0, 1} and k ∈ [0, θ].

Node i stores about node j a list of marks encoding the results of its θ (where θ is bounded by θ max to avoid side effects) last interactions with the node. If the interaction was a success (the transmitted data was the expected data), the mark is 1, and it is 0 otherwise. The marks included in this list correspond to the own experience of node i about node j and are used to compute trust mark for future interactions. It is initialized to ∅.

We will call positive (resp. negative) an attribute (valued in [0, 1]) that is greater or equal to (resp. less than) 1 2 .

Protocol dynamics

We describe in this subsection the updating methods of the protocol attributes described above. First, let us note that the list of external trust knowledges (et) and the list of external trust indices in trust knowledges (iel) are simply transmitted through the network by nodes at each interaction. Then, the update of the list of past interactions is done as follows: when node i asks node j for data transmission, if the list of past interactions of i is not full, the mark affected to this new interaction is placed at the end of the list; otherwise, the first mark in the list is removed and the new mark is placed at the end. Now, let us focus on other attributes that need more complicated updating functions.

• Trust marks

When node i wishes to interact with node j at time t, it has not only to rely on its own knowledges but also on the external knowledges about node j it has already received from the other nodes. Thus, to update its trust mark about node j, node i needs to evaluate two different data. First, we define the personal knowledges PK as the simple average of the results of past interactions node i had with node j. So:

PK t+1 i,j = θ k=1 pil i,j (t)[k] θ
Now, let us also define the external knowledges. Let S * be the set of nodes k such as

{k = i, j, iet i,k (t) > 1 2 , et i,k (t)[j] < 1 2 , iel i,k (t)[j] > 1 2 }. At time t + 1,
the external knowledges EK that node i has about node j are computed as:

EK t+1 i,j = 1 |S * | n k=1 S * et i,k (t)[j].
The choice made here is to average only the bad opinions other nodes have about node j. Of course, we will restrict to nodes in which we trust the opinion and who themselves trust their own opinion about node j. Note that this choice leads the protocol to be efficient rather than reliable. Now, we update the trust mark by simply taking the weighted sum of these two knowledge values:

tm i,j (t + 1) = ε.[θ × PK t+1 i,j ]+ (E) ζ.[θ max × EK t+1 i,j ] ε.θ+ζ.θ max iff (E) =0
where ε and ζ are two parameters managing the weight of personal and external knowledges. So, the more ζ has a high value compared to ε, the more the diffusion will be fast. However, ε and ζ have to be well balanced to prevent the fact a node is easily disturbed if it has a majority of dishonest neighbours. Let us note that, if these two parameters are equal, external knowledges weight is greater than personal knowledges weight when the θ max first interactions have not already been executed, because θ max is the upper bound of θ. It seems to be reasonable to think that, when a node has few personal knowledges about an other with which it wants to interact, external knowledges are more important in its decision.

This trust marks evolution method shows that this protocol can not be less efficient than another one which does not implement trust diffusion. Indeed, if an interaction between two nodes i and j is possible, when i updates its trust mark about j, it can not obtain a higher new mark than one obtained only with its past interactions with j since it only uses negative external knowledges.

• Trust indices in personal trust marks

The updating method of the trust mark index (when i asks j for an interaction at time t+1) is the following. After having computed E(pil i,j (t)) and σ(pil i,j (t)) where E and σ refer respectively to the arithmetic average and the standard deviation, the new value of the trust mark index is obtained by a threshold algorithm. So, we fix a threshold named ω itm and we verify if the new trust mark tm i,j (t + 1) is in

]E(pil i,j (t)) -ω itm ; E(pil i,j (t)) + ω itm [
The index itm grows if α is in this range of values and decreases otherwise. In each case, updates are proportional to σ(pil i,j (t)); however, in order to go one step further in the "efficiency" direction, we have chosen the decreasing function to be faster than the increasing one.

• Trust indices in external trust marks

Node i will use the newly calculated trust mark on node j to adjust the trust index it has about the external trust knowledges stored for each other node. In other words, for each other node k, node i compares its new trust mark about j to the opinion k has about j. The more different (resp. close) it is, the more node i decreases (resp. increases) its trust index on the list of external knowledges given by k. A threshold algorithm (where the threshold is named ω iet ) is used here too. When i interacts with j, we consider the belonging of each et i,k (t)[j] to the following range of values:

]tm i,j (t + 1) -ω iet ; tm i,j (t + 1) + ω iet [ If et i,k (t)[j]
is in this range of values, the index iet i,k (t) grows, otherwise it decreases. There also, the decreasing goes faster than the growth.

Decision to interact

At time t + 1, if node i wishes to interact with node j, it verifies its trust mark about node j and decides to realize the interaction only if this trust mark is positive.

Simulations

In this section, we present the simulations that we have done on the trust diffusion protocol presented above (called Trudi) and compare the results of the latter to the so-called β-protocol in which nodes do not transmit their knowledges about other nodes. We are going to evaluate the efficiency of the protocols by measuring three parameters: (i) the number of honest nodes succeeding in identifying the dishonest nodes; (ii) the total number of bad interactions (transmitted data are not the expected ones) that the dishonest nodes succeeded to do; (iii) in some cases, we also measure the number of interactions needed for the honest nodes to identify the dishonest ones.

Let us precise the meaning of honest and dishonest we will use in the following. We will call action-dishonest (resp. action-honest) a node that wishes to give to every other unexpected data (resp. expected data) and opinion-dishonest (resp. opinion-honest) a node that wishes to give wrong opinions (resp. correct opinions) about some nodes to many others. A node is called dishonest if it is either action-dishonest or opinion-dishonest and totally-dishonest if it is both action-dishonest and opinion-dishonest.

Moreover, a trust protocol is called coherent if, when a little proportion of the system actors is action-dishonest, the honest ones succeed in the detection of those latter.

We will start by describing the simulation protocol used in the study before exposing some of the obtained results.

Simulation protocol

In order to run simulations not to far from real interactions in ad hoc networks, we are going to make the hypothesis that if we observe the set of interactions in the network after a sufficient amount of time, we will observe a network close either to a social network or to a complete network (where each node has interacted at least once with any other one). A social network is the personal or professional set of relationships between individuals. It has been shown [START_REF] Albert | Statistical Mechanics of Complex Networks[END_REF][START_REF] Mark | The Structure and Function of Complex Networks[END_REF] that these networks very often present the following characteristics: their degree distribution follows a power law (the probability that a node has a degree k is a decreasing power of k); they present a clustering coefficient (the probability for two nodes having a common neighbour to be connected) significantly higher than for classical random networks [START_REF] Strogatz | Collective Dynamics of 'Small-World' Networks[END_REF]. Of course, there is no evidence that these social networks precisely corresponds to the reality, but it makes sense to suppose that interactions in ad hoc networks will share some features with interactions in social networks.

In order to run simulations, that will lead at the end to such interaction networks, we determine the interaction graph a priori and we then realize the interactions on this graph. Thanks to this method, we control the final shape of the interactions superposition. Because of the above hypothesis, the chosen graphs are either power law degree distribution graphs with high clustering coefficient (such a network will be called "social network" in the following) or complete graphs.

Simulations results

The results presented in this subsection come from simulations of 1000000 interactions on 100 nodes networks.

In the first simulation, there is only one action-dishonest node which sends bad data with probability 1 5 and no opinion-dishonest node. In the social network, this node is given the highest degree. We measure the average and maximum number of interactions realized by the other nodes with this action-dishonest one. In Table 1, which presents the results of this simulation, we observe that Trudi is significantly better than the β-protocol. Indeed, the perception of dishonesty is about 45 times faster with Trudi than with the β-protocol. This result can be explained by the nature of trust diffusion in Trudi. Only negative marks and opinions are used to realize the updates of the protocol attributes. Therefore, as the weight of external knowledges is always more important than this of personal knowledges, honest nodes perceive quasi-directly the dishonesty.

Another simulation has been run by placing the dishonest node on the lowest degree node. In this case, there is no real difference between the two protocols because of the too small connectivity of the node. Thus, the interest of the protocol depends on the probability dishonest nodes have to give bad interactions and on the place of dishonesty in the network, an effect that we will call topology dependence. Let us now assume that some nodes make a coalition in order to attack the network, i.e. to use their knowledge about the protocol to break its efficiency or reliability. We are going to present in the following some results describing the effects that such an attack can have, with and without the Trudi protocol.

• Simple coalition attacks

A simple coalition attack manages a group of nodes in which the members are honest with each other but totally-dishonest with nodes that are not a part of the group, i.e. they protect themselves by giving positive opinions about the group members outside the group.

We have run simulations of this type of attack on complete networks. In each simulation, we have increased the size of the group of dishonest nodes (which send bad data with probability 1 5 ) to focus on the minimal size this group must have to destabilise the network, i.e. at least one honest node does not succeed in perceiving all the dishonest ones.

√ × 50 √ × 60 √ × 70 × × 80 × ×
Table 2 shows that Trudi is extremely robust against simple coalitions contrary to the β-protocol. With Trudi, the system can actually be destabilised only if the rate of good interactions is much less than the rate of bad interactions provided by dishonest nodes. Let us illustrate these results on a simple example which models a complete network where every dishonest node has a probability of 1 5 to provide a bad interaction. At each interaction, an edge has a probability of 1 9900 to be chosen (the interactions from i to j and from j to i are differentiated). So, when only thirty nodes are honest, the probability for two honest nodes to interact is 30×29 9900 ≈ 0.0878. When forty nodes are honest, this probability is 40×39 9900 ≈ 0.1575. So, in the second case, the system is not destabilised because the rate of good knowledges diffusion is close to the rate of bad interactions provided by dishonest nodes. This explanation has been verified by other simulations.

Moreover, note that the results presented in Table 2 are also valid for social networks despite the topology-dependence.

Let us now considerer more sophisticated attacks.

• Trojan attacks

A Trojan attack (by analogy to the mythical war of Troy) manages a group of nodes where one of them (Ulysses) is honest inside the group and totally-dishonest outside the group. Ulysses is protected by the other members of the group, i.e. they are action-honest but they always provide positive opinion to the honest nodes about Ulysses.

For this attack, two simulations have been executed (illustrating protocols robustness both on complete and social networks) to see the number of bad interactions (nbi) realized by Ulysses (which sends bad data with probability 1 5 ) and the maximum number of interactions (mni) executed by a honest node with Ulysses before perceiving its dishonesty. In the second simulation, Ulysses is the highest degree node and is protected by 9 high degree nodes. Table 3 validates the robustness of Trudi against Trojan attacks on complete networks as well as on social networks. However, other simulations where the dishonest group takes place on nodes of lowest degrees have shown that Trudi efficiency is not really significantly more interesting than one of the β-protocol. These results confirm that topology dependence is the main factor of the interest of the protocol in social networks.

• Detonator attacks A detonator attack is a Trojan attack where the dishonesty of Ulysses is started only after a determined time which is called the detonation time.

Our objective is to determine how many bad interactions (nbi) are provided by Ulysses, after its detonation, to be recognised as a mistrustable node and how many honest nodes (nhn) succeed in perceiving the dishonesty of Ulysses. Two simulations have been run to show the robustness of the two protocols against detonator attacks in complete networks where Ulysses sends bad data with probability 1 5 and is protected by 9 other nodes. The first (resp. the second) simulation illustrates this type of attack where Ulysses becomes dishonest after the 500000th (resp. 100000th) interaction. Table 4 shows that a too little number of interactions executed after the detonation of Ulysses node prevents many honest ones to perceive its dishonesty. Indeed, as the interactions are chosen randomly in the n (n -1) possible, in a 100 nodes system, 9900 (≈ 10000) different interactions exist. Consequently, an edge between a given node and Ulysses is potentially chosen fifty times in the first simulation. It is not sufficient for the given node to detect the dishonesty of Ulysses because of two reasons. Firstly, the given node has already a positive opinion about Ulysses and, secondly, the latter has a too small probability to provide bad interactions (p = 1 5 ). The second line of Table 4 validates this hypothesis. If we focus on the number of bad interactions, we can conclude our protocol is interesting because it provides a mean of 11.68 of bad interactions to each honest node.

Then, another simulation with the same features than the second one has been run to show the robustness of the two protocols against these attacks in social networks. Ulysses is here the node with the highest degree (with 15 honest neighbours). The last line of Table 4 shows that, with the two protocols, all Ulysses neighbours succeed in perceiving its dishonnesty. However, with Trudi, Ulysses node gives less bad data than with the βprotocol. Moreover, if we focus on the total number of interactions realized by the 15 neighbours of Ulysses after the detonation time, we see that this number is approximately three times more important with the β-protocol than with Trudi.

Conclusion and perspectives

The results presented in this paper show that such a trust diffusion protocol could have interesting efficiency properties. Besides, it could apply not only in ad hoc networks but also in other networks. It could be of interest in Internet and P2P systems. For example, it could be useful to measure trust users give to web sites and could broadcast sites dishonesty. However, to obtain such protocols, improvements have to be done.

In particular, future research should focus on the reduction of the size of the transmitted messages because the protocols "lightness" is a primordial point to avoid the network surcharge. So, we should try to find a good balance between reducing the transmitted messages size and the interest of the diffusion principle. Different solutions could be imagined. For example, since the goal of the protocol is to ensure efficiency, we could only diffuse the information about bad behaving nodes. Another solution would be to diffuse only a proportion (function of the network size) of these negative knowledges. Such hypotheses should be tested via new simulations processes to be validated.

Table 3 :

 3 Comparison of Trudi and the β-protocol after Trojan attacks. Type of network and Trudi β-protocol Ulysses error probability Complete network -20% nbi = 52 ; mni = 8 nbi = 1233 ; mni = 130 Social network -20% nbi = 14 ; mni = 25 nbi = 344 ; mni ∼ 900

Table 1 :

 1 Comparison of Trudi and the β-protocol after a simple attack.

	Trudi	β-protocol
	Avg ; Max 3.38 ; 10 153.04 ; ∼ 420

Table 2 :

 2 Comparison of Trudi and the β-protocol after coalition attacks. Number of dishonest nodes in the group Trudi β-protocol 40

Table 4 :

 4 Comparison of Trudi and the β-protocol after detonator attacks.

	Number of interaction	Trudi	β-protocol
	after the detonation
	Complete network	nbi = 836 ; nhn = 19	nbi = 837 ; nhn = 12
	500000
	Complete network	nbi = 1051 ; nhn = 75 nbi = 1411 ; nhn = 27
	900000
	Social network	nbi = 419 ; nhn = 15	nbi = 657 ; nhn = 15
	900000

Trust diffusion protocolBecause of the constraints inherent to ad hoc networks, the proposed model has to be based on three fundamental concepts: information decentralisation, diffusion and management of trust histories.