Heterogeneous Viral Environment in a HIV Spatial Model

We consider the basic model of virus dynamics in the modeling of Human Immunodeficiency Virus (HIV), in a 2D heterogenous environment. It consists of two ODEs for the non-infected and infected $CD_4^+$ $T$ lymphocytes, $T$ and $I$, and a parabolic PDE for the virus $V$. We define a new parameter $\lambda_0$ as an eigenvalue of some Sturm-Liouville problem, which takes the heterogenous reproductive ratio into account. For $\lambda_0<0$ the trivial non-infected solution is the only equilibrium. When $\lambda_0>0$, the former becomes unstable whereas there is only one positive infected equilibrium. Considering the model as a dynamical system, we prove the existence of a universal attractor. Finally, in the case of an alternating structure of viral sources, we define a homogenized limiting environment. The latter justifies the classical approach via ODE systems.

1. Introduction. The acute infection by the Human Immunodeficiency Virus (in short HIV) is characterized by a huge depletion of the CD + 4 T -lymphocytes (CD 4 ) and a peak of the virus load [7]. After few weeks, these two components reach a steady state which characterizes the asymptomatic phase of the infection. Before the availability of highly active antiretroviral therapy, this later phase lasted after 10 years in median with an accelerated decrease of CD 4 and an increase of virus load. A substantial number of nonlinear ODE systems have been suggested by Perelson et al. (see [21,22]) to model the complex dynamics of HIV-host interaction. For instance, such models have been used to estimate the infected cell half-life and the viral clearance during antiretroviral therapy [11,20,29], or to understand the dynamics during acute infection [23].  The common basic model of viral dynamics [2] includes three variables: T , the non-infected CD 4 , I, the infected CD 4 , and V , the free virus: (1.1a) This model, describing the interaction between the replicating virus of HIV and host cells, is based on some simple hypotheses. Non-infected target cells are produced by the thymus at a constant rate α and die at a rate µ T . By contact with the free virus particles (virions) they become infected at a rate proportional to their abundance, γV T (see Fig. 1). These infected cells die at a rate µ I and produce free viruses during their life-time at a rate N . Free particles are removed at a rate called the clearance µ V . All these parameters are generally positive constants. This simple model study has led to interesting results (see [18,28]) and suggested a treatment strategy (see [2]).

