Existence, uniqueness, stability and asymptotic behavior of solutions for a mathematical model of atherosclerosis

We study an atherosclerosis model described by a reaction-diffusion system of 
three equations, in one dimension, with homogeneous Neumann boundary conditions. 
The method of upper and lower solutions and its associated monotone iteration (the monotone iterative method) are used to establish existence, 
uniqueness and boundedness of global solutions for the problem. 
Upper and lower solutions are derived for the corresponding steady-state problem. 
Moreover, solutions of Cauchy problems defined for time-dependent system are presented as alternatives upper and lower solutions. 
The stability of constant steady-state solutions and the asymptotic behavior of the time-dependent solutions are studied.


1.
Introduction. Mathematical modeling of the atherosclerosis processes leads to complex systems of flow, transport, chemical reactions, interactions of fluid and elastic structures, movement of cells, coagulation and growth processes and to additional complex dynamics of the vessel walls.
Partial differential equations have been used to model this complex process. As an example found in the literature, we can cite [9], where is presented a model consisting of reaction-diffusion equations describing the concentration of macrophages and cytokines in the intima (a vessel layer), which leads to an inflammatory disease. A model leading to the atherosclerotic plaque formation and the early atherosclerotic lesions was suggested by [2] and [3], respectively. Systems of convectionreaction-diffusion equations were used to describe the transport and the concentration of oxidized low density lipoproteins (LDL), macrophages, foam cells and signal in the intima. Recently, in 2014, a more complex and realistic model was presented in [8]. The authors used reaction-diffusion equations to describe the distribution of substances, such as LDL, high densities lipoproteins (HDL), oxidized LDL, free radicals among others, in the intima, and convection-reaction-diffusion equations for each species of cells, such as macrophages, T cells or foam cells.
In 2007, results of stability analysis of constant equilibrium solutions for a system of two reaction-diffusion equations in 1D (not considering the intimal thickness) were presented in [9]. In there, the existence of traveling waves was also proved. The existence of traveling waves for a 2D model was proved later in 2009 (see [10]), where the authors suggested a model consisting of a reaction-diffusion system in a strip with nonlinear boundary conditions. The existence and uniqueness of global solutions in Hölder spaces was proved in 2012, [11].
Our main contribution in this work is to prove the existence, uniqueness and boundedness of global solutions, as well as the analysis of the stability properties, and the long time behavior of the solutions for a system of three reaction-diffusion equations in 1D, provided with homogeneous Neumann boundary conditions. Such system describes the atherosclerosis inflammatory process by representing the evolution of the concentrations of oxidized LDL, macrophages and signal path between cells inside the arterial wall. This 1D model is a simplified version of the model presented in [2].
For reaction-diffusion equations of the type where the nonlinear function F does not depend on the gradient ∇U, a priori estimates of the uniform norm are sufficient to prove global existence, [19]. Ladyzenskaja [13], Friedman [5], Krylov [12], among others, studied the existence and uniqueness of classical solutions of reaction-diffusion equations for smooth initial data local in time, in Hölder spaces. Results for a more general setting in Banach spaces can be found in [6]. Based on the monotone sequences method, Sattinger [20] and Pao [15] proved the existence, uniqueness and boundedness of global solutions, whenever the function F (x, U ) is differentiable. They also characterized the stability and instability in terms of upper and lower solutions.
Basically, the monotone iterative method consists in using an upper or a lower solution as the initial iteration in a suitable iterative process to obtain a monotone sequence, which converges to a solution of the problem.
In the present work, we apply the general monotone sequence method as presented in [15]. To this end, we need to verify that our model fulfills several non-trivial technical conditions, before concluding our analytical results. This paper is organized as follows. In section 2, we recall some preliminaries and notations that will be useful in the sequel. Section 3 is devoted to the description of the mathematical model we are interested on. First, in subsection 3.1, we present a 1D mathematical model, which describes the atherosclerosis inflammatory process, together with some mathematical assumptions. Then, in subsection 3.2, the model is rewritten as a parabolic problem, in order to be consistent with the notation used later on. Section 4 comprises the description of the monotone iterative method in the frame of our problem (subsection 4.1), as well as the proof of existence, uniqueness and boundedness of global solutions, for a given pair of upper and lower solutions (subsection 4.2).
Section 5 addresses the construction of those upper and lower solutions. For this purpose, it is shown that: the solution of the steady-state problem exists and it is bounded by the corresponding upper and lower solutions (subsection 5.1); the timedependent solution is bounded by the upper and lower solutions of the steady-state problem (subsection 5.2). From the steady-state problem, positive upper and lower solutions are determined in subsection 5.3. Alternatively, in subsection 5.4, another pair of upper and lower solutions is constructed from the Cauchy problems. Finally, results about stability and asymptotic behavior of the solutions are presented in section 6.

