COEXISTENCE AND COMPETITIVE EXCLUSION IN AN SIS MODEL WITH STANDARD INCIDENCE AND DIFFUSION

. This paper investigates a two strain SIS model with diﬀusion, spatially heterogeneous coeﬃcients of the reaction part and distinct diﬀusion rates of the separate epidemiological classes. First, it is shown that the model has bounded classical solutions. Next, it is established that the model with spa- tially homogeneous coeﬃcients leads to competitive exclusion and no coexistence is possible in this case. Furthermore, it is proved that if the invasion number of strain j is larger than one, then the equilibrium of strain i is unstable; if, on the other hand, the invasion number of strain j is smaller than one, then the equilibrium of strain i is neutrally stable. In the case when all diﬀusion rates are equal, global results on competitive exclusion and coexis- tence of the strains are established. Finally, evolution of dispersal scenario is considered and it is shown that the equilibrium of the strain with the larger diﬀusion rate is unstable. Simulations suggest that in this case the equilibrium of the strain with the smaller diﬀusion rate is stable.


1.
Introduction. Trade-off mechanisms for coexistence of pathogen strains have been a long standing question in the evolution of pathogens. In the basic SIR model with multiple strains, the main outcome is competitive exclusion [3] but multiple mechanisms have been found to cause strain coexistence [16]. One of the most natural trade-off mechanisms is coexistence generated by a spatial segregation. It is intuitively clear that two or more strains can coexists if they occupy different spatial niches, a process called in ecology niche partitioning.
Competition and dispersal have been more thoroughly studied mathematically with Lotka-Volterra competition models or chemostat models with diffusion, rather than epidemic models. In the context of chemostat models, coexistence in the presence of diffusion was demonstrated by Waltman [20]. Competition in a Lotka-Volterra model with diffusion was studied extensively (see a review of models in [17]).
To demonstrate mathematically that multiple strains can coexist when spatial component is involved, one may consider either a multi-patch or a diffusion model. Article [18] considered a multi-patch multi-strain model and found coexistence caused by the spatial component. Article [10] considered an avian influenza model with multiple strains, diffusion and non-local transmission terms. They found coexistence, but it was most likely generated by a different trade-off mechanism, namely mutation. Tuncer and Martcheva [19] considered a two-strain version of a SIS model with diffusion, first introduced by Allen et al [2]. What is distinct about the model in [2] is that in contrast to most diffusion models studied which have constant coefficients (e.g. [10]), the model in [2] has spatially heterogeneous coefficients. Article [2] defined R 0 in this case, and showed that besides the disease-free equilibrium, the model has a unique endemic equilibrium. Tuncer and Martcheva [19] introduced the invasion numbers of the two-strain version and showed analytically and numerically that if both invasion numbers are larger than one, then there is a coexistence equilibrium.
Here we consider the same SIS epidemic model with diffusion and two strains, considered in [19]. In general the model has spatially heterogeneous coefficients in the reaction part and distinct constant diffusion coefficients. Similar to [2] we define the reproduction numbers of the strains and the invasion numbers of the strains [19]. Although the model was also studied in [19], here we concentrate on some aspects that were not discussed in [19]. In particular, we want to study analytically under what conditions the model leads to competitive exclusion, and what characteristics of the model really imply coexistence. We begin by showing that a classical solution exists globally and it is bounded. We show rigorously that if the coefficients of the reaction part are constant, then coexistence does not occur. This rules out the possibility that diffusion by itself may be a factor and shows that the assumption of spatially heterogeneous coefficients is necessary for coexistence. We revisit the proofs of the fact that if the strain j invasion number is larger than one, then strain i equilibrium is unstable. We further show that if the invasion number of strain j is smaller than one, then the equilibrium of strain i is linearly neutrally stable. We state a result, obtained in [19], that if both invasion numbers are grater than one, then there is a coexistence equilibrium. In the case when all three diffusion rates are equal, we establish some global results. In particular, we show that the boundary equilibria are globally stable under certain conditions, and that the system is permanent, if both reproduction numbers and both invasion numbers are larger than one. Finally we discuss the evolution of dispersal. Evolution of dispersal has been a topic of interest in the ecological modeling literature (see e.g. [5,13,14,15]). Prior results from ecological models suggest that competitors evolve towards a smaller diffusion rate [6,9]. Here we show that if the transmission and recovery rates are spatially heterogeneous and identical for both strains while diffusion rates for both strains are distinct then there is no coexistence equilibrium. In addition the equilibrium of the competitor with the larger diffusion rate is unstable. Numerical simulations suggest that the equilibrium of the competitor with smaller diffusion rate is stable but we were unable to show that rigorously.
In the next section, we introduce the model and the reproduction numbers and invasion numbers. In section 3 we state the main results. In section 4 we provide some numerical simulations to complement the analytical results. In section 5 we prove the main results.
2. The model. Suppose that a species is living in a bounded domain Ω ⊆ R m with smooth boundary ∂Ω, and it has the possibility to get infected by a disease with two pathogen strains. Let S(x, t) be the density of the susceptible individuals and I i (x, t) be the density of the infected individuals with strain i (i = 1, 2) at position x and time t. The movement of susceptible and infected populations is modeled by random walk described by the Laplacian operator. The positive constants d s , d 1 and d 2 denote the corresponding diffusion rates for the susceptible and infected populations. Environmental effects are incorporated in the model via spatially varying transmission (β 1 (x), β 2 (x)) and recovery (γ 1 (x), γ 2 (x)) rates. The infection rate of susceptible individuals is described by the standard incidence. The diffusive two-strain SIS model [19] is formulated as x ∈ Ω, t > 0, with no-flux boundary condition Neumann boundary conditions mean that there is no population flux though the boundary ∂Ω, and that susceptible and infected individuals remain in the domain Ω. Here the Hölder continuous functions β i (x) and γ i (x) (i = 1, 2) are the spatially dependent disease transmission rates and recovery rates, respectively. We assume that the initial data S(x, 0) and I i (x, 0), i = 1, 2, are continuous, and the number of infected individuals is positive initially: Let where N is the total number of individuals at t = 0. Adding the three equations in (1) and integrating over Ω, we get ∂ ∂t Ω S(x, t) + I 1 (x, t) + I 2 (x, t)dx = 0, which indicates that the total population remains a constant for all time t > 0.
When the movement of infected individuals are affected by the disease we set d s = d 1 = d 2 . If the infected individuals continue to move despite the disease then we set d s = d 1 = d 2 .
Next, we will briefly summarize the results regarding to the dynamics of the model (1)-(4), but for more details one can refer to [19]. To study the dynamics of the model, we first obtain the equilibrium solutions. The model (1)-(4) has four feasible equilibria. The unique disease free equilibrium (DFE), E 0 := (S 0 , I 0 1 , I 0 2 ) = (N/|Ω|, 0, 0), where |Ω| denotes the measure of the domain Ω, and E 0 represents the state in which both strains die out in the population. We consider two types of endemic equilibria (EE): strain one/two dominance equilibria and coexistence equilibria. The I 1 -only equilibrium E 1 = (S * , I * 1 , 0) represents the state with I * 1 (x) > 0 for some x ∈ Ω. Similar statement is true for I 2 -only equilibrium, E 2 = (Ŝ, 0,Î 2 ) withÎ 2 (x) > 0 for some x ∈ Ω. The coexistence endemic equilibrium E 3 = (S,Ĩ 1 ,Ĩ 2 ) if it exists (it may not be unique) represents the state in whichĨ 1 andĨ 2 are both nonnegative and nontrivial.
The reproduction number of each strain can be characterized as a variational problem as stated in [19]. For i = 1, 2, the reproduction number of strain i is This definition is related to the principal eigenvalue λ * i of the following eigenvalue problem d i ∆ψ + (β i − γ i ) ψ + λ i ψ = 0 in Ω, and ∂ψ ∂n = 0 on ∂Ω.
It is known that R i > 1 if and only if λ * i < 0 and R i < 1 if and only if λ * i > 0. Whether a strain persists when circulating alone is related to the strain specific reproduction numbers R 1 and R 2 as stated in the following Lemma 2.1.
Lemma 2.1. The following statements hold: The basic reproduction number is defined as R 0 = max{R 1 , R 2 }. It can be proved that the DFE is globally asymptotically stable if R 0 < 1 and it is unstable if R 0 > 1. Moreover, a unique I i -equilibrium E i exists if and only if R i > 1 for i = 1, 2.
If the I 2 -equilibrium E 2 exists, i.e. R 2 > 1, the ability of strain one to invade strain two is described by the invasion numberR 1 , which is defined bŷ : ϕ ∈ H 1 (Ω) and ϕ = 0 , whereN =Ŝ +Î 2 . As before,R 1 is greater than one if and only if the principal eigenvalue of the following problem in Ω, and ∂ψ ∂n = 0 on ∂Ω is negative. Similarly, if the I 1 -equilibrium exists, one can define a corresponding invasion number aŝ where N * = S * + I * 1 , whose magnitude is determined by a related principal eigenvalue.
As it has been shown in [19], the question of the existence of a coexistence equilibrium can be reduced to study the positive solution (Ĩ 1 ,Ĩ 2 ) satisfying d 1Ĩ1 + d 2Ĩ2 ≤ 1 of the following elliptic problem with boundary condition The main objective of the present work is to investigate the impact of the coefficients and diffusion rates on the existence and stability of the equilibria. We show that there is no coexistence equilibrium if the coefficients are homogenous (while coexistence equilibrium may exist when the coefficients are heterogeneous by [19]), which demonstrates that spatial heterogeneity promotes coexistence. We study the stability of the boundary equilibria, and show that competitive exclusion and coexistence are both possible when the diffusion rates are the same for susceptible and infected individuals. We also give evidence that smaller diffusion rate is an advantage to the pathogens. In the following section we state the main results of this work.
3. Main results. The following result states that the solution of the model is bounded and thus exists globally. This actually implies that (1)-(4) has a global attractor.
where K is some positive constant independent of the initial condition.
We first assume that the coefficients are homogeneous, i.e β i and γ i are positive constants for i = 1, 2. We show that there is no coexistence equilibrium in this case, which indicates that we are likely to expect either extinction or competitive exclusion for the two strains (and that also suggests the necessity to consider spatial inhomogeneous coefficients). Notice that and the basic reproduction number is Basically, we show that when the transmission and recovery rates are assumed to be constants, the dynamics of the model (1)-(4) resembles the corresponding ODE model. For the dynamics of the corresponding ODE model, see [19].
The previous theorem states that if all the coefficients are homogeneous, then there is no coexistence steady state. Is this still the case when the coefficients β i and γ i , i = 1, 2, are spatially-dependent? We will see that this is not true. We begin our further study with the stability of the boundary equilibria.
One expects that the I 1 -equilibrium is stable ifR 2 < 1, but we are unable to prove it here. Instead, we show that: Theorem 3.4. Assume thatR 2 < 1 and d S = d 1 . Then the I 1 -equilibrium (S * , I * 1 , 0) is linearly neutrally stable or stable (we assume R 1 > 1 such that the I 1 -equilibrium exists).
If the two boundary equilibria exist and are unstable, then problem (1)-(4) has a coexistence equilibrium (this result has been proven in [19]).
We then study the global attractivity of the equilibria in the special case if all the diffusion rates are equal. The global dynamics of (1)-(4) can be largely determined by the magnitudes of the reproduction and invasion numbers. We present the following result (consult [1] for the proof) on the global dynamics of the model (The readers may need to see Section 5.4 for the definition ofǏ 1 andǏ 2 ).
Then the following statements hold.
exists at least one coexistence EE of (1)-(4). Moreover, Problem (1)-(4) is permanent in the sense that there exist positive constants m and M such that for all x ∈ Ω and t ≥ T , where T is dependent on the initial conditions. According to Theorem 3.6, if the coefficients β i and γ i are heterogeneous, then problem (1)-(4) may have a rich dynamical behavior. It is possible for one strain to exclude the other one, and it is also possible for the two strains to coexist. It is then natural to ask whether all cases in the theorem are possible. We will see that they are indeed possible as shown in the following result (For h ∈ C(Ω), we definē h = Ω h(x)dx/|Ω| and h + (x) = max{h(x), 0} for any x ∈ Ω).
Then the following statements hold.
then there exists a positive number d * such that both boundary equilibria E 1 and E 2 exist while E 1 is globally attractive for any d > d * ; Ifβ 2 /γ 2 >β 1 /γ 1 > 1 and 1 +γ 2 /β 2 < 2γ 1 /β 1 , then there exists a positive number d * such that both boundary equilibria E 1 and E 2 exist while E 2 is globally attractive for any Then there exists a positive number d * such that the model (1)-(4) has a coexistence equilibrium and it is permanent for any d < d * .
If β i ≥ γ i for i = 1, 2, then the condition in (c) can be simplified as which can be interpreted ecologically as niche partitioning. Finally, we consider the role of the diffusion rates in the model. In the classic diffusive logistic competition model, it is well-known that the species with smaller diffusion rate has competitive advantage over the one with larger diffusion rate. So a natural conjecture is that this is still true for the diffusive epidemic models. To investigate that, we assume that the coefficients for the two strains are the same, i.e. β 1 (x) = β 2 (x) = β(x) and γ 1 (x) = γ 2 (x) = γ(x), and strain one is different from strain two only by the diffusion rate, i.e. d 1 = d 2 . Our result does support the conjecture that the strain with smaller diffusion rate has competitive advantage.
For the next experiment, we set  shows that the solution converges to a coexistence equilibrium whose spatial structure is shown in Figure 5. Since strain two has a smaller diffusion rate, we expect that it has competitive advantage over strain one. Indeed, Figure 6 shows that the total population of strain one is decreasing to zero while strain two persists in the long run. This observation is consistent with what we expected, although it takes longer time to see the trend compared with the previous cases.