Figure 1. The HIV life cycle
It is easily seen that the system has two equilibria: (i) the non-infected steady state which corresponds to a non-negative equilibrium in case of no infection; (ii) the infected steady state also called seropositivity steady state, corresponding to a positive equilibrium in case of infection. Some authors (see e.g., [2,18]) have considered the basic reproductive ratio R 0 : a dimensionless parameter defined by epidemiologists as the average number of infected cells that derive from any one infected cell in the beginning of the infection [18, p. 16]. Stability properties of the two steady states are usually studied around this quantity: if R 0 < 1 the non-infected steady state is stable, if R 0 > 1 the infected steady state has a biological meaning and it is stable, while at R 0 = 1 both steady states coincide. So, R 0 = 1 is a bifurcation point (see Fig. 2). Further   models have been used involving other populations present in the immune system (see [17,18,4]). However, these models assume that the populations T, I, V are homogeneous over the space for all time, which is a common, but not a very realistic, assumption. Actually, the interaction between the virus and the immune system (either as a target with CD 4 or as an agent for controlling infection) is localized according to the type of tissues [3] and also in a given tissue (e.g. lymph nodes). To examine the effects of both diffusion and spatial heterogeneity, Funk et al. [8] introduced a discrete model based on (1.1a)-(1.1c). These authors adopted a twodimensional square grid with 21 × 21 sites and assumed that the virus can move to the eight nearest neighboring sites. They pointed out that the presence of a spatial structure enhances population stability with respect to non-spatial models. However, our analysis does not confirm this observation (see Section 4 below).
Recently, Wang et al. [26] generalized Funk et al.'s model. They assumed that the hepatocytes can not move under normal conditions and neglected their mobility, while viruses can move freely and their motion follows a Fickian diffusion. They proposed the following system of two ODEs coupled with a parabolic PDE for the virus: where d V is the diffusion coefficient. They assumed that the domain is the whole real line and proved the existence of traveling waves. Wang et al. [27] introduced a delay to take into account the time between infection of a target cell and the emission of viral particles [6]. They considered (1.4a)-(1.4c) in a one-dimensional interval with Neumann boundary conditions. In the spirit of the above works, we intend to study System (1.4a)-(1.4c) in a two-dimensional spatial domain (0, ℓ) × (0, ℓ) with periodic boundary conditions. There are two main situations: (i) the environment is homogeneous and, hence, all the parameters in (1.4a)-(1.4c) are constant. Therefore, the system with diffusion has the same equilibria as System (1.1a)-(1.1c); (ii) the environment is heterogeneous, therefore certain parameters become positive functions of the space variable. Then, the virus is spatially structured. For simplicity, we assume throughout the paper that only the rate α varies while the other parameters are fixed positive constants. In fact, it is biologically plausible to assume that the arrival of new CD 4 may vary according to local areas. More precisely, α is piecewise continuous and periodic in each variable with period ℓ. Then it is convenient to define the heterogeneous reproductive ratio: The sites where R 0 (x) < 1 are called sinks while the sites where R 0 (x) > 1 are called sources [8]. The paper is organized as follows. In Section 2 we are interested in the stationary problem associated with (1.4a)-(1.4c) and its non-negative equilibria. The virus equilibrium verifies the elliptic semilinear equation with periodic boundary conditions. A first issue is to define a parameter which will play the role of the bifurcation parameter R 0 in the case the latter is constant. A candidate for this role is the largest eigenvalue λ 0 of the operator: which is the linearization around V ≡ 0 of (1.6). For the reader's convenience, we recall some basic facts about two-dimensional Sturm-Liouville eigenvalue problems with periodic boundary conditions such as (1.7) and give some proofs in Appendix A.
It is clear that, whenever R 0 is a constant, λ 0 = µ V (R 0 − 1). Therefore, we distinguish two cases, depending upon the sign of λ 0 : (i) λ 0 ≤ 0: the trivial non-infected solution V u ≡ 0 is the only solution of (1.6); (ii) λ 0 > 0: (1.6) has exactly two non-negative solutions, namely the trivial noninfected solution V u ≡ 0 and the positive infected solution V i . Section 3 is devoted to the study of the evolution problem (1.4a)-(1.4c). In the case λ 0 < 0, we prove that the trivial non-infected solution (T i , I i , V i ) is asymptotically stable. Then, we turn our attention to the biologically relevant case λ 0 > 0. First, we prove the non-infected solution becomes unstable. Second, we consider (1.4a)-(1.4c) as a dynamical system and prove the existence of an universal (or maximal) attractor. Since the system is only partly dissipative, we use a result of Marion [16]. The following Section 4 is devoted to some special cases where the positive infected solution V i is stable. Particular attention is paid to the case when R 0 is a constant: in this case discrete Fourier transform can be applied.
We point out that our proof can be extended to further models in HIV literature. It is not difficult to take a logistic term into account in the T equation [22], although the steady equation (1.6) will be more involved. Adding such a term, Hopf bifurcations have been observed numerically in ODE systems (see [20]). Therefore, proving the stability of the infected solution may be, in general, challenging.
In the last section (Section 5), we consider the case when a heterogeneous environment is formed of sinks and sources alternating very rapidly, with a heterogeneous reproductive ratio R 0 ( x ε ). We determine the homogenized limiting medium as ε → 0. It is fully characterized by a constant reproductive ratio, the mean value of R 0 . Therefore, the classical approach of HIV dynamics via ODEs in a homogenous environment can be a posteriori justified in this respect.
Notation. Throughout this paper, for any ℓ > 0, we denote by L 2 the usual space of functions f : (0, ℓ) 2 → R such that f 2 is integrable. The square (0, ℓ) 2 will be simply denoted by Ω ℓ . By H k we denote the Sobolev space of order k, i.e., the subset of L 2 of all the functions whose distributional derivatives up to k-th order are in L 2 . Both L 2 and H 2 are endowed with their Euclidean norm. Finally, by H k ♯ we denote the closure in H k of the space C m ♯ of all m-th continuously differentiable functions f : R 2 → R which are periodic with period ℓ in each variable. The space H k ♯ is endowed with the norm of H k . Finally, we denote by Id the identity operator.