Preliminaries and notations.
Let Ω be an open domain in R d . We denote by ∂Ω the boundary of Ω and by Ω its closure. For each T > 0, let Ω T = Ω × (0, T ] be a domain in R d+1 , ∂Ω T = ∂Ω × (0, T ] and Ω T the closure of Ω T . As usual, the set of all continuous functions whose partial derivatives up to the m th order are continuous in Ω is denoted by C m (Ω), which is a Banach space for the norm is the seminorm in C m (Ω) and α is the multi-index of order m. The space C 0 (Ω) ≡ C (Ω) denotes the set of continuous functions in Ω.
We denote by C m,l (Ω T ) the space of functions whose space and time derivatives up to the order m and l, respectively are continuous in Ω T . In particular, we will consider the Banach space C 2,1 (Ω T ) with the the norm defined by where U x and U xixj denote partial derivatives of first and second order, respectively. Similar notations are used for C m (Ω) and C m,l (Ω T ).
The space of continuous functions in Ω T , equipped with the norm is denoted by C(Ω T ). A function U ∈ C (Ω) is said to be Hölder continuous of order δ ∈ (0, 1) if the seminorm is finite. When δ = 1, we obtain the usual definition of Lipschitz continuity.
For 0 < δ < 1 and m = 0, 1, 2, ..., the Hölder space C m+δ (Ω) for which the norm is finite, is a Banach space. Here [U ] C m+δ (Ω) is defined as The set of all Hölder-continuous functions in Ω with finite norm is denoted by C δ (Ω). Let F : Ω T × R n → R. The function F satisfies the Hölder condition in Ω if there exists a constant H (the Hölder constant) such that and in Ω T , We denote by C k (Q) the product function space C k (Q) n , where n ∈ N, k is any nonnegative constant and Q can be Ω T or Ω T . For any vector function U = (U 1 , ..., U n ) in C k (Q) the components U 1 , ..., U n are all in C k (Q). The norm of U is defined by where |U | k is the norm in C k (Q). When k is equal to zero, m or δ, where δ ∈ (0, 1) and m is a positive integer, the corresponding function spaces are denoted by C(Q), C m (Q) and C δ (Q), respectively. These function spaces equipped with the above norm corresponding to k = 0, m and δ are all Banach spaces. When Q = Ω T we denote by C 2,1 (Ω T ) the product space of C 2,1 (Ω T ) taken n times. In this work we will use the special case n = 3 for the above notations of product spaces.
3. The mathematical model. In this section we will describe in detail the mathematical model for the atherosclerosis inflammatory process. Moreover, we will rewrite the problem as a parabolic system, according to the notation that will be used in the next sections.
3.1. A 1D mathematical model of atherosclerosis. Based on the model proposed by Calvez and co-authors in [2] and [3], we define a 1D simplified model, consisting of the following system of three reaction-diffusion equations along with the boundary conditions and initial conditions The functions C ox , M and S are, respectively, the concentrations (or densities) of oxidized LDL, macrophages and signal (a generic chemoattractant which gathers the cytokynes), which are continuously differentiable in t and twice continuously differentiable in x.
The second term on the left-hand side of each reaction-diffusion equation of the system (1), represents the diffusion term, and d ox , d M and d S are the diffusion coefficients, which are positive constants. The first term on the right-hand side is the reaction between ox-LDL and macrophages, where β is a positive constant of proportionality.
The function τ in (1a) is the permeability of the blood vessel, and C LDL is a given LDL-cholesterol concentration. Here we are assuming that τ depends on the wall shear tress (WSS), which is the mechanical force imposed on the endothelium by the flowing blood. As we know, low WSS favors the penetration of both LDL and macrophages [18]. The permeability function τ is defined to be a smooth and nonnegative function of W SS 1 .
In the equation (1b), we are assuming that the incoming monocytes immediately differentiate into macrophages and the recruitment of new monocytes depends on a general pro-inflammatory signal which gathers both chemokines and cytokines. This signal acts through the function g, which is defined, in order to impose a limit in the macrophages recruitment, as g (S) = S/ (1 + S) and g(S) > 0 for S > 0, g(0) = 0, and g(S) → 1 as S → ∞.
The second term in the equation (1c) denotes the natural death of the cells and λ is a degradation rate. The starting point of the signal emission is assumed to be a too high oxidized LDL concentration. This is described by the third term γ C ox − C th ox , where C th ox corresponds to a given ox-LDL quantity and γ is an activation rate. In order to have an inflammatory response C ox should be greater than C th ox . The physical meaning of the boundary condition in (1d) is that the boundary surface is completely isolated so that there is no flow across the boundary.
The functions C ox0 (x), M 0 (x) and S 0 (x) defined in (1e) are smooth and nonnegative (to ensure the nonnegativity of the solutions, for physical reasons, see [4]) satisfying the boundary condition (1d) at t = 0.
Thus, the system (1) with the boundary condition (1d) and initial condition (1e) can be written as for i = 1, 2, 3. Mathematically, the function τ is defined in such a way that F 1 (x, C 1 , C 2 ) is Hölder-continuous in Ω. The initial condition C i,0 (x) is nonnegative Höldercontinuous in Ω and satisfies the boundary condition at t = 0 and the concentration function C i (x, t) ∈ C(Ω T ) ∩ C 2,1 (Ω T ). The system is closed with the homogeneous Neumann boundary condition.
The time-dependent system (3) is often referred to a weakly coupled parabolic system (since each equation contains derivatives of just one component and the system is coupled only through the terms which are not differentiated) [17].
Since the initial conditions C i,0 (x) for i = 1, 2, 3 are nonnegative, and for all ox ≥ 0, so, the nonnegativity of the solutions of (3) is preserved with time ( [4] and [15]).
Using the split notation of the vector function C, we rewrite the function F i as where a i , b i are two nonnegative integers such that a i + b i = 2 and [C] ai , [C] bi denote, respectively, the a i and b i -components of the vector C.
The split form of C varies with respect to i and is determined by the quasimonotone property of F i .

