ASYMPTOTIC BEHAVIOUR OF THE SOLUTIONS TO A VIRUS DYNAMICS MODEL WITH DIFFUSION

. Asymptotic behaviour of the solutions to a basic virus dynamics model is discussed. We consider the population of uninfected cells, infected cells, and virus particles. Diﬀusion eﬀect is incorporated there. First, the Lyapunov function eﬀective to the spatially homogeneous part (ODE model without diﬀusion) admits the L 1 boundedness of the orbit. Then the pre-compactness of this orbit in the space of continuous functions is derived by the semigroup estimates. Consequently, from the invariant principle, if the basic reproductive number R 0 is less than or equal to 1, each orbit converges to the disease free spatially homogeneous equilibrium, and if R 0 > 1, each orbit converges to the infected spatially homogeneous equilibrium, which means that the simple diﬀusion does not aﬀect the asymptotic behaviour of the solutions.


Introduction.
A considerable number of studies have been made on virus dynamics. A lot of research results are obtained by mathematical models that assume spatial homogeneity. To consider the interaction among many agents, it is natural to start with ordinary differential equations, which assume the spatial homogeneity. However, diffusion effect may cause a change of qualitative behaviour of solutions. Therefore, models incorporating diffusion effect also should be discussed. For examples, Prüss, Zacher, and Schnaubelt [9] and Wang, Yang, and Kuniya [14] incorporate the spatial movement of some agents into their models.
In this paper, we consider the asymptotic behaviour of the solutions of the virus dynamics model with diffusion 2010 Mathematics Subject Classification. Primary: 35B40, 35K57; Secondary: 37B25, 92B99. Key words and phrases. Asymptotic behaviour, Lyapunov functions, reaction-diffution equations, virus dynamics model.

TORU SASAKI AND TAKASHI SUZUKI
in a bounded domain Ω ⊂ R n , where u 1 (x, t), u 2 (x, t), and u 3 (x, t) stand for the population density of uninfected cells, that of infected cells, and that of viruses, respectively, at the time t and at the point x ∈ Ω. We assume that the boundary ∂Ω is smooth, and impose the Neumann boundary condition where ν is the outer unit normal vector on ∂Ω, and u = (u 1 , u 2 , u 3 ). In the initial condition we assume that u 0 = (u 1,0 , u 2,0 , u 3,0 ) ∈ C 2 (Ω; R 3 ), ∂u 0 /∂ν| ∂Ω = 0, and each component of u 0 is non-negative and not identically equal to zero. Equation (1) is obtained by incorporating the diffusion terms to the ordinary differential equation (ODE) model with which Nowak and Bangham [8] considers the persistent infection of HIV. Thus, uninfected cells are produced at a constant rate λ, have a death rate m, and become infected at a rate βu 3 . The death rates of infected cells and viruses are denoted, respectively, by a and b. Each infected cell releases r virus particles when it bursts. Equation (4) has two equilibria u = (u 1 , 0, 0) and u * = (u * 1 , u * 2 , u * 3 ), where Here, R 0 = rβλ/(bm) is the basic reproduction number [1]. The equilibrium u is the disease free state. The equilibrium u * is the interior equilibrium, whose elements are all positive, if and only if R 0 > 1. We put Φ(s) = s − log s − 1 (s > 0).
We consider the asymptotic behaviour of the solution to (1), (2), and (3), mainly using the method in Latos, Suzuki, and Yamada [7]. Note that the spatially homogeneous steady solutions to (1) are u(x, t) ≡ u and u(x, t) ≡ u * . We can use the Lyapunov functions (6) and (7) to show that, as t → ∞, each solution u(·, t) to (1) converges in C(Ω; R 3 ) to the spatially homogeneous equilibrium solution whereû = u , the disease free spatially homogeneous equilibrium, if R 0 1, and u = u * , the infected spatially homogeneous equilibrium, if R 0 > 1.
Here, the Lyapunov functions play two roles. First, they establish the L 1 bound of the orbit. Second, they provide the omega limit argument by the pre-compactness of the orbit in C(Ω; R 3 ). For this property to show, we use the semigroup estimates and the standard parabolic regularity.
This paper is composed of seven sections and two appendices. Sections 2, 3, and 4 are devoted to the L 1 , L p , and W 2,p a priori estimates, respectively, and Theorem 1.1 is proven in Section 5. Section 6 shows some results of numerical simulations. Section 7 is a discussion, and the positivity of the solution and the semigroup estimate used for the proof are stated in Appendixes A and B, respectively.
2. L 1 bound of the solution. Having assumed the smooth initial value, we have unique existence of the classical solution local-in-time, which is extended to globalin-time if we can derive a priori bound on its L ∞ norm. In this section we show the following proposition.