2.
A semilinear equation for the virus steady states: existence and uniqueness of the equilibria. We start from the system for the virus dynamics: set in Ω ℓ . Periodic boundary conditions for T, I and V are prescribed. We are interested in the existence of steady state solutions to the equations (2.1a)-(2.1c) which belong to the space L 2 × L 2 × H 2 ♯ . Clearly, any steady state solution to Problem (2.1a)-(2.1c) is a solution to the following stationary system: Hence, the function V turns out to solve the equation as a (trivial) steady state solution. We will call the triplet (T u , I u , V u ) the noninfected solution.
We are interested in studying the uniqueness of the non-infected solution in the class of all the non-negative steady state solutions to Problem (2.1a)-(2.1c). Of course, in the case when uniqueness does not hold (a situation which can actually occur, look for instance at the case when R 0 > 1 and α is constant, discussed in the introduction) we want to characterize all biological relevant steady state solutions to Problem (2.1a)-(2.1c).
For this purpose, we need to recall the following results about Sturm-Liouville eigenvalue problems in dimension two with periodic boundary conditions. Theorem 2.1. Let d and µ be, respectively, a positive constant and a bounded measurable function. Further, let A : H 2 ♯ → L 2 be the operator defined by A u = d∆u − µu for any u ∈ H 2 ♯ . Then, the spectrum of A consists of eigenvalues only. Moreover, its maximum eigenvalue λ max is given by the following formula: Finally, the eigenspace corresponding to the eigenvalue λ max is one dimensional and contains functions which do not change sign in Ω ℓ .
This is a rather classical result. Nevertheless, for the reader's convenience, we give a proof in Appendix A.
In view of Theorem 2.1, we can define the constant λ 0 to be the maximum eigenvalue of the operator ϕ → d V ∆ϕ + µ V (R 0 − 1)ϕ, which is the linearization According to (2.5), As we are going to show, the uniqueness of the non-infected steady state solution is related to the value of λ 0 . Lemma 2.2. If λ 0 ≤ 0, then the non-infected solution (T u , I u , V u ) (see (2.4)) is the only non-negative solution of (2.1a)-(2.1c).
Proof. We argue by contradiction. Let us suppose that Problem (2.2a)-(2.2c) admits another solution (T * , I * , V * ) different from (T u , I u , V u ). Then, the function V * ∈ H 2 ♯ does not identically vanish in Ω ℓ and it solves the equation (2.3). Multiplying both the sides of this equation by V * and integrating by parts in Ω ℓ , we get: Since V * does not identically vanish in Ω ℓ , the last integral term is positive, implying that Hence, the infimum in (2.5) is negative which contradicts our assumption −λ 0 ≥ 0.
Remark 2.3. From formula (2.6) it is immediate to check that, when the maximum of R 0 in Ω ℓ is less than or equal to 1, the constant λ 0 is non-positive. Hence, in this situation the non-infected solution is the only relevant steady state solution to Problem (2.1a)-(2.1c) in complete agreement with the case when R 0 is constant (see the Introduction).
The result in Lemma 2.2 is very sharp as the following theorem shows.
Theorem 2.4. Suppose that λ 0 > 0. Then, there exists a steady state solution 1c) whose components are all positive in Ω ℓ , with .
Moreover, (T u , I u , V u ) and (T i , I i , V i ) are the only steady state solutions whose components are non-negative in Ω ℓ .
Proof. It is clear, that we can limit ourselves to dealing with the equation (2.3).
Being rather long, we split the proof into two steps.
Step 1: (existence). To prove the existence of a positive solution to the equation (2.3) in H 2 ♯ , we use the classical method of upper and lower solutions. To simplify the notation, we denote by R the sup-norm of the function R 0 . We look for an upper solution v 0 of (2.3) as a constant C > 0. It is immediate to check that the best choice Note that, by Remark 2.3, R is strictly greater than 1. To determine a lower solution, in the spirit of [14, Chapt. 13, Sec. 3], we take as a candidate to be a lower solution the function v 0 = cϕ 0 with c > 0 to be fixed. Here, ϕ 0 is the unique (positive) solution to the equation Since λ 0 > 0 is fixed, the last side of the previous chain of inequalities is non-negative as soon as the function v 0 = cϕ 0 turns out to be a positive lower solution to the equation Hence, the classical method of upper and lower solutions provides us with a positive solution to the equation (2.3). It is enough to define the sequence (v n ) by recurrence in following way: for any fixed n ∈ N, v n is the unique solution in H 2 ♯ of the equation Since the function t → h(t) := t γt+µT is increasing in [0, +∞), the maximum principle (see Proposition B.1) shows that the sequence (v n ) is pointwise non-increasing. Moreover, v 0 ≤ v n ≤ v 0 for any n ∈ N. Hence, the sequence (v n ) is bounded. Moreover, by very general results for the heat equation, there exist positive constants C 1 and C 2 , independent of n, such that Since it converges pointwise in Ω ℓ , we can now infer that v n converges strongly in H 1 ♯ and weakly in H 2 ♯ to a function v ∈ H 2 ♯ which, of course, turns out to be a solution to the equation (2.3). For further details on the method of lower and upper solutions, we refer the reader, e.g., to the monograph [19].
Step 2: (uniqueness). To prove the uniqueness of the nontrivial non-negative solution to the equation (2.3), we adapt to our situation a method due to H.B. Keller [13].
Let us suppose that u is another nontrivial non-negative solution to (2.3). Then, the function w := v − u belongs to H 2 ♯ and solves the equation x, y ≥ 0.

HETEROGENEOUS VIRAL ENVIRONMENT IN A HIV SPATIAL MODEL 9
Let us denote by λ max (u, v) and λ max (u, 0), the maximum eigenvalues in L 2 of the operators Id, respectively. By Theorem 2.1, they are given by the formula Since v is non-negative and it does not identically vanish, q(u, v) ≤ q(u, 0) and there exists an open subset of Ω ℓ where q(u, v) < q(u, 0). Hence, Let us now rewrite the equation satisfied by u in the following way: Fredholm alternative implies that ϕ should be orthogonal to ζ, where by ζ we have denoted the function which spans the eigenspace associated with the eigenvalue λ max (u, 0). As it has been already remarked (see Theorem 2.1), the function ζ does not change sign in Ω ℓ . Similarly, since µ V < λ max (u, 0), ϕ is non-positive and it does not identically vanish since u does not. Hence, the function ϕ cannot be orthogonal to ζ and this leads us to a contradiction. The proof is now complete.