4.
Existence, uniqueness and boundedness of solutions. We use the method of upper and lower solutions and its associated monotone iteration to establish an existence theorem. The analytical studies in this section are a result of the application of the general theory of parabolic problems presented in [20] and [15]. We will therefore need to verify a certain number of assumptions. Before doing so, we recall some definitions and analyze their meaning in our case.
4.1. Monotone iterative method. By using an upper or a lower solution as the initial iteration in a suitable iterative process the resulting sequence of iterations is monotone and converges to a solution of the problem. In fact, this corresponds to the monotone iterative method that can be used to establish an existence theorem for the initial-boundary-value problem (3). The definition of upper and lower solutions and the construction of monotone sequences depend on the quasimonotone property of the reaction function in the system. Therefore, we need to verify if the reaction function .., F n ) is said to possess a quasimonotone nondecreasing property (see [15]).
In our case, for i = 1, the partial derivative ∂F1 ∂C2 = −βC 1 ≤ 0 and since F 1 is independent of C 3 , we have Finally, for i = 3, ∂F3 ∂C1 = βC 2 + γ ≥ 0 and ∂F3 ∂C2 = βC 1 ≥ 0 and consequently, So, we conclude that the reaction function Now the definition of upper and lower solutions follows: The differential inequality (8a) is written explicitly as and we rewrite (8b) as well For a given pair of coupled upper and lower solutions C, C, the sector C, C (in the componentwise sense) is defined as the functional interval where the inequalities between vectors should be interpreted in the componentwise sense.
To ensure the existence of a solution for (3) it will be necessary to impose some condition on the reaction function, such as an Hölder or Lipschitz condition. Indeed, the reaction functions (2) are continuously differentiable in (R + ) 3 with respect to C, therefore they satisfy the one-sided Lipschitz condition where Ci ≤ C i ≤ Ci ≤ Ci, Ri = sup ∂F i ∂C i : Ci ≤ Ci ≤ Ci, for i = 1, 2, 3 and the Lipschitz condition, whith Lipschitz constant L i = R i , (because F i is C 1 -function) and consequently they are Hölder-continuous in C, C .
In the next paragraphs we describe the construction of monotone sequences. By choosing C (0) = C 1 , C 2 , C 3 and C (0) = C 1 , C 2 , C 3 as two initial iterations, the upper and lower sequences and for k = 1, 2, ..., are obtained through the iterative process with boundary and initial conditions given by i (x, 0) at t = 0 (17b) for k = 1, 2, ... and i = 1, 2, 3, where R i , the Lipschitz constants, are given by Since for each k the two systems (15) and (16) are uncoupled scalar linear problems (that heve a unique solution Ω T ) and by the properties of F i , the sequences and C (k) are well defined (for the proof see [15], page 58). The following positivity lemma for the linear scalar problem, which is a consequence of the maximum principle, is one of the bases for the monotone iterative method.
Lemma 4.2 (Positivity). Let W ∈ C(Ω T ) ∩ C 2,1 (Ω T ) and let R be a bounded function in Ω T such that in Ω T , Then W (x, t) ≥ 0 in Ω T . Moreover W (x, t) > 0 in Ω T , unless it is identically zero.
for every k. Moreover, for each k, C (k) and C (k) are coupled upper and lower solutions of (3). (15a) and (9a), From the boundary and initial conditions (8c), (8d), we have ∂ Following the same argument, from (15b) and (15c) we find that U 2 ≥ 0 and U 3 ≥ 0 in Ω T . This proves the relation Similarly, using the property of a lower solution and the sequences (16), we have since, by the one-sided Lipschitz condition (12), By a similar argument The function g C and g is monotone nondecreasing.
This shows that in Ω T .
By induction, we prove that the sequence C (k) is monotone nonincreasing and C (k) is monotone nondecreasing. The monotone property follows for every k = 0, 1, 2, . . .. By the monotone iteration process (15a)-(16c) for each fixed k, and by the onesided Lipschitz condition (12), we have