TORU SASAKI AND TAKASHI SUZUKI
and along the classical solution u = u(x, t) to (1), (2), and (3), defined on Ω × [0, T ] for some T > 0. This use is justified by the positivity of the solution assured by the following proposition.
We note that the spatially homogeneous steady solutions to Equation (1) coincide with the equilibria of the corresponding ODE model (4). First we consider the case R 0 1.
Proof. For U(u) defined by (12), the derivativeU (1) (u) along the solution to (1) satisfiesU Hence it holds thatU (1) by Green's formula and the Neumann condition (2). Here,U (4) stands for the right-hand side of (9), and hence the last term of the right-hand side on (14) is non-positive. We thus end up witḣ Hence there exists a constant C 3 that satisfies This proves the proposition because Φ(s) ≥ 0 and u 2 , u 3 > 0.
In the case R 0 > 1, we have a similar result.
Proposition 2.4. Suppose R 0 > 1, and let u = (u i (x, t)) be the solution of (1) and (2), and (3). Then there exists a positive constant C 4 such that Proof. For V(u) defined by (13), its derivative along the solution of (1) becomeṡ by Green's formula and the Neumann condition (2). Here,V (4) stands for the right-hand side of (10), and hence the last term on the right-hand side on (16) is non-positive. We thus end up witḣ Hence there exists a constant C 5 that satisfies Since Φ(s) ≥ 0, this proves the proposition.
To complete the proof of Proposition 2.1, we note the following proposition.
where |Ω| denotes the volume of Ω, and s 0 is the root of Φ(s/2) = log 2 with s > 2.
Given v = v(x) > 0, let Then (18) implies We are ready to give the following proof, regardless of the value R 0 .
. Then Propositions 2.3 and 2.4 prove the proposition.
3. L p bound of the solution. The L 1 bound of the solution is improved to L p bound of the solution for 1 p < ∞, including its derivatives by the method of Latos, Suzuki, and Yamada [7]. In this section, we show the following proposition.
The first observation is the following proposition.
we have a constant C 9 such that Proof. Letū 1 =ū 1 (x, t) be the solution to by u 1 , u 3 > 0, and hence it follows that v 1 (x, t) > 0 from the strong maximum principle. Thus we obtain Henceforth, ∆ is provided with the Neumann boundary condition. Let L = for p 1. Here, e tL u 1,0 is the solution to ∂u/∂u = d 1 ∆u − (m/2)u with the initial value u 1,0 . Therefore, we have e tL u 1,0 p u 1,0 p . Now we use the semigroup estimate (Proposition B.1 in Appendix), which is equivalent to (19), the last integral on the right-hand side of (22) is bounded in t > 0. Since λ − (m/2)u 1 q λ|Ω| 1/q + (m/2) u 1 q , the assumption, u 1 (·, t) q < C 8 for t > 0, and inequalities (20) and (21) prove the proposition.
Assume that for an integer l (0 l n−2), there exists a constant C 11 such that u 1 (·, t) n n−l C 11 (t > 0).
If l < n − 2, then there exists a constant C 12 such that If l = n − 2, then for any p (1 p < ∞), there exists a constant C 13 such that
Then it holds that and hence, by Proposition 3.2, there exists a constant C 12 such that u 1 (·, t) p C 12 for t > 0. If l = n − 2 we take q = n/(n − (n − 2)) = n/2, which means 1/q − 2/n = 0. Then the conclusion of Proposition 3.2 holds for any p 1, and we obtain the result.
Proof. By Proposition 2.1, the assumption in Corollary 3.3 holds for l = 0. Thus repeating the use of Corollary 3.3 proves the corollary. Now we turn to the estimates on u 2 and u 3 . t)) be the solution to (1), (2), and (3), and 1 q < ∞. If there holds for a constant C 15 , then for p q with (19), i.e., we have a constant C 16 such that Similarly, if there holds for a constant C 17 , then for p q with (19) we have a constant C 18 such that Proof. The proof is similar to that of Proposition 3.2, because we have by (1), (2), and (3).
Corollary 3.6. Assume that for an integer l (0 l n−2), there exists a constant C 19 such that If l < n − 2, then there exists a constant C 20 such that If l = n − 2, then for any p (1 p < ∞), there exists a constant C 21 such that u 2 (·, t) p C 21 , u 3 (·, t) p C 21 (t > 0).
From the other semigroup estimate (Proposition B.2 and Remark B.3 in Appendix), on the other hand, it follows that ∇e r(d1∆−m) w p C 32 e −mr max{1, (d 1 r) −1/2 } w p , Here, by the Hölder inequality we obtain Then, by Proposition 3.1 (replaced p by 2p), (26) proves the proposition for i = 1. Using (25), we can similarly prove it for i = 2, 3.
Now we consider the derivatives of the second order. Here we may assume 1 < p < ∞. The first term of the right-hand side is estimated as ∆e t(d1∆−m) u 1,0 p = e t(d1∆−m) ∆u 1,0 p C 34 (t > 0), using u 1,0 ∈ D(∆). For the second term, we use Proposition B.2 to deduce Then, by Remark B.3 in Appendix, the Hölder inequality, and Proposition 3.1 and 4.1, we have for some C 36 and C 37 , This proves ∆u i (·, t) p C 38 (t > 0) for i = 1. Then inequality (27) for u 1 follows from the elliptic estimate. The other inequalities for u 2 and u 3 can be shown similarly.