2.2.
Numerical illustration (steady state). In accordance with Funk et al. [8], the domain Ω ℓ is a discrete square grid with n×n sites of equal dimension ℓ/n×ℓ/n. We assume in this numerical part that all parameters vary randomly from site to site in such a way that 0.1 ≤ R 0 (x) ≤ 5.0. We deal with two cases: (i) In the first case, see Fig. 5 (left), the distribution of R 0 is as in Fig. 4a. The sources represent only 26% of the sites. We compute λ 0 = −1.70. Solving numerically the equation (2.3), we find the non-infected solution V u ≡ 0. We also represent T u according to Formula (2.4).   Note that V i is smoothly structured in space although R 0 is not. We also represent T i according to Formula (2.7).
3. Study of the dynamical system. We recall the evolution problem for the virus dynamics: set in Ω ℓ with periodic boundary conditions. We consider (3.1a)-(3.1c) as a dynamical system S (t), which has two equilibria: the non-infected trivial solution (T u , I u , V u ) and the infected, positive solution (T i , I i , V i ), the latter for λ 0 > 0 only. At first, we prove that the non-infected solution is stable for λ 0 < 0 and unstable for λ 0 > 0. By stable we mean asymptotically stable. For λ 0 > 0 the instability of (T u , I u , V u ) does not usually imply the stability of (T i , I i , V i ). Our aim is to prove the existence of a universal (or maximal) attractor which attracts all the orbits (see e.g., [25]). Since System (3.1a)-(3.1c) is only partly dissipative, we will use a result of Marion [16]. Some special cases where the stability of the infected solution is granted will be discussed afterwards.
Let us introduce the following notations: D is the domain in R 3 defined by: where M 1 , M 2 are positive constants which will be fixed throughout the proof of the next theorem. We also set for u := (u 1 , u 2 , u 3 ): Our main result is the following theorem.
Theorem 3.1. The following properties are met: ) possesses a universal attractor that is connected in H.

