Two-dimensional stability analysis in a HIV model with quadratic logistic growth term

We consider a Human Immunodeficiency Virus (HIV) model with a logistic growth term and continue the analysis of the previous article [6]. We now take the viral diffusion in a two-dimensional environment. The model consists of two ODEs for the concentrations of the target T cells, the infected cells, and a parabolic PDE for the virus particles. We study the stability of the uninfected and infected equilibria, the occurrence of Hopf bifurcation and the stability of the periodic solutions.


Introduction
Over the past thirty years, there has been much research in the mathematical modeling of Human Immunodeficiency Virus (HIV), the virus which causes AIDS (Acquired Immune Deficiency Syndrome). The research directions have been twofold: (i) the epidemiology of AIDS; (ii) the immunology of HIV as a pathogen. We are interested in the latter approach.
The major target of HIV infection is a class of lymphocytes, or white blood cells, known as CD4 + T cells. When the CD4 + T-cell count, which is normally around 1000 mm −3 , reaches 200 mm −3 or below in an HIV-infected patient, then that person is classified as having AIDS.
Mathematical models have been proved valuable in understanding the in vivo dynamics of the virus. A gamut of models have been developed to describe the immune system, its interaction with HIV, and the decline in CD4 + T cells. They have contributed significantly to the understanding of HIV basic biology.
Recently, the effect of spatial diffusion has been taken in account in HIV modeling. Funk et al. [7] introduced a discrete model: they adopted a two-dimensional square grid with 21 × 21 sites and assumed that the virus can move to the eight nearest neighboring sites. K. Wang et al. [18] generalized Funk's model. They assumed that the hepatocytes can not move under normal conditions and neglected their mobility, whereas virions can move freely and their motion follows a Fickian diffusion. In [1], two of the authors considered a two-dimensional heterogenous environment: the basic reproductive ratio is generalized as an eigenvalue of some Sturm-Liouville problem. Furthermore, in the case of an alternating structure of viral sources, the classical approach via ODE systems is justified via a homogenized limiting environment.
In this article, we consider a HIV model which takes the viral diffusion into account in a homogeneous two-dimensional environment, and includes a quadratic logistic growth term as previously proposed in [16,17] to consider the homeostatic process for the CD4 + T-cell count. The model reads: The spatial domain is denoted by Ω ℓ = (0, ℓ) × (0, ℓ), periodic boundary conditions are prescribed for V . Since the system (1.1)-(1.3) defines a dynamical system or semiflow, we will also use the abstract notation X(t) = (T (t), I(t), V(t)).
Our aim is to continue the analysis of the previous paper [6], where we studied the system (1.1)-(1.3) when d V = 0. We refer to [6] for an extended introduction to the biological issues. In brief, we recall that T and I denote the respective concentrations of uninfected and infected CD4 + T cells. The concentration of free virus particles, or virions, is V (for the sake of simplicity, we call V the virus). In (1.1), r is the average specific T-cell growth rate obtained in the absence of population limitation. The term 1 − T/T max shuts off T-cell growth as the population level T max is approached from below. Here µ T is the natural death rate of CD4 + T cells, the term γVT models the rate at which free virus infects a CD4 + T cell. The infected cells die at a rate µ I and produce free virus during their life-time at a rate N. In addition, µ V is the death rate of the virus. According to the literature (see e.g., Table 1 or [2] where µ T = 0.01, µ I = 0.39), we assume the following biologically relevant hypothesis: µ I > µ T .
(1.4) Note that the quantity r − µ T , the net T-cell proliferation rate, needs not to be positive (see [16, p. 86] Production rate for uninfected CD4 + T cells 1.5 day −1 mm −3 γ Infection rate of uninfected CD4 + T cells 0.001 day −1 mm 3 T max Maximal population level of CD4 + T cells at 1500 mm −3 which the CD4 + T-cell proliferation shuts off µ T Death rate of uninfected CD4 + T-cell population 0.1 day −1 µ I Death rate of infected CD4 + T-cell population 0.5 day −1 µ V Clearance rate of free virus 10 day −1 Derived variable T 0 CD4 + T-cell population for HIV negative persons mm −3 From a mathematical viewpoint: (i) N > 0 and r ≥ 0 are parameters; (ii) the quantities (with associated dimension) α, γ, µ I , µ T and µ V are fixed positive numbers throughout the paper; (iii) T max is a large perturbation parameter, larger than any finite combination of α, γ, µ I , µ T and µ V of the same dimension (mm −3 ). In particular, this hypothesis contains the condition T max > α/µ T of [16, p. 85].
The paper is organized as follows: in Section 2, we prove that System (1.1)-(1.3) admits, for any value of the parameters r and N, the uninfected steady state X u = (T u , 0, 0) and that, in a region of the space of parameters, there exists also another steady-state solution, the so-called infected steady state X i = (T i , I i , V i ), where T i , I i and V i are positive. In the parameter space, we define the regions U and I (this latter being the region where the infected steady-state exists), respectively for uninfected (the reproductive ratio is such that R 0 < 1) and infected (R 0 > 1). We recall that the basic reproductive ratio denotes the average number of infected T cells derived from one infected T cell ( [4]). We prove that the uninfected steady state is asymptotically stable in U, and unstable in I.
In [6], we have exhibited an unbounded subdomain P in I in which the positive infected equilibrium becomes unstable whereas it is asymptotically stable in the rest of I. In this unstable region, the levels of the various cell types and virus particles oscillate, rather than converging to steady values. This subdomain P may be biologically interpreted as a perturbation of the infection by a specific or unspecific immune response against HIV. In Section 3, we consider the linearization around Jacobian matrix L i . A modal expansion of the resolvent equation enables us to construct a finite number of subdomains P k (k = 0, . . . , K 2 ) in I, that form a monotone non-increasing sequence (for the inclusion) with P 0 = P. It turns out that the infected equilibrium X i = (T i , I i , V i ) is asymptotically stable for (N, r) ∈ I \ P and unstable in the interior of P. Therefore the stability issue is governed by the 0-th mode, hence similar to the case without viral diffusion (d V = 0). As a matter of fact, we are unable to confirm Funk et al. [7], who suggested that the presence of a spatial structure enhances population stability with respect to non-spatial models (see also [1]).
In Section 4, we take the logistic parameter r as bifurcation parameter and prove the existence of Hopf bifurcations at the boundary ∂P. Since the system is only partially dissipative, the resolvent operator associated to the realization L i of L i is not compact and therefore the proof demands more attention: it relies on the analyticity of the semigroup exp(tL i ) (see e.g., [10], [15]). Next, we perform a nonlinear analysis at the Hopf points via the Center Manifold theorem. It turns out that the bifurcating periodic solutions are independent of the space variables.
Numerical illustrations are presented in Section 5. Finally, for the sake of completeness, we recall in an Appendix some basic facts about the eigenvalues of the two-dimensional Laplace operator with periodic boundary conditions and some Sturm-Liouville operators.
Notation. For any ℓ > 0 we denote by L 2 the usual space of square-integrable functions f : (0, ℓ) 2 → R. The square (0, ℓ) 2 will be simply denoted by Ω ℓ . By H k ♯ (k = 1, 2, . . .) we denote the closure in H k (the subset of L 2 of all the functions whose distributional derivatives up to k-th order are in L 2 ) of the space C k ♯ of all k-th continuously differentiable functions f : R 2 → R which are periodic with period ℓ in each variable. The space H k ♯ is endowed with its Euclidean norm. If X is any of the previous spaces, we write X C to denote the space of complex-valued functions f such that Re f and Im f are in X. The norm in X C is defined in the natural way: If v is a vector of C 3 , we denote by v 1 , v 2 and v 3 its components. Similarly, if f is a function defined in Ω ℓ with values in R 3 (resp. C 3 ), we denote by f 1 , f 2 and f 3 its components. If the components of the vector v are complex numbers, we denote by v the vector whose components are the conjugates of the components of v. The Euclidean inner product in L 2 for any f, g ∈ L 2 C × L 2 C × L 2 C . Finally, we denote by Id the identity operator, and by (·) + the positive part of the number in brackets.

Equilibria
In this section we are devoted to determine the non-negative equilibria of System (1.1)-(1.3), i.e., the solutions (T, To state the first main result of this section, let us introduce some functions and a few notation.
By X u and X i we denote, respectively, the function whose entries T u , I u and V u are given by and the function whose entries are given by T i , I i and V i , where We further introduce two sets which will play a fundamental role in all our analysis, namely the uninfected and infected regions U and I in the parameter space, which are defined by is the reproduction ratio. The interface R 0 (N, r) = 1 between the two regions U and I is the graph of the mapping which is decreasing by virtue of the condition T max > αµ −1 T . Its image is the interval Inverting the roles of N and r, it is useful to define the inverse mapping: As it has been stressed in the introduction, throughout the paper we assume that To prove the following theorem we assume also that where T 0 max is fixed and large, and depends only on the quantities in brackets.
Theorem 2.1. The following properties are satisfied: (i) in U the uninfected steady state X u = (T u , I u , V u ) is the only non-negative equilibrium; (ii) in I there exist two non-negative equilibria, respectively X u and the infected steady Replacing the expression of I given by (2.2) in (2.3) and, then, using (2.7), we obtain the following self-contained nonlinear equation for V: As it is easily seen any solution to (2.8) in H 2 ♯ leads to a solution to System (1.1)-(1.3). Moreover, from any non-negative solution to Equation (2.8) we can obtain an equilibrium to System (1.1)-(1.3) will all the components non-negative in Ω ℓ . Hence, we can limit ourselves to looking for non-negative solutions V ∈ H 2 ♯ to Equation (2.8). Clearly, (2.8) admits the trivial function V ≡ 0 as a solution. This solution leads to the equilibrium X u .
(ii) Let us look for other positive constant solutions to Equation (2.8). We are thus lead to look for solutions to the equation Φ(V) − µ V V = 0 which are non-negative.
A straightforward computation reveals that V i is the unique solution to such an equation. Moreover, for any fixed r > 0, V i is positive if and only if N > N crit (r) (see (2.4)) i.e., if and only if (r, N) ∈ I. In this case, replacing V = V i into (2.7) and (2.2), we immediately conclude that the function X i is an equilibrium of System (1.1)-(1.3).
(iii) Showing that T u < T max is just an exercise. On the hand, the inequality T i < T max in I follows from the definition of T i observing that, if (r, N) ∈ I, then (iv) To prove that X u and X i are the only equilibria of System (1.1)-(1.3) with all the components being non-negative, we adapt to our situation a method due to H.B. Keller [13]. We argue by contradiction. We suppose that X = (T, I, V) is a solution to System (2.1)-(2.3) with V non-negative in Ω ℓ and not identically vanishing, and such that V V i . Let us set W := V − V i . Since both V and V i are solutions to (2.8), clearly W ∈ H 2 ♯ and solves the equation for any x, y ≥ 0.
Let λ max (V, V i ) and λ max (V, 0) denote the maximum eigenvalues in L 2 of the operators d V ∆ + Λ(V, V i )Id and d V ∆ + Λ(V, 0)Id, respectively. By Corollary A.2, we know that We now observe that Here, we have taken advantage of the Sobolev embedding theorem to infer that V ∈ H 2 ♯ is continuous in Ω ℓ . Note that the constant C is positive. From this remark we can easily infer that Since W satisfies (2.9) and it does not identically vanish in Ω ℓ , µ V ≤ λ max (V, V i ) < λ max (V, 0).
To get to a contradiction, we now rewrite the equation satisfied by V in the following way: Fredholm alternative implies that Z should be orthogonal to the function ψ which spans the eigenspace associated with the eigenvalue λ min (V, 0) of the operator d V ∆ + Λ(V, 0)Id. But this can not be the case. Indeed, by Corollary A.2 the function ψ does not change sign in Ω ℓ . Moreover, since V is non-negative in Ω ℓ and it does not identically vanish in Ω ℓ and, in addition, µ V < λ max (V, 0), Z is non-positive and it does not identically vanish in Ω ℓ . Hence, ψ is not orthogonal to Z.

Stability of the equilibria
In this section we are going to study the stability of the equilibria X u and X i . We begin by studying the stability of the uninfected equilibrium X u . (i) in the domain U the uninfected equilibrium X u is asymptotically stable; (ii) in the domain I the uninfected equilibrium X u is unstable.
Proof. To avoid cumbersome notation, throughout the proof we do not stress explicitly the dependence of the functions and operators on r and N. We prove the statement showing that the linearized stability principle (see e.g., [11, Chpt. 5, Cor. 5.1.6]) applies to our situation. For this purpose, we begin by observing   Table 1, N crit decreases 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 [14, Prop. 2.4.1(i)] and conclude that L u is sectorial. Since H 2 ♯,C is dense in L 2 C , the associated analytic semigroup is strongly continuous.
To complete the proof, we need to study the spectrum of the operator L u . We fix λ ∈ C and consider the resolvent equation If λ −µ I we can use the second equation to write I in terms of V. Substituting it in the last equation we get the following self-contained equation for V: We recall that the spectrum of the realization A of the Laplace operator in L 2 C , with H 2 ♯,C as a domain consists of eigenvalues only and it is given by for any λ with non-negative real part, and Equation (3.2) is uniquely solvable.
We can now uniquely determine I ∈ L 2 C from the second equation in (3.1). Finally, from the first equation in (3.1), observing that we can uniquely determine T ∈ L 2 C . We have so proved that any λ ∈ C with non-negative real part is in the resolvent set of the operator L u , if (N, r) ∈ U. In view of the linearized stability principle this implies that the trivial uninfected solution to System (1.1)-(1.3) is asymptotically stable.
The issue of the stability of the infected solution is obviously much more complicated. For notational convenience we sort the eigenvalues of the realization A of the Laplacian in L 2 C , with H 2 ♯,C as a domain (see Appendix A), into a non-increasing sequence {−λ k }. Similarly, we denote byẽ k the eigenfunction associated with the eigenvalue λ k . This allows ut to expand any function f ∈ L 2 ♯,C into the Fourier series where f k denotes the k-th Fourier coefficient (with respect to the system (ẽ k )) of f . As it is observed in the proof of Theorem A.1, λ 0 is simple and all the other eigenvalues are semisimple, and their multiplicity can be computed explicitly.
We are going to prove the following results.

Theorem 3.2.
Under the hypothesis (2.6), it holds: 3.1. The resolvent equation. As in the proof of Theorem 3.1, we do not stress explicitly the dependence on r and N of the operators and the sets that we consider in what follows. For (N, r) ∈ I, the linearization around generates an analytic strongly continuous semigroup. Moreover, the spectrum of L i is given by where, for any k ∈ N, σ k is the spectrum of the matrix Proof. The same arguments as in the proof of Theorem 3.1 show that L i generates an analytic strongly continuous semigroup in (L 2 C ) 3 . Let us determine its spectrum. For this purpose we use the discrete Fourier transform.
On the other hand, if λ σ k for any k = 0, 1, . . ., then all the coefficients where and as k → +∞. Thus, the sequences {v 1,k }, {v 2,k } and {λ k v 3,k } are square-summable. This shows that the series whose Fourier coefficients are v 1,k , v 2,k and v 3,k , respectively, converge in L 2 C (the first two ones) and in H 2 ♯,C (the latter one).
then the series having v 1,k , v 2,k and v 3,k as Fourier coefficients do not, in general, converge in L 2 (the first two ones) and in H 2 ♯ (the latter one). Hence, these values of λ belong to the essential spectrum of L i . Thus, (3.3) is proved.
The Routh-Hurwitz criterion enables us to determine whether the elements of σ k have negative real parts. The latter holds if and only if d 1,k , d 3,k and the leading Hurwitz determinant D 2,k = d 1,k d 2,k − d 3,k are positive. The case k = 0 corresponds to the system of ODEs considered in [6].
Obviously, d 1,k > 0. As far as d 3,k is concerned, we remark that d 3,k (N, r) > d 3,0 (N, r), which is positive in I and vanishes at N crit (r).
For (N, r) ∈ I, we compute the Hurwitz determinant D 2,k (N, r) and we get where . In (3.13), as a function of λ k , the denominator vanishes at which is positive and generically does not meet any of the λ k 's, k ≥ 1. There are two cases (see Fig. 2): , which in turn is of the sign of a k N 2 + b k N + c k . The coefficients a k , b k and c k read: Let us examine the signs of these coefficients.
(i) As a function of λ k , the coefficient a k vanishes at Clearly, Λ 2 is positive thanks to (2.6) and, as Λ 0 , generically does not meet any of the λ k 's for k ≥ 0. Then, a k is positive if 0 ≤ λ k < Λ 2 and negative otherwise. (ii) b k < 0 due the hypothesis (2.6). (iii) c k > 0 under the biologically relevant hypothesis µ I > µ T (see (2.5)). Next we compute: Again, thanks to the hypothesis µ I > µ T , δ k is always positive. Therefore, the roots of a k N 2 + b k N + c k = 0, namely the ones of ∆ k (N) = 0, are: (3.14) Remark 3.4. Since B k (N 0,k ) = 0 and A k (N 0,k ), C k (N 0,k ) are both positive, a k (N 0,k )) 2 + b k N 0,k + c k < 0. This provides us with a useful information regarding the position of N 0,k with respect to N 1,k and N 2,k , according to (3.14), i.e., N 1,k < N 0,k < N 2,k if a k > 0 and N 0,k < N 2,k or N 0,k > N 1,k if a k < 0. In particular, N 0,k can meet neither N 1,k nor N 2,k .
Now we are in a position to begin the discussion, depending on the position of λ k . Let us distinguish three cases. Case I: 0 ≤ λ k < Λ 2 . In this situation a k > 0, b k < 0 and c k > 0 and k ranges in a finite set of indexes. It is an extension of the case k = 0, see [6]. It follows from Remark 3.4 that B k vanishes between N 1,k and N 2,k which are both positive. In particular, for T max large enough (as we are assuming), N 1,k > µ V /(γT max ). Indeed, if T max is large enough. It follows that either µ V /(γT max ) < N 1,k or µ V /(γT max ) > N 2,k . But as it is immediately seen, N 2,k > µ V /(γT max ). Indeed, as T max → +∞. Hence, µ V /(γT max ) < N 1,k as it has been claimed. We consider four subcases depending on the position of N with respect to N 1,k and N 2,k .
admits the two real and positive roots: Observe that as T max → +∞. Consequently, the Hurwitz determinant D 2,k (N, r) is positive for r crit (N) ≤ r < r 1,k (N) and r > r 2,k (N), it vanishes at r = r 1,k (N) and r = r 2,k (N), and is negative for r 1,k (N) < r < r 2,k (N). (iv) Assume N ∈ {N 1,k , N 2,k }. In such a case, ∆ k (N j,k ) = 0 and the polynomial A k r 2 + B k (N j,k )r + C k (N j,k ) has the double root r k (N j,k ) = −B k (N j,k )/2A k . However, this solution makes sense only if B k (N j,k ) < 0. Hence, only the case N = N 2,k is relevant, and we have D 2,k (N 2,k , r) > 0 for r > r crit (N 2,k ) except at r(N 2,k ) = −B k (N 2,k )/2A k , where it vanishes.
We are now in a position to define the subdomain P k of I by see (3.14) and (3.15), and at N = N 2,k , r 1,k (N 2,k ) = r 2,k (N 2,k ) = −B k (N 2,k )/2A k . In the domain I, d 1,k and d 3,k are positive, and D 2,k (N, r) is positive except in P k . More precisely, the Hurwitz determinant D 2,k (N, r) is negative in the interior of P k and it vanishes on the boundary of P k . Case II: Λ 2 < λ k < Λ 0 . Now, a k < 0, b k < 0, c k > 0. Hence, N 2,k < 0, N 1,k > 0. Since N 0,k > 0, it holds that N 0,k > N 1,k according to Remark 3.4. Moreover, as it is immediately seen, N 0,k > (γT max ) −1 µ V . There are two possibilities: (i) (N, r) ∈ I satisfies N < N 0,k . Then, B k (N) > 0. Therefore, A k r 2 + B k (N)r + C k (N) > 0 for any r ≥ r crit (N) since the coefficients are all positive. It thus follows that D 2,k (N, r) > 0. (ii) (N, r) ∈ I satisfies N ≥ N 0,k . Then, N > N 1,k and a k N 2 + b k N + c k < 0. Therefore ∆ k (N) < 0 and A k r 2 + B k (N)r + C k (N) has the sign of A k which is positive, so D 2,k (N, r) > 0.
Case III: λ k > Λ 0 . Here, N 0,k < 0 and, therefore, B k (N) > 0 for all N > 0. The conclusion is the same as in Case II (i). We summarize our results in the following proposition.

Remark 3.6. To give an idea, with the numerical values of
To conclude this subsection we prove the following proposition which gives a much clearer picture of how the sets P k are ordered in the space of the parameters. Proposition 3.7. Let K 2 be as in the statement of Proposition 3.5. Then, the following set inclusions hold: Proof. To begin with we claim that N 2,k < N 2,k+1 (see (3.14)) for any k = 0, . . . , K 2 − 1. To prove the claim we observe that δ k = δ(λ k ), where We compute the derivative of the function δ and get Under hypothesis (2.6) this function is positive and, consequently, k → δ k is non-decreasing.
Similarly, a k = a(λ k ) and b k = b(λ k ), the functions a and b being strictly decreasing. Hence, the sequences {a k } and {b k } are non-increasing. Since a k > 0 and b k < 0 we now easily get the claim.

Proof of Theorem 3.2.
The proof follows from Propositions 3.5, 3.7, Routh-Hurwitz criterion and the linearized stability principle.
By Proposition 3.7 it holds that P k P 0 , k = 1, 2, . . .. Hence, we conclude that D 2,k (N, r) > 0 for any k ∈ N, if (N, r) ∈ I \ P 0 . Since the other two Hurwitz determinants are positive in the whole of I, it follows from the Ruth-Hurwitz criterion that, if (N, r) ∈ I \ P, then all the element of k∈N σ k have negative real part. Hence, Re σ(L i ) < 0 (see (3.4)). It remains to invoke the linearized stability principle as in the proof of Theorem 3.1.
(ii) The instability of X u can be deduced from [6] which deals with System (1.1)-(1.3) in the case when d V = 0 and shows that, in this situation, the infected equilibrium X i is unstable.

Hopf bifurcation and instability
For fixed N > 0 we take the logistic parameter r > r crit (N) as a bifurcation parameter. We recall that at fixed (N, r) ∈ I, System (1.1)-(1.3) has two equilibria: the uninfected trivial solution X u (N, r) and the infected, positive solution X i (N, r). At X i (N, r), the Jacobian matrix is L i = L i,N,r , see (3.3). As we already observed in Proposition 3.3, the realization L i,N,r of the operator L i,N,r in ( generates an analytic strongly continuous semigroup that we denote by e tL i,N,r . In this section we are interested in proving that Hopf bifurcation occurs on the boundary of the set P (i.e., at the points (N, r 1,0 (N)) and (N, r 2,0 (N)) with N ≥ N 2,0 , where r 1,0 (N), r 2,0 (N) and N 2,0 are given by (3.14) and (3.15)) and in analyzing the stability of the bifurcated periodic solutions.
Note, that for T max large, Hence, N 2,0 is positive if T max is large enough, let us say, if T > T 1 max (µ T , µ I , µ V , α, γ) > T 0 max (µ T , µ I , µ V , α, γ) (see (2.6)). We assume hereafter that Here, differently from the previous sections, to avoid confusion we stress explicitly the dependence of the operators, numbers and sets that we consider on r. We do not stress the dependence on N since in the following discussions only the parameter r varies, N is (arbitrarily) fixed. In particular, we simply write r 1 and r 2 instead of r 1 (N) and r 2 (N).
Proof. We limit ourselves to considering the case when r = r 1 , the case r = r 2 being completely similar.
For r in some neighborhood of r 1 , we set u = X − X i (r), s = r − r 1 and write System where Clearly, by the Sobolev embedding theorem, F is a smooth function defined in L 2 ×L 2 ×H 2 ♯ . Note that the derivative F u (0, 0) is the operator L i,r 1 in Proposition 3.3. More precisely, By Proposition 3.3, operator L i,r 1 is the generator of a strongly continuous analytic semigroup in (L 2 C ) 3 . Let us prove that σ(L i,r 1 ) consists of eigenvalues with negative real part and a pair of purely imaginary and conjugate eigenvalues λ 1 (r 1 ) and λ 2 (r 1 ), which are simple eigenvalues and satisfy the transversality condition. Once checked, these properties will yield the assertion in view of [14,Thm. 9.3.3] (which deals with fully nonlinear problems but, of course, it applies also to the semilinear case).
Being rather long, we split the proof into four steps.
Step 1. Here, we prove that σ(L i,r 1 ) consists of eigenvalues with negative real part and a pair of purely imaginary and conjugate eigenvalues. For this purpose, we observe that, since P k is properly contained in P 0 for any k = 1, . . . , K 2 (see Proposition 3.7), the pair (N, r 1 (N)) belongs to I \ P k for any k = 1, . . . , K 2 . Therefore, from the results in Subsection 3.2 and Proposition 3.3, it follows that σ k,r 1 is contained in the halfplane {λ ∈ C : Re λ < 0}.
Since 1≤k<K 3 σ k,r 1 consists of finitely many eigenvalues with negative real part, up to replacing M with a smaller constant if needed, we can assume that 1≤k σ k,r 1 ⊂ {λ ∈ C : Re λ < −M}.
First, we prove that the resolvent operator R(λ, L i,r 1 ) has a simple pole at λ j (r 1 ) ( j = 1, 2). We limit ourselves to proving this property for the eigenvalue λ 1 (r 1 ), since for the other one the proof is completely similar.
As it is immediately seen the coefficient in front of λ k does not vanish at λ = λ 1 (r 1 ). Hence, there exist a neighborhood U of λ 1 (r 1 ), k 0 ∈ N and a positive constant χ such that |D k,r 1 (λ)| ≥ χλ k for any k ≥ k 0 and any λ in U. From (3.5)-(3.7) we thus deduce that for any k ≥ k 0 , j = 1, 2, 3 and λ ∈ U.
Since D k,r 1 (λ 1 (r 1 )) 0 for any k 0, the previous estimate can be extended to any k ≥ 1. Hence, where the second term in the previous splitting defines a function with values in L((L 2 C ) 3 ) which is bounded in U.
To conclude that it is, actually, a simple eigenvalue, we have to show that the eigenspace associated with λ 1 (r 1 ) is one dimensional. This property follows from recalling that D k,r 1 (λ 1 (r 1 )) 0 if k ≥ 1. Hence, any eigenfunction associated with λ 1 (r 1 ) is a constant.
We can now prove the following theorem: max depends on α, γ, µ I , µ T and µ V (see the proof). Then, the following properties are satisfied.
(i) If N < N * (where N * is the first positive zero of the function H in (4.13)), then the periodic solution X 1 # (ε) is orbitally asymptotically stable with asymptotic phase. (ii) For any N > 0, the periodic solution X 2 # (ε) is orbitally asymptotically stable with asymptotic phase.
Proof. The arguments in Henry's book [11] (see also [3] in a more general situation) show that the stability of the bifurcated periodic solutions can be read on a Center Manifold. This allows to reduce our problem, which is set in a infinite dimensional Banach space, to a problem in a finite dimensional space.
To obtain this finite dimensional problem, we first need to determine the spectral projection associated to the eigenvalues −ω j i and ω j i ( j = 1, 2). As a general fact, such a projection is the sum of the spectral projections P j,+ , associated to the eigenvalue iω j , and P j,− , associated to the eigenvalue −iω j . Since iω j and −iω j are simple eigenvalues (see Theorem 4.1), there exists a unique projection on the eigenspace relative to iω j which commutes with L i . Similarly, there exists a unique projection of the eigenspace relative to −iω j which commutes with L i . Using these facts it is easy to check that In what follows we set ϕ j = (ϕ j,1 , ϕ j,2 , ϕ j,3 ) and ψ j = (ψ j,1 , ψ j,2 , ψ j,3 ).
As it is well known, P j allows to split (L 2 C ) 3 into the direct sum of the two subspaces In particular, P j maps (L 2 ) 3 into itself and allows us to split the space (L 2 ) 3 into the direct sum of the two subspaces P j ((L 2 ) 3 ) = {zϕ + z ϕ : z ∈ C} and (I − P j )((L 2 ) 3 ) = {u ∈ (L 2 ) 3 : (u, ψ) 2 = 0}.
Since an explicit computation of these coefficients for any value of T max is uneasy, and we are interested in large (enough) values of T max , as in [6,Sec. 4.3] we determine the sign of Re c 1 (r j ) via an asymptotic analysis as T max → +∞. We get where D(N) and H(N) are respectively given by Whereas D(N) is always positive, the sign of H(N) depends on N. Since H(0) < 0 and lim N→+∞ H(N) = +∞, the function H has at least a positive zero. We define by N * the (first) positive zero of H. Therefore, H(N) < 0 for 0 ≤ N < N * . We thus conclude that, for T max large enough (let us say T max > T (3) max > T (2) max , which depends on α, γ, µ I , µ T and µ V ), Re c 1 (r 1 ) < 0 for any 0 < N < N * , whereas Re c 1 (r 2 ) < 0 for any N > 0. This completes the proof.

Numerical results
In order to show the stability of the infected steady state numerically, we can fix N = 300, start from the value r crit,0 = 0.05625, increase the logistical parameter r monotonically until a critical condition is reached such that any further change would result in instability, other parameters can be found in Table 1. We present the graphs of numerical solution of the system (1.1)-(1.3) and the trajectory of the solution in the three-dimensional T -V-I space. Initial data are T 0 = T u + ε(sin x cos y), I 0 = 0.0, V 0 = 0.0185. Some figures assure that this solution approaches the limit cycle in the instability subdomain P. In Figure 7 corresponding to the subdomain P, the solution approaches the periodic orbit.    Let A be the realization of the Laplace operator ∆ in L 2 C , with H 2 ♯,C as a domain. The following is a well-known result. Nevertheless, for the reader's convenience we provide a short proof. Denote by e h the function defined by for any h ∈ Z. Then, the functions (x, y) → e h (x)e k (y) are an orthogonal basis of L 2 C . Hence, any function g ∈ L 2 C can be expanded into a Fourier series as follows: g(x, y) = u(x, y)e −k 1 (x)e −k 2 (y)dxdy e k 1 (x)e k 2 (y) =: for almost every (x, y) ∈ Ω ℓ . Multiplying both sides of (A.2) by e k 1 (x)e k 2 (y) and integrating over Ω ℓ , it thus follows that, if u ∈ H 2 ♯,C , then the Fourier coefficients of u solves the infinitely many equations λℓ 2 u k 1 ,k 2 + 4π 2 (k 2 1 + k 2 2 )u k 1 ,k 2 = ℓ 2 f k 1 ,k 2 .
The following classical result on Sturm-Liouville problems is the key tool to prove Theorem 2.1(iv).
Corollary A.2. Let d and µ be, respectively, a positive constant and a bounded measurable function. Further, let B : H 2 ♯ → L 2 be the operator defined by Bu = d∆u − µu for any u ∈ H 2 ♯ . Then, the spectrum of B 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 Ω ℓ .