The monotone property (24) implies
By the same argument Consider now (17a) and (17b), we conclude that C exist and satisfy the relation Due to the quasimonotone property of F, the one-sided Lipschitz condition (12), the Lipschitz condition (13) and the monotone property (24) we can conclude, by applying the existence-comparison theorem proved in [15] (page 429), that the limits of C (k) and C (k) coincide and yield to a unique solution of (3). The uniqueness is ensured since the F is a C 1 -function for all C i ≥ 0. The solution is global, since T is arbitrary. These conclusions can be summarized in the following result.
Theorem 4.4. Let C = C 1 , C 2 , C 3 and C = C 1 , C 2 , C 3 be a pair of nonnegative coupled upper and lower solutions of (3) and let F = (F 1 , F 2 , F 3 ) be quasimonotone satisfying the global Lipschitz condition (13). Then the two sequences C (k) and C (k) given by (15)-(17b) converge monotonically to a unique solution C = (C 1 , C 2 , C 3 ) with

5.
Construction of upper and lower solutions. As we have seen above, the existence and uniqueness of the global solution to (3) is ensured whenever we can provide a pair of coupled upper-lower solutions of (3). In order to construct such pair of solutions, we will first turn our attention to the steady-state problem.

5.1.
Existence and boundedness of steady-state solutions. Since F is a quasimonotone function independent of t, a function C s = C(x) is said to be an equilibrium solution or a steady-state solution of (3) if it satisfies the corresponding steady-state problem, which is a system of ordinary differential equations.
Following the same approach used for the time-dependent problem (3) it is possible to construct two monotone sequences for the steady-state problem, using a similar iterative process. By neglecting the initial condition (17b) and dropping the time derivative terms in (15) and (16), for i = 1, 2, 3, a monotone argument shows that the corresponding steady-state sequences, C (k) s and C (k) s converge to some functions C s and C s in the same monotone way as for the time-dependent system.
The limits C s = C 1 , C 2 , C 3 and C s = (C 1 , C 2 , C 3 ) satisfy the boundary condition in (27) and the coupled equations However, neither C s nor C s is a solution of (27) unless C s = C s (in the componentwise sense). These limits C s , C s , are called quasisolutions of (27). The existence proof is based on the Schauder fixed point theorem, see [15] (page 440) for more details.
The relation is a result of the monotone property for the problem (27). Proof. Since the limits C s , C s are themselves upper and lower solutions of (27), there exists at least one solution C s in C s , C s . In fact, every solution of (27) in C s , C s is actually in C s , C s .