3.1.
Proof of (i). We begin the proof observing that the linearization (around generates an analytic strongly continuous semigroup. Indeed, L u is a bounded perturbation of the diagonal operator which is clearly sectorial since all its entries are. Hence, we can apply [15, Prop. 2.4.1(i)] and conclude that L u is sectorial. Since H 2 ♯ is dense in L 2 , the associated analytic semigroup is strongly continuous.
Let us prove that all the elements of the spectrum of L u have negative real part. In view of the linearized stability principle (see e.g., [10, Chapt. 5, Cor. 5.1.6]) this will imply that the trivial non-infected solution to Problem (3.1a)-(3.1c) is stable.
To study the spectrum of the operator L u , we fix f = (f 1 , f 2 , f 3 ) ∈ (L 2 ) 3 and consider the resolvent system where we look for a triplet of functions ϕ 1 , ϕ 2 ∈ L 2 and ϕ 3 ∈ H 2 ♯ . Suppose that λ differs from both −µ I and −µ T (which belong to the essential spectrum). Then, we can use equations (3.2a) and (3.2b) to make ϕ 1 and ϕ 2 explicit in terms of ϕ 3 . In particular, replacing the expression of ϕ 2 in terms of ϕ 3 into (3.2c) and using the very definition of the function R 0 (see (1.5)), we can transform Problem (3.2a)-(3.2c) into the equivalent equation for ϕ 3 only: Adding and subtracting −µ V R 0 from the left-hand side of (3.3), we can rewrite the equation (3.3) into the equivalent form: Note that, for any λ ∈ C, the operator A λ defined by the left-hand side of (3.4) has compact resolvent. Hence, its spectrum consists of eigenvalues only. We are going to prove that, for λ ∈ C with non-negative real part, 0 is not an eigenvalue of A λ . For this purpose, we observe that −λ 0 can be equivalently characterized as the infimum of the ratio ♯ . This shows, in particular, that for any function ψ as above. Let now ϕ 3 be a complex-valued solution to (3.4). Multiplying both the sides of such an equation by the conjugate of ϕ 3 and integrating by parts, we easily see that Taking the real part of both the sides of (3.6) and using (3.5), we obtain Since, by assumptions, λ 0 < 0 and R 0 is a positive-valued function, the only solution to the previous inequality, when Reλ ≥ 0, is the trivial function ϕ 3 ≡ 0. Hence, for these values of λ, 0 is not an eigenvalue of the operator A λ , i.e., any λ with non-negative real part belongs to the resolvent set of the operator L u .

3.2.
Proof of (ii). Again in view of the linearized stability principle, to prove the instability of non-infected solution (T u , I u , V u ) we can limit ourselves to showing that L u admits a positive eigenvalue. Hence, we are led to the study of Problem Since the parameters µ T , µ I and µ V are all positive and we are looking for positive eigenvalues λ, we can limit ourselves, as in the proof of (i), to studying the equation For any s ≥ 0, let us consider the operator d V ∆ + µ V µI R0 s+µI − 1 Id defined in H 2 ♯ . By Theorem 2.1 its spectrum consists of eigenvalues only, and the largest one is given by the following formula: As it is immediately seen, λ(0) = λ 0 . Moreover, since the function s → µI R0 s+µI is continuous, decreasing in [0, +∞) and it tends to 0 as s → +∞, the function s → λ(s) is continuous, decreasing and tends to , as s → +∞, i.e., it converges to the largest eigenvalue of the operator , it is immediate to check that the fixed point equation s = λ(s) has a positive solution. Of course, this fixed point is the positive eigenvalue of the operator L u we were looking for. This accomplishes the proof.

3.3.
Proof of (iii). Let us show that [16, Thm. 5.1] applies. In this respect the assumption α ∈ W 1,∞ (Ω ℓ ) is enough. To avoid conflict with notations, throughout the proof, T denotes time as usual, whereas the triplet (T, I, V ) is denoted by (u 1 , u 2 , u 3 ). We split the proof into several steps.
To prove the existence of a classical solution to Problem (3.7), let us fix β ∈ (1/2, 1) and introduce, for any T > 0, the space C β (T ) consisting of all functions v : (0, T ) → D ∆2 (β, ∞) such that v C β (T ) := sup t∈(0,T ] t β v(t, ·) D∆ 2 (β,∞) < +∞. Clearly, C β (T ) is a Banach space when endowed with the above norm. Moreover, D ∆2 (β, ∞) is continuously embedded into D ∆2 (β − ε, 2) for any ε ∈ (0, β). We now fix ε small enough such that θ = β − ε > 1/2. From the above results, it is immediate to infer that C β (T ) is embedded into the set of all continuous functions f : (0, T ] × Ω ℓ → R and there exists a positive constant C 1 , independent of T , such that sup Let us solve the Cauchy problem for u 1 , taking u 3 as a parameter. The (unique) solution to such a problem in C 1 ((0, T ]; L 2 ) ∩ C([0, T ]; L 2 ) is the function u 1 defined by for any t ∈ (0, T ]. If u 0,1 is bounded and continuous in Ω ℓ this result is straightforward. In the general case, we approximate u 0,1 ∈ L 2 by a sequence of smooth functions u converges to D t u 1 in C([ε, T − ε]; L 2 ) for any ε > 0. It follows that the function u 1 in (3.9) is a solution to the Cauchy problem for u 1 also in the case when u 0,1 is in L 2 .
Let us now denote by Λ 1 the operator defined in C β (T ) by the right-hand side of (3.9). A very easy computation shows that Λ 1 maps C β (T ) into C 1−θ ([0, T ]; L 2 ). Moreover, We now consider the equation for u 2 . Replacing u 1 = Λ 1 (u 3 ) in the right-hand side of this equation and using the same argument as above, we easily see that the (unique) solution in C([0, T ]; L 2 ) ∩ C 1 ((0, T ]; L 2 ) is the function u 2 defined by Let us denote by Λ 2 the operator defined in C β (T ) by the right-hand side of (3.11). Taking (3.8), (3.10a) and (3.10b) into account, one can easily show that for any v, w ∈ C β (T ). Let us now observe that (u 1 , u 2 , u 3 ) is a classical solution to Problem (3.7) if and only if u 3 is a fixed point of the operator Γ, formally defined by We are going to prove that the operator Γ is a contraction in Y β (T ) = {u ∈ C β (T ) : u C β (T ) ≤ M } provided that T and M are properly chosen. As a first step, let us prove that Γ maps Y β (T ) into itself if T, M are suitably chosen. For this purpose, we set Taking (3.12a) into account, we can estimate for any v ∈ Y β (T ). Hence, if we fix M > K u 0,3 L 2 , we can then choose T small enough such that Γ(v) C β (T ) ≤ M for any v ∈ Y β (T ). Moreover, taking (3.12c) into account, we can estimate for any v, w ∈ Y θ (T ). This estimate shows that Γ is a 1/2-contraction provided that T is sufficiently small. We can thus apply the Banach fixed point theorem and conclude that there exist T > 0 and a unique function  (u 1 , u 2 , u 3 )) enjoys the following properties: Step 2. Here, we prove that if u 0 ≥ 0 (where the inequality is meant componentwise) then the maximal defined solution u to Problem (3.7) is non-negative as well in (0, T * ). Clearly, using formulae (3.9) and (3.11) it is immediate to check that u 1 and u 2 are both non-negative whenever u 0,1 is.
Let us now consider the problem for u 3 , which we rewrite here: The heat semigroup is positive in C ♯ by the maximum principle. Since the heat semigroup in C ♯ is the restriction to C ♯ of the heat semigroup {e t∆2 } in L 2 , by density it follows that {e t∆2 } is non-negative as well. This is enough for our aims. Indeed, the function u 3 is given by the variation of constants formula (3.14) u 0,3 and u 2 being non-negative, the function u 3 is non-negative as well.
We have so proved that any solution to Problem (3.7), corresponding to an initial datum in the first octant, is confined to the first octant for any t ∈ (0, T * ). In such a case we can forget the absolute value in (3.7).
Step 3. Here, we prove that any solution u to Problem (3.7), corresponding to an initial datum u 0 ∈ D, exists for any positive time and it stays bounded. Here, For this purpose, it is convenient to introduce the so-called Svab-Zeldovich variable v = u 1 + u 2 . As it is immediately seen, the function v satisfies the Cauchy problem Since u 1 and u 2 are both positive, then (3.16) Multiplying both the sides of (3.16) by a non-negative function ϕ ∈ L 2 and integrating over Ω ℓ , one obtains that the function w(t) = Ω ℓ v(t, x)ϕ(x)dx is in C 1 ([0, T * )) and solves the differential inequality Hence, or, equivalently,