5.
Proof of the main results.

Proof of Theorem 3.1.
It is easy to see that the solution is positive in Ω × (0, T max ) by the maximum principle, where T max is the maximum existence time of solution. Adding the three equations in (1) and integrating over Ω, we find which together with the positivity of solution implies that  We first consider the case that β 1 /γ 1 > β 2 /γ 2 > 1. Let X = C(Ω) × C(Ω) be the ordered Banach space with positive cone K = C(Ω) + × (−C(Ω) + ). Let ≤ K be the order generated by K. Denote by X + = C(Ω) + × C(Ω) + the usual positive cone, which generates the usual order ≤ on X.
Let T : [0, ∞) × X + → X + be the semiflow induced by the solution of the following parabolic system subject to Neumann boundary condition The system (9)-(10) has only three boundary equilibria: F 0 = (0, 0), To see these three equilibria are the only boundary equilibria, we letĨ 2 = 0. Then the first equation of (9) becomes

It has been shown in Lemma 3.3 of [2] that this equation with Neumann boundary
condition has a unique positive equilibrium, which implies the uniqueness of I 1equilibrium. Moreover, one can easily check that F 1 is an I 1 -equilibrium, and hence it is the unique one. Similarly, we can see that F 2 is the unique I 2 -equilibrium. Let and We can check that We claim that D is invariant for T (t). To see this, let (Ĩ 1,0 ,Ĩ 2,0 ) ∈ D and we want to show T (t)(Ĩ 1,0 ,Ĩ 2,0 ) = (Ĩ 1 (·, t),Ĩ 2 (·, t)) ∈ D for all t > 0. For simplicity, we only consider the caseĨ 1,0 = 0 andĨ 2,0 = 0. Then we haveĨ 1 (x, t) > 0 and I 2 (x, t) > 0 for all x ∈ Ω and t > 0. Assume to the contrary that T (t) first leaves D at t 0 . By (11), we have that for i = 1, 2. LetJ i be the solution of the problem Then by the comparison principle, we haveĨ i <J i on Ω × (0, t 0 ] for i = 1, 2. Moreover, one can check that 1/d i is an upper solution of problem (12), for i = 1, 2, which contradicts the assumption that T (t)(Ĩ 1,0 ,Ĩ 2,0 ) leaves D at t 0 . By (11), T (t) : [0, ∞) × D → D is a strictly monotone semiflow with respect to the order ≤ K . Moreover, T (t) satisfies all the conditions of Proposition 2.4 in [8](also see the remark after this proposition), and thus exactly one of the following holds: (a) There exists a positive equilibrium of T in D; (b) There is a monotone orbit in D from F 1 to F 2 ; (c) There is a monotone orbit in D from F 2 to F 1 .
We claim that (a) is not possible. To see this, we consider the ordinary differential equation system The solution of this system induces a semiflow T : [0, ∞) × D → D , which is strictly monotone with respect to the order ≤ K (Here K is the positive cone . So similar to the semiflow T , T also satisfies exactly one of (a)-(c). For T , it is easy to check that there is no positive equilibrium in D . Hence there is either an orbit of T from F 1 to F 2 or from F 2 to F 1 , which can be viewed as an orbit for T if we take a real number as a constant function on Ω. Hence (a) is not possible for T . Therefore, the model (1)-(4) has no coexistence equilibrium. We then consider β 1 /γ 1 > 1 > β 2 /γ 2 or 1 ≥ β 1 /γ 1 > β 2 /γ 2 . In this case, we have R 2 < 1, which leads to I 2 → 0 as t → ∞ by Lemma 2.1. Hence there is no coexistence equilibrium.

Proof of Theorem 3.3 and 3.4.
We need the following well-known eigenvalue comparison result.
Lemma 5.1. Suppose that h ∈ L ∞ (Ω) and d > 0. Let λ 1 (d, h) be the principal eigenvalue of If with equality holds if and only if h 1 = h 2 a.e. in Ω. If h is non-constant, then We now prove that the I 1 -equilibrium is unstable ifR 2 > 1. Note thatR 2 > 1 if and only if λ * < 0, where λ * is the principal eigenvalue of the following problem We now prove Theorem 3.3.
We now prove Theorem 3.4.
Proof. We proceed analogously as the proof of Theorem 3.3, and this leads to the following eigenvalue problem By the principle of linearized stability, we need to show that all the eigenvalues of problem (20) have nonnegative real parts. Suppose that (λ, φ, ψ 1 , ψ 2 ) is an eigenpair. If ψ 2 = 0, then λ is an eigenvalue of problem (15). Hence, λ is real and λ ≥ λ * . SinceR 2 < 1, we have λ * > 0 and so λ > 0. Therefore, we suppose that ψ 2 = 0 and then problem (20) can be simplified as Adding the first two equations in (21) and by d S = d 1 , we have If φ + ψ 1 = 0, λ is an eigenvalue of −d S ∆, and it is well-known that it must be non-negative. Therefore, we suppose φ + ψ 1 = 0. Multiplying both sides of the second equation of (21) by ψ 1 and integrating over Ω, we get Therefore, λ is positive and the proof is completed.

5.4.
Equal diffusion rates case. In this section, we suppose that the susceptible and infected individuals have the same diffusion rates, i.e. d S = d 1 = d 2 = d.
depends continuously on h ∈ L ∞ (Ω) and d, and it satisfies that whereh is the spatial average of h, i.e.h = Ω h(x)dx/|Ω|. We now prove Theorem 3.7.
Proof. (a) Suppose thatβ 1 /γ 1 > 1 >β 2 /γ 2 . To show E 1 is globally attractive, it suffices to check R 1 > 1 > R 2 by Theorem 3.6. Note that the magnitudes of the basic reproduction numbers are determined by the signs of the principal eigenvalues of the corresponding eigenvalue problems. Since λ 1 (d, The second part of the statement can be proved analogously. (b) Supposeβ 1 /γ 1 >β 2 /γ 2 > 1 and 1 +γ 1 /β 1 < 2γ 2 /β 2 . Similar to (a), there exists d > 0 such that R 1 and R 2 are both greater than one, and so both E 1 and E 2 exist for all d > d . SinceÎ 2 satisfies we have thatÎ We claim that To see this, let > 0 be given. For the sake of convenience, denote λ 1 (d, for any d > d . By Lemma 5.1, we have for any d > d . It then follows from (24) that Then the claim follows since is arbitrary.
We know thatR 1 > 1 if and only if λ 1 (d, andǏ 1 depends onÎ 2 continuously, we havě Then similar to the previous claim, we can show that Hence, there exists d * >d such thatR 1 > 1 and λ 1 (d, β 2 − γ 2 − β 2Ǐ1 |Ω|/N ) ≥ 0 for any d > d * . So by Theorem 3.6, the proof of the first part is complete. The second part of the statement can be proved similarly.
(c) Since λ 1 (d, β i − γ i ) → min{−β i (x) + γ i (x) : x ∈ Ω} < 0, as d → 0, there exists d > 0 such that R i is greater than one for i = 1, 2 and for any d < d . So the two boundary equilibria exist for all d < d and , as d → 0.