An invariant set of the time-dependent problem.
A sector V, W between two functions V, W with V ≤ W is called an invariant set of (3) if for any C 0 ∈ V, W the corresponding solution C (x, t) of (3) remains in V, W for all t > 0.

Remark 1.
If C 0 in an attractor (see section 6) then C 0 is also an invariant set, the reverse is not true.
The following theorem ensures that C s , C s is an invariant set of (3).
Theorem 5.3. Let C s and C s be a pair of nonnegative coupled upper and lower solutions of (27) and let C s , C s be its quasisolutions. Then the solution C(x, t) of (3) remains in C s , C s for all t > 0 whenever it holds at t = 0. Moreover it satisfies the relation Proof. If the initial function C 0 is in C s , C s , any coupled upper and lower so- in Ω, where R i is defined as in (18)- (20), with the boundary conditions V i (x) = 0 for i = 1, 2, 3. This problem is defined from the system (27), since The scalar boundary value problem (32a) has a unique positive solution V 1 . The reaction functions in the equations (32b) and (32c) are quasimonotone nondecreasing. Therefore, knowing V 1 , the system of these two equations has a unique solution (V 2 , V 3 ).
Consider the eigenvalue problem in one dimensional space with Neumann boundary conditions where µ i is the principal eigenvalue and φ i are the corresponding normalized eigenfunctions, (max φ (x) : x ∈ Ω = 1), for i = 1, 2, 3.
Proof. Let (V 1 , V 2 , V 3 ), be an upper solution of (27). Clearly it satisfies the upper inequality (28c) and for (28a)-(28b) we have which are always true, since for i, j = 1, 2, with i = j, is a lower solution of (27), then the lower inequality (29c) is directly satisfied. For the inequalities (29a)-(29b) we have which are always satisfied as well.
The boundary condition (29d) is directly satisfied by the definition of (32) and (33).