From this integral inequality, we can infer that
for any t ∈ [0, T * ) and almost any x ∈ Ω ℓ .
Since v = u 1 + u 2 and u 1 and u 2 are both non-negative, it follows that u 1 and u 2 can be estimated by the right-hand side of (3.17).
Finally, let us consider the function u 3 . From (3.14) and the above results, we can infer that Note that the function u 3 ≡ M 2 , with M 2 being given by (3.15), satisfies the previous inequality. Hence, if u 0,3 ≤ M 2 , then, the solution to Problem (3.13) is bounded from above by M 2 . With the previous choices of M 1 and M 2 , we see that the solution to Problem (3.7) which corresponds to u 0 ∈ D, stays in D for any t ∈ (0, T * ). By virtue of [15,Prop. 7.1.8], u can be extended to all the positive times.
Step 4. Here, we show that, for any u 0 ∈ D, the solution u that we have determined in the previous steps is, in fact, the unique weak solution to Problem (3.7) which belongs to L 2 ((0, T ); L 2 )× L 2 ((0, T ); L 2 )× L 2 ((0, T ); H 1 ♯ ) for any T > 0. Even if the following arguments are standard, for the reader's convenience we go into details.
Let us now consider the Cauchy problem for u 3 (i.e., problem (3.13)). Since u 2 is Lipschitz continuous in [0, T ] with values in L 2 , for any T > 0, by [15,Thm. 4 . Now, we turn back to the equations for u 1 and u 2 and conclude that D t u 1 and D t u 2 are in C([0, T ]; L 2 ) for any T > 0, this implying that u 1 and u 2 are in C 1 ([0, T ]; L 2 ). Hence, any weak solution to Problem (3.7) with data in D is such that u 1 , u 2 ∈ C 1 ([0, T * ); L 2 ) and u 3 ∈ C([0, T * ); L 2 ) ∩ C 1 ((0, T * ); L 2 ) ∩ C((0, T * ); H 2 ♯ ). Since we have proved uniqueness of the solution in this class of functions, uniqueness of the weak solution follows as well.
Since for non-negative solutions the Cauchy problem (3.7) coincides with problem (3.1a)-(3.1c), we have, thus, established the following: Step 3. We are now in a position to apply [16, Thm. 5.1]. The set Obviously the matrix G(·, u 3 ) has positive eigenvalues whenever u 3 ≥ 0, which are bounded from below by positive constants. Hence, for any ξ ∈ R 2 . Thus, condition (4.6) in [16] is satisfied. The proof of Theorem 3.1 is completed.

4.
A gamut of some special cases.

Numerical illustration (evolution).
We continue the discussion of Subsection 2.2 in the framework of the evolution problem (3.1a)-(3.1c). In the first case (λ 0 < 0), only the non-infected steady state V u ≡ 0 exists. We solve (3.1a)-(3.1c) under particular initial conditions: we start the infection at the center of the grid with an inoculum of one viral unit, assuming that T and I are at their uninfected steady state. One observes that the virus vanishes very rapidly and the target cells return to their initial level (see Fig. 6 left) in accordance with the stability of the uninfected equilibrium V u . In the second case (λ 0 > 0), two equilibria exist, V u and V i . Starting with the same initial conditions, the virus population grows while the population of target cells decreases. Both of them achieve an equilibrium corresponding to the positive infected solution (see Fig. 6 right).   Fig. 4a, λ 0 < 0 (no infection); right: R 0 as in Fig. 4b, λ 0 > 0 (infection). The solid, the long and the short dashed lines depict respectively the dynamics at (10,10), at the site of viral inoculum (20,20) and at a border site (1,1). Here d V = 1, ℓ = 1.
We are now going to review some particular cases where the stability of the infected solution is granted.

4.2.
Homogeneous environment with diffusion. We consider the case when α is a constant such that R 0 (see (1.3)) is a constant, as well, greater than 1, together with d V > 0 as in [8].
The linearization around (T i , I i , V i ) (see (1.2)) of Problem (3.1a)-(3.1c) is associated with the linear operator The same arguments as in the proof of Theorem 3.1(i) show that the realization L i of the operator L i in (L 2 ) 3 with domain D(L i ) = L 2 × L 2 × H 2 ♯ generates an analytic strongly continuous semigroup. We are going to determine its spectrum. For this purpose we use the discrete Fourier transform.
We claim that the spectrum of the operator L i is given by where σ k is the spectrum of the matrix To check the claim, let us observe that, if a function where f i,k denotes the k-th Fourier coefficient of the function f i (i = 1, 2, 3). Clearly, any eigenvalue of M k (k = 0, 1, . . .) is an eigenvalue of L i . Therefore, On the other hand, if λ ∈ σ k for any k = 0, 1, . . ., then all the coefficients (v 1,k , v 2,k , v 3,k ) are uniquely determined through the formulae Note that, if λ differs from both −µ I and −µ T R 0 , then Hence, for any λ / as k → +∞. It follows that the sequences {v 1,k }, {v 2,k } and {λ k v 3,k } are in ℓ 2 . This shows that the series whose Fourier coefficients are v 1,k , v 2,k and v 3,k , respectively, converge in L 2 (the first two ones) and in H 2 ♯ (the latter one). The inclusion follows. We now observe that a straightforward computation shows that λ = −µ T R 0 is in the essential spectrum of L i . Also λ = −µ I belongs to the essential spectrum of L i and, in the case when µ I = µ T , it belongs also to the point spectrum. The set equality (4.1) is proved.
Clearly σ k has three eigenvalues (counted with their multiplicity), either all real, or one real and two complex conjugates. Routh-Hurwitz criterion enables us to determine whether the elements of σ k have negative real parts. The latter holds if and only if b k , d k and b k c k − d k are positive, which is clearly true whenever R 0 > 1. Remark 4.1. As it is easily seen σ 0 is the spectrum of the operator  which is associated with the linearization at (T i , I i , V i ) of Problem (1.1a)-(1.1c).
Since, as we have already remarked, the spectrum of L i is the union of the sets σ k (k ∈ N ∪ {0}) and the points −µ I , −µ T R 0 , the scenario is one of the following: (a) the diffusion does not improve the stability of the solution (T i , I i , V i ). Therefore, the stability issue is identical to that of the ODE system (d V = 0); (b) the diffusion worsen the stability of the solution (T i , I i , V i ). Therefore, we are unable to confirm Funk et al. [8], who pointed out that the presence of a spatial structure enhances population stability with respect to nonspatial models. Only some smoothing effect can be credited to the diffusion.

4.3.
Death rates µ I ≤ µ T , λ 0 > 0. This case leads to a mathematically interesting framework although it has little biological relevance, since µ I > µ T in the literature (e.g., µ T = 0.01, µ I = 0.39 in [5]). For the latter reason we will not elaborate the case extensively.
Using the Svab-Zeldovich variable v = u 1 + u 2 we can transform Problem (3.1a)-(3.1c) into the following equivalent one for the unknowns u 1 , v and u 3 : 2) It is not difficult to see that the mapping u 3 → u 1 is non-increasing, so is the mapping u 1 → v thanks to the hypothesis µ T − µ I > 0. Hence, the mapping u 3 → v is non-decreasing. Finally, the mapping u 3 → u 2 = v − u 1 is non-decreasing. Based on these observations, following [19] it is possible to construct two sets of monotone sequences which converge to the solution of (4.2). These sequences start respectively from upper and lower solutions defined as in Section 2, to stay away from the trivial solution. It is well-known that a solution of an evolution problem constructed via such a monotone sequence scheme, with suitable initial conditions between upper and lower solutions, achieves a stable equilibrium (see [19]). Therefore, the infected solution is asymptotically stable. Numerical computations in the phase plan (see Fig. 7) illustrate the difference in the virus dynamics when µ I < µ T (monotonicity) and µ I > µ T (spirals).