5.
Proof of Theorem 1.1. First, we derive an L ∞ bound of the solution.
Proposition 5.1. For each solution u = (u i (x, t)) to (1), (2), and (3), there exists a constant C 39 such that Proof. By Morrey's theorem (Gilbarg and Trudinger [4], Corollary 7.11 and its notes), we have W 1,p (Ω) ⊂ C 1−n/p (Ω) with the continuous imbedding. By Propositions 3.1 and 4.1 for p > n, it holds that By Proposition 4.2, u i (·, t) (i = 1, 2, 3) is bounded in W 2,p (Ω) for any p 1 (See Remark B.3), and hence ∇u i (·, t) is bounded in W 1,p (Ω). Thus, by Morrey's theorem again, we obtain From the above a priori estimate, the solution u = (u i (x, t)) to (1), (2), (3) is global-in-time. Now we can use the invariance principle to study its asymptotic behaviour under the presence of the Lyapunov function (Henry [5,Theorem 4.3.4]), noting the following proposition. Proof. By Proposition 5.1, the family {u i (·, t)| t 0} is uniformly bounded and equicontinuous on Ω. Therefore, the Arzelá-Ascoli theorem implies the result. Proposition 5.2 implies that the omega limit set ω((u i,0 )) is nonempty, compact, invariant under the flow defined by (1)-(2), and is connected (Henry [5,Theorem 4.3.3]). Furthermore, the Lyapunov function is constant on this set. We are thus ready to complete the following proof.
6. Numerical simulations. In this section, we show some results of numerical simulations. The simulations were carried out with the Crank-Nikolson method, and confirmed with the explicit method (Smith [11,Chapter 2]). We consider the case n = 1 and Ω = {x| 0 < x < 10}. We set To depict the solutions tending to the spatially homogeneous equilibrium, we set the initial functions that fluctuate: where ψ 1 (x) = {2 + cos(4πx/10)}/3, which oscillates between 1/3 and 1. The graph of u i (x, t)'s, for this initial condition, are shown in Figure 1. For each u i 's, the spatial heterogeneity becomes smaller as time elapses, and u i (x, t) is almost spatially homogeneous at t = 8. Once the solution of (1) becomes almost spatially homogeneous, we can expect that its dynamics is governed by (4) and that the solution converges to u . We consider a situation where a small amount of viruses invade the middle of Ω. Hence we set where ψ 2 (x) = exp(−1/(1 − x 2 )) for −1 < x < 1, and ψ 2 (x) = 0 otherwise.
In Figure 2, each graph of u i (x, t)'s is shown from two view angles. The values of u i 's drastically change from t = 3 to t = 5, and during this time interval, the spatial heterogeneity also changes drastically. After the period, each u i becomes almost spatially homogeneous, and u i (·, t) converges to u * i as t → ∞ (i = 1, 2, 3). In Figure 2, the graphs of u i (x, t)'s seem flat for 0 t 2 and for 6 t 15, and it might be expected that u i (x, t)'s are spatially homogeneous during these periods.
However it may be due to the large ranges of values u i (x, t)'s. In Figure 3, we divide the the graph of u 3 (x, t) (Figure 2 (c)) into three time intervals, [0, 2.1], [2.1, 6], and [6,15]. Figure 3 (a) shows that u 3 (x, t) is not spatially homogeneous at least for 1.5 t 2. On the other hand, Figure 3 (c) shows that u 3 (x, t) is almost spatially homogeneous for 6 t 15. Then we can expect that the dynamics is governed by (4) and the solution converges to the spatially homogeneous equilibrium u * . 7. Concluding remarks. In Section 2 we confirmed that Lyapunov function used for the ODE model induces L 1 bound of the solution for the corresponding PDE system incorporating the diffusion. This structure is quite common to virus dynamics models. The proof of the compactness of orbits (from L 1 bound to W 2,p bound), on the other hand, depends on the reaction terms in our argument, which, however is applicable also to many other models. Then, the invariance principle of the omega limit set implies the convergence of the spatially inhomogeneous solution to the spatially homogeneous equilibrium. We shall discuss these points in a forthcoming paper, treating other virus dynamics models.
The function Φ(s) in (5) is different from the entropy density effective to the thermodynamical model, Φ(s) = s(log s − 1) + 1, which guarantees the Csizár-Kullback inequality valid to 0 < g, G ∈ L 1 (Ω) with g 1 = G 1 = M . It is now well-known that this inequality can be a basis to infer the exponential convergence of the spatially inhomogeneous solution to the spatially homogeneous equilibrium in L 1 norm (see [2]). Without (29), however, the convergence (11) is exponential, thanks to the theory of linearization. In fact, in the case of R 0 > 1, for example, the ODE theory has revealed that the stationary solution u * = (u * 1 , u * 2 , u * 3 ) is linearly asymptotically stable for the ODE part. This meas that the eigenvalues of the linearized matrix where σ(A) denotes the spectrum of A and δ > 0. If u = (u i (x, t)) is the solution to (1), (2), and (3), we obtain Here we recall that ∆ is provided with the Neumann boundary condition.
Appendix B. Semigroup estimates. We state estimates for the semigroup generated by the realization of the Laplacian in L p (Ω) (1 p < ∞) under the Neumann boundary condition. We assume that the spatial domain Ω is bounded in R n , and its boundary is smooth.
The following estimate follows from Lemma 3 in Part I in Rothe [10].