5.4.
Solutions of Cauchy problems as upper and lower solutions. The construction of upper and lower solutions when the boundary condition is of the homogeneous Neumann type can be done by considering the following systems of ordinary differential equations where By the mixed quasimonotone property of (F 1 , F 2 , F 3 ) and from the definition of upper and lower solutions, the spatially homogeneous upper and lower solutions W = (W 1 , W 2 , W 3 ) and Q = (Q 1 , Q 2 , Q 3 ) are required to satisfy the coupled relation where W ≥ Q.
Corollary 3. Let W(t) and Q(t) be the respective solutions of (34) and (35), with Proof. Clearly W(t) and Q(t) are coupled upper and lower solutions of (3). Therefore, by theorem 4.4 the conclusion of the theorem holds.
If F i is independent of x then problems (34) and (35) are reduced to the same form: 6. Stability and asymptotic behavior of solutions. For a given steady state solution C s , it is important to know what is the set of initial functions, for which the corresponding time-dependent solution converges to C s as t → ∞. This leads to the discussion of stability (Lyapunov stability) and asymptotic stability of a steady-state solution and its stability region. The steady-state solution C s = (C s 1 , C s 2 , C s 3 ) is said to be stable if for any constant > 0 there exists a δ > 0 such that for all t > 0, whenever then C s = (C s 1 , C s 2 , C s 3 ) is said to be asymptotically stable. The solution C s is unstable if it is not stable.
We recall that, for an asymptotically stable steady-state solution C s , the set of initial functions C 0 for which the corresponding solution C (x, t) of (3) converges to C s , is called the stability region (or domain of attraction) of the steady-state solution C s . When this domain is the whole space of all initial functions, C s is called a global attractor.
In this section we will analyze the stability of constant stationary solutions C s and the asymptotic behavior of the time-dependent solution C(x, t).
6.1. Stability of constant steady-state solutions. When F is independent of x, problem (3) becomes We note that, since we are dealing with homogeneous Neumann boundary conditions, studying the stability of a constant steady state solutions is equivalent to study the stability of equilibrium solutions to a Cauchy type problem In fact, if E is a stable equilibrium solution of (43), then for any > 0 there exists δ > 0 such that Let W 1 (t), W 2 (t) be the solutions of (43) with Then W 1 , W 2 are upper and lower solutions of (42), whenever Hence the solution C(x, t) satisfies the relation |C(x, t) − E| < when |C 0 − E| < δ.
6.1.1. Constant steady-state solutions. If we assume that the wall shear stress is uniform, then we are allowed to also assume that the permeability function τ is constant, in Ω. Therefore, F is independent of x. For simplicity of notation, let h = τ C LDL and m = C th ox . To find the equilibria, we need to solve the following system: If all parameters are positive, then system (44) has a unique positive solution E = (C 1 , C 2 , C 3 ), with , Let us consider now C LDL = m = 0. Then C 1 = C 3 = 0 and therefore, we have infinitely many equilibria represented by E * = (0, A, 0), where A can be any nonnegative value.
In this situation, since we have no concentration of oxidized LDL in the intima, which implies that no signal is sent to recruit the monocytes, we have no activation of the inflammatory process.
Since the starting point of the signal emission is assumed to be a too high oxidized LDL, C 1 > m, we have no recruitment of macrophages, and consequently no activation of the inflammatory process.
Therefore, we have a unique equilibrium solution (m, 0, 0), depending on the parameter m, with m > 0, which is equal to E * * . We conclude that, when we have no penetration of LDL in the intima, h = 0, the activation of the inflammatory process will not take place.
We are interested in studying the case when the system leads to activation of the inflammatory process, more precisely when all parameters are positive. This is the case of the positive equilibrium solution E, which exists and is unique, depending on the values of the parameters. Indeed, an equilibrium solution E l is stable if the real part of each eigenvalue of J F (E l ) is negative [1].
When the equilibria are E * and E * * , we have which implies that E * and E * * are unstable, since the eigenvalue µ = 0 belongs to the spectrum of J F (E * ) and J F (E * * ).
For the equilibrium solution E = (C 1 , C 2 , C 3 ), let the matrix B = J F (E).
To analyze the spectrum of we need to consider the characteristic equation where a 1 = βG 1 + βG 2 + λ, The coefficients a 1 and a 3 are always positive.
Assuming that λ ≥ (1 − h) 2 , the coefficient a 2 is strictly positive. Therefore, the eigenvalues of the matrix B are all real and negative (Descartes method) and consequently, the equilibrium solution E is stable.
6.2. Asymptotic behavior of the time-dependent solution. We want to prove that when the initial function in the problem (3) is an upper or a lower solution of the corresponding steady-state problem, the time-dependent solution is monotone in t and converges to the steady-state solution.
6.2.1. The uniqueness of steady-state solutions. The uniqueness of solutions to problem (27) is not guaranteed by theorem (5.2) even restricted to a sector C s , C s .
In fact, the quasisolutions, C s , C s , are not necessarily the same even if F is a C 1 -function.
To ensure uniqueness it is necessary to impose some additional condition on F. When the reaction function is quasimonotone nondecreasing, both C s , C s are true solutions of (27) and they are called maximal and minimal solutions in C s , C s . Since F is mixed quasimonotone, we can embed the system (3) into an extended system of 2 × 3 equations, such that, its reaction function has the quasimonotone nondecreasing property.
Using the splitting form for F i , (5), the extended problem for the steady-state system (27)  The components [C] ai and [C] bi , were computed in (7a)-(7c). Therefore, the extended problem (45) becomes in Ω T , with the boundary conditions (45c), for i = 1, 2, 3. Let is quasimonotone nondecreasing in The condition C s = C s is ensured if only if the extended problem has a unique solution. Monotone sequences can be obtained from the iteration process for the transformed system using the same argument as in (27) Applying the convergence results for quasimonotone nondecreasing systems we obtain the uniqueness of the extended problem. These results can be found in [16] for a general parabolic system. 6.2.2. An attractor of the time-dependent problem. The following theorem gives a sufficient condition for ensuring that, the sector C s , C s , between the two quasisolutions C s and C s , is an attractor of the time-dependent problem (3) for all C 0 in C s , C s . Theorem 6.1. Let C s , C s be a pair of coupled upper and lower solutions of (27), and let C s and C s be its quasisolutions. Then, for any C 0 ∈ C s , C s : (i) the corresponding solution C (x, t) of (27) satisfies the relation where x ∈ Ω; (ii) if C s = C s ≡ C s , then C s is the unique solution of (27) in C s , C s and the corresponding solution C (x, t) of (3) converges to C s (x) as t → ∞.