4.4.
Quasi-steady problem. In this part we assume that T and I are at their equilibrium. In such a case, (3.1a)-(3.1c) reads which is equivalent to the scalar parabolic equation for V only, with periodic boundary conditions: The latter is the natural evolution problem associated with (2.3). It is clear that (4.3) has the same non-positive equilibria, namely V u = 0 and V i > 0, in the case when λ 0 > 0. The stability of V i can be proved according to [10,Sec. 5.3] by constructing a Lyapunov function.

5.
Homogenization. This section is concerned with the case wherein the environment is heterogeneous and is formed of rapidly alternating sinks and sources. For a fixed integer k, we imagine that Ω ℓ is divided into a network of k 2 periodic squares Ω ε , where ε = ℓ/k. The heterogenous reproductive ratio R 0 will depend upon ε, see Fig. 8b. The idea of homogenization is to let ε → 0 and find the equivalent homogenized medium. Therefore, such a heterogenous environment can be replaced by its homogenized limit for easier computations and analysis. More precisely, we introduce a normalized periodic function R 0 as a function of the variable y = (y 1 , y 2 ), of period (1, 1), and we define: Remark 5.1. x = (x 1 , x 2 ) is the macroscopic variable while y = (y 1 , y 2 ) is the microscopic one.
We consider the problem: on Ω ℓ with periodic boundary conditions as above. The idea is to find the limiting homogenized equation as ε → 0. We start with the following lemma: Lemma 5.2. Let R 0 : R 2 → R be a periodic (with period 1 in each variable), piecewise continuous function. For x ∈ R 2 and ε > 0, set R ε 0 (x) = R 0 ( x ε ). Then, the following properties are met.
(i) R ε 0 tends to M y (R 0 ) as ε → 0, weakly in L p loc (R 2 ) for any p ≥ 1; (ii) For any ε > 0, set Then, λ ε 0 → µ V (M y (R 0 ) − 1) as ε → 0. Proof. Property (i) follows straightforwardly from e.g., [12, p. 5]. Thanks to (i), one can take the limit as ε → 0 in the right-hand side of (5.2) and show that λ ε 0 tends to the largest eigenvalue of on Ω ℓ with periodic boundary conditions. Let us prove this claim. As a first step, we observe that there exist two constants C 1 and C 2 such that C 1 ≤ λ ε 0 ≤ C 2 , since the function R ε 0 is bounded. Next, λ ε 0 is the largest eigenvalue of the Sturm-Liouville eigenvalue problem with periodic boundary conditions we documented in Theorem 2.1, associated with the eigenfunction ψ ε . We may assume that Ω ℓ (ψ ε ) 2 dx = 1. As we pointed it out in Theorem 2.1, ψ ε does not change sign. It is clear that ψ ε is bounded in H 2 ♯ . Then, there exists an infinitesimal sequence {ε n } such that λ εn 0 → λ 0 0 , ψ εn → ψ 0 weakly in H 2 ♯ , strongly in H 2−η ♯ and (hence) uniformly in Ω ℓ . Note that Ω ℓ (ψ 0 ) 2 dx = 1 and ψ 0 does not change sign. Since d V ∆ψ εn − µ V (1 − R εn 0 )ψ εn = λ εn ψ εn , it is not to difficult to pass to the limit in the above equation as n → +∞ in the distributional sense and see that with periodic boundary conditions. Therefore, λ 0 0 is an eigenvalue of the operator d V ∆ − µ V (1 − R ε 0 )Id, associated with the eigenfunction ψ 0 . Since ψ 0 does not change sign, λ 0 0 is the largest eigenvalue of (5.3). Obviously, ψ 0 is a constant, hence λ 0 0 is explicit. Finally, checking that all the sequence λ ε 0 converges to λ 0 0 , as ε → 0, is an easy task.
Proof. As a first step, we observe that, since M y (R 0 ) > 1, from Lemma 5.2(ii) it follows that there exist δ, ε 0 > 0 such that λ ε 0 > δ for ε ∈ (0, ε 0 ], and (5.1) has a unique positive solution V ε . Clearly, ∆V ε is bounded in L ∞ (Ω ℓ ) and this implies that V ε is bounded in H 2 ♯ . Arguing as in the proof of Lemma 5.2, one can extract an infinitesimal sequence {ε n } such that V εn converges strongly in H 2−η to a limit V 0 ∈ H 2 ♯ , which verifies the equation Because of the periodic boundary conditions, Equation (5.4) has only constant solutions. Therefore, and it only remains to prove that V 0 is not the trivial solution. With obvious notations, we recall (see the proof of Theorem 2.4) that V ε ≥ v ε 0 = c ε ϕ ε 0 where ϕ ε 0 is the positive eigenfunction associated with the largest eigenvalue λ ε 0 and Since max ϕ ε 0 = 1, it is clear that V ε remains bounded away from 0 whenever ε ≤ ε 0 . Finally, checking that V ε itself converges to V 0 as ε → 0 is immediate. This concludes the proof.
It is easy to compute the mean value M y (R 0 ) = 1.16 and the homogenized viral density V 0 = 16.1. Fig. 9 shows how the virus V ε at ε = 0.1 (left) oscillates slightly around its homogenized limit V 0 (right).