Global stability of a multistrain SIS model with superinfection and patch structure

We study the global stability of a multistrain SIS model with superinfection and patch structure. We establish an iterative procedure to obtain a sequence of threshold parameters. By a repeated application of a result by Takeuchi et al. [Nonlinear Anal Real World Appl. 2006 7:235-247], we show that these parameters completely determine the global dynamics of the system: for any number of patches and strains with different infectivities, any subset of the strains can stably coexist depending on the particular choice of the parameters. Finally, we return to the special case of one patch examined in [Math Biosci Eng. 2017 14:421-35] and give a correction to the proof of Theorem 2.2 of that paper.


Introduction
Several viruses have different genetic variants (subtypes) called strains which may differ in their infectivity and virulence. Stronger strains might superinfect an individual already infected by another strain and there can be a coexistence of different virus strains with different virulence. Nowak [1] considered a model to provide an analytical understanding of the complexities introduced by superinfection. In our earlier work [2], we considered a multistrain SIS model with superinfection with n infectious strains and showed that it is possible to obtain a stable coexistence of any subgroup of the n strains. We established an iterative method for calculating a sequence of reproduction numbers, which determine the strains being present in the globally asymptotically stable coexistence equilibrium.
Marvá et al. [11] considered a spatially distributed periodic multistrain SIS epidemic model with patches of periodic migration rates without superinfection. Considering global reproduction numbers in the non-spatialized aggregated system that serve to decide the eradication or endemicity of the epidemic in the initial spatially distributed nonautonomous model, and comparing these global reproductive numbers with those corresponding to isolated patches, they showed that adequate periodic fast migrations can in many cases reverse local endemicity and get global eradication of the epidemic.
Motivated by our earlier work on multistrain models and by the recent results on spatial spread of diseases, we extend our previous model [2] to the general case of p patches. In Section 2, we establish a multistrain SIS model with superinfection with n infectious strains and patch structure. In Section 3, we establish an iterative procedure to determine the globally asymptotically stable equilibrium of the multipatch model introduced in Section 2. In Section 4, we turn to the case p = 1, studied in Dénes, Muroya and Röst [2] and give a correction of the proof of Theorem 2.2 of that paper.

The model
We consider a heterogeneous virus population with n virus strains having different infectivities and virulences. We will assume that superinfection is possible, and more virulent strains outcompete the less virulent ones in an infected individual taking over the host completely, i.e. we assume that an infected individual is always infected by only one virus strain. Let n denote the number of strains with different virulences while p stands for the number of patches. On each patch, the population is divided into n + 1 compartments depending on the presence of any of the virus strains: the susceptible class of patch ℓ is denoted by S ℓ (t) and on each patch ℓ, there are n infected compartments T ℓ 1 , . . . , T ℓ n where a larger index corresponds to a compartment of individuals infected by a strain with larger virulence, so for i < j, T j individuals superinfect T i individuals. Let B ℓ denote the birth rate and b ℓ the death rate on the ℓth patch. We denote by β ℓ kj the transmission rate on patch ℓ by which the kth strain infects those who are infected by the jth strain. The transmission rates from susceptibles to strain k on patch ℓ will be denoted by β ℓ kk . Recovery rate on patch ℓ among those infected by the kth strain will be denoted by θ ℓ k . By m ℓi we denote the travel rate from patch i to ℓ, which, on a given patch is equal for all compartments on that patch. Using these notations, we consider the following multistrain SIS model with superinfection and patch structure: where δ kj denotes the Kronecker delta such that δ kj = 1 if k = j and δ kj = 0 otherwise, and where β ℓ kj = β ℓ kk , 1 ≤ j ≤ k, and β ℓ kj = −β ℓ jj , k + 1 ≤ j ≤ n, k = 1, 2, . . . , n, ℓ = 1, 2, . . . , p.
Note that for n = 2 and p = 1, (2.1) corresponds to the model by A. Dénes and G. Röst describing the spread of ectoparasites and ectoparasite-borne diseases [12,13], while for p = 1, it corresponds to the multistrain SIS model by A. Dénes, Y. Muroya and G. Röst [2].

Main result
Let us introduce the notation Then, by (2.3), we have β ℓ kj = −β ℓ jk for k = j and hence, Thus, (2.1) is equivalent to (3.2c) The equations (3.2b)-(3.2c) are clearly independent from the rest of the equations. In particular, the equations (3.2c) are also independent from the equations (3.2b). As the coefficient matrix A of the linear system of equations is a strictly diagonally dominant Z-matrix, it is nonsingular and its inverse is positive, hence, this algebraic system has a unique, positive solution Let us define P ℓ (t) := N ℓ (t) − N ℓ * , ℓ = 1, . . . , p, then for P ′ ℓ (t), we have the equation From the properties of the matrix −A, applying the Gershgorin circle theorem, we obtain that P ℓ (t) → 0 exponentially as t → ∞, ℓ = 1, . . . , p. Hence, for the equations (3.2c), there exist positive constants N ℓ * n , ℓ = 1, 2, . . . , p such that exponentially and (3.2b) has the following limit system: which is a p-dimensional Lotka-Volterra system with patch structure, in the form as Equation (2.1) in Takeuchi et al. [14] We introduce the notationm (1 − δ iℓ )m iℓ , i = 1, 2, . . . , p, and define the connectivity matrix Let us denote by s(L) the stability modulus of a p × p matrix L, defined by s(L) := max{Re λ : λ is an eigenvalue of L}. If L has nonnegative off-diagonal elements and is irreducible, then s(L) is a simple eigenvalue of L with a (componentwise) positive eigenvector (see, e.g., Theorem A.5 in Smith [15]). Note that we may take that the populations go extinct in every patch not only if s(M n ) < 0 but also if s(M n ) = 0 (see Theorem 2.2 of Faria [16]).
Let E * n = (T 1 * n , T 2 * n , . . . , T p * n ) be the unique equilibrium of (3.5) which is globally asymptotically stable. Then, Therefore, in the first case, the unique equilibrium of (3.5), is globally asymptotically stable on {(T 1 n , T 2 n , . . . , T p n ) ∈ R p + }, while in the second case, the unique positive equilibrium E * n = (T 1 * n , T 2 * n , . . . , T p * n ) with T ℓ * n > 0, ℓ = 1, 2, . . . , p is globally asymptotically stable with respect to . This way, substituting T i * n , 1 = 1, . . . , p into the place of T i n (t) in (3.1) and (3.2), we may consider the following reduced system of (3.2) for the global stability of (2.1): It is easy to see that (3.6) is of similar structure as (3.2), but with dimension p(n−1)+1. The positivity of the new parameters follows from the conditions (2.3). This means that by repeating the above steps, namely, substituting the limit of the total populations in the patches and then substituting the limit of the Lotka-Volterra system for the strongest strain, we can further reduce the dimension by substituting the values of the equilibrium which is globally asymptotically stable, of the decoupled p dimensional Lotka-Volterra system into the remaining equations.
In general, after performing the above steps q times, we arrive at the system From the equations (3.7c), similarly as before, there exist positive constants N ℓ * n−q , ℓ = 1, 2, . . . , p such that lim and (3.7) has the following reduced limit system: Let us define Again, (3.9b) can be the decoupled from the rest of the equations as a p dimensional Lotka-Volterra system with patch structure: (3.10) Similarly as before, assuming the irreducibility of M n−q , this system has a globally attractive equilibrium (T 1 * n−q , T 2 * n−q , . . . , T p * n−q ), which is either the trivial equilibrium if s(M n−q ) ≤ 0 or a positive equilibrium if s(M n−q ) > 0.

and the new variables
We obtain the system which again, is a system with the same structure. In the end, we arrive at a p dimensional Lotka-Volterra system, the dynamics of which can be determined in a similar way as in the above case. This final system will give us an equilibrium value for S 1 (t) and (T 1 1 (t), T 1 2 (t), . . . , T 1 p (t)). Thus, by the above discussion, we can reach a conclusion by induction to the global dynamics of the model (2.1) and we formulate the following theorem. Proof of Theorem 3.2. The main part of the proof consists of the above description of the steps of the procedure. There is one point left to be shown: we have to prove that in each step, when we substitute the limits N ℓ * k , resp. T ℓ * k into the equations, the dynamics of the resulting system is indeed equivalent to that of the preceding one.
We summarize the steps of the procedure in the following.
For the validity of Step 3 in the qth cycle, we need to verify that M n−q is irreducible. Since M n−q = M + diag[c 1 n−q , . . . , c p n−q ] and we assumed that M is irreducible, M n−q is also irreducible. To obtain that in each case, the limit of the solutions of the resulting system after the substitution will be the same equilibrium as the limit of the solutions of the original system, we will apply Theorem 4.1 of Hirsch and Smith [17]. To apply this theorem, we recall the quasimonotone condition [17] for a differential equation x ′ (t) = f (t, x(t)): we say that the time-dependent vector field f : J × D → R n (where J ⊂ R and D ⊂ R n ) satisfies the quasimonotone condition in D if for all (t, y), (t, z) ∈ J × D, we have y ≤ z and According to Theorem 4.1 of Hirsch and Smith [17], if f, g : J × D → R n are continuous, Lipschitz on each compact subset of D, at least one of them satisfies the quasimonotone condition, and f (t, y) ≤ g(t, y) for all (t, y) ∈ J × D, then y, z ∈ R n , y ≤ z implies x(t; t 0 , y) ≤ x(t; t 0 , z) for all t > t 0 , where x(t; t 0 , y) denotes the solution of x ′ (t) = f (t, x(t)) started from y at t = t 0 .
To show that the limits T ℓ * k obtained during the procedure by substituting the limits of (3.10) into (3.9a) are the same as the limit of the variables T ℓ k , k = 1, . . . , n, ℓ = 1, . . . , p in the original system, we will use an induction argument. It is clear from the above that the claim is true for k = n. Let us now suppose that the claim is not true for all T ℓ k (t), then there exists a largest index 1 ≤ r ≤ n − 1 such that T m * r is not equal to the limit of T m r (t) in the original system for some 1 ≤ m ≤ p. The limits T ℓ * r are obtained by first substituting the limits T ℓ * r+1 into the equations for T ℓ j (t), 1 ≤ j ≤ r and then substituting the limits N ℓ * r into the equations for T ℓ r (t), hence, we have to compare the limits of the two systems (3.13) ℓ = 1, 2, . . . , p.
As we have assumed that for all larger indices, the limits of the compartments of the original system (3.2) are equal to the limits obtained during the procedure, using the equations for T 1 r (t), . . . , T p r (t) after n − r + 1 cycles of the procedure satisfy the quasimonotone condition and the comparison (3.14), the limits obtained for these have to coincide with those of the original system (for r = n, the statement follows directly).
To prove that not only attractivity, but also global asymptotic stability holds, we will again use induction. Let E = (S 1 ,T 1 1 , . . . ,T 1 n , . . . ,S p ,T p 1 , . . . ,T p n ) denote the equilibrium obtained at the end of the procedure, whereT j i = 0 orT j i > 0 depending on the stability moduli (s(M 1 ), s(M 2 ), . . . , s(M n )) and let E k = (S 1 ,T 1 1 , . . . ,T 1 k , . . . ,S p ,T p 1 , . . . ,T p k ) be the equilibrium of the p(k + 1)-dimensional system obtained during the procedure, consisting of the first p(k + 1) coordinates of E. Let us suppose that E k is a stable equilibrium of the p(k + 1)-dimensional reduced system for some k ≤ n. We will show that in each step, E k+1 is a stable equilibrium of the p(k + 2)-dimensional reduced system. Suppose this does not hold, i.e. E k+1 is unstable. In this case there exists an ε > 0 and is a sequence {x m } → E k+1 , |x m −E k+1 | < 1/m such that the orbits started from the points of the sequence leave B(E k+1 , ε) Let us denote by x ε m the first exit point from B(E k+1 , ε) of the solution started from x m , reached at time τ m . There is a convergent subsequence of the sequence x ε m (still denoted by x ε m ) which tends to a point denoted by We will show that the E k+1 ∈ α(x * ε ). For this end, let us consider the set S(E k+1 , ε 2 ). Clearly, all solutions started from the points x m (we drop the first elements of the sequence, if necessary) will leave the set B(E k+1 , ε 2 ). We denote the last exit point of each trajectory from this set before time τ m , respectively, by x ε/2 m . Also this sequence has a convergent subsequence (still denoted the same way), let us denote its limit by x * ε/2 . We will show that the trajectory started from x * ε/2 goes through x * ε . As E k+1 is globally attractive, this trajectory will eventually enter S(E k+1 , ε 4 ) at some time T > 0. Let us suppose that the trajectory started from x * ε/2 does not go through x * ε and let us denote by d > 0 the distance of this trajectory from x * ε . For continuity reasons, there is an N ∈ N so that for any m > N , |x * ε/2 t − x ε/2 m t| < max{ d 2 , ε 8 } for 0 < t < T . This means that for m large enough, the trajectory started from x ε/2 m will enter again S(E k+1 , ε 2 ) without getting close to x ε * which contradicts either x ε m being the first exit point from B(E k+1 , ε) or x ε/2 m being the last exit point before τ m from B(E k+1 , ε 2 ). Hence, we have shown that the trajectory started from x * ε/2 goes through x * ε . Proceeding like this (taking neighbourhoods of radius ε/4, ε/8 etc.) we obtain that the backward trajectory of x * ε enters any small neighbourhood of E k+1 as t → −∞, hence, E k+1 ∈ α(x * ε ), while it follows from the global attractivity of E k+1 that the ω-limit set of the trajectory is {E k+1 }. Let us denote this trajectory by γ(x * ε ) We know that the equations for T 1 k+1 (t), . . . , T p k+1 (t) and N 1 k+1 (t), . . . , N p k+1 (t) can be decoupled from the rest of the equations and using the exponential stability of the limits (3.4) and Proposition 3.1 we obtain thatT 1 k+1 , . . . ,T p k+1 is a stable equilibrium of the system consisting of the equations for d dt T 1 k+1 (t), . . . , d dt T p k+1 (t). Therefore, the equilibrium E k+1 is stable in the coordinates T 1 k+1 , . . . , T p k+1 in the sense that for anyε > 0 there exists aδ(ε) > 0 such that for any initial value x with |x − E k+1 | <δ, |T ℓ k+1 (t)−T ℓ k+1 | <ε for all t > 0 and ℓ = 1, . . . , p. Thus, the trajectory γ(x * ε ) obtained above lies entirely in the subspace {T 1 k+1 =T 1 k+1 , . . . , T p k+1 =T p k+1 }. On the other hand, the current p(k + 2)-dimensional system coincides with the p(k + 1)-dimensional system on this subspace. For the latter system, stability of the equilibrium E k follows from the induction assumption. However, the existence of an orbit whose ω-limit set is {E k+1 } and whose α-limit set contains E k+1 contradicts the stability of the equilibrium E k . This implies the global asymptotic stability of the equilibrium of the p(k + 2)-dimensional system.
For k = 1, the assertion holds trivially, hence, repeating the inductive step we obtain global asymptotic stability of the equilibrium E. In this section, we consider the special case of one patch examined in Dénes, Muroya and Röst [2] and give a correction to the proof of Theorem 2.2 of that paper. First, we recall this theorem about the globally asymptotically stable equilibrium of the multistrain SIS model (1 − δ kj )β kj T j (t) − (b + θ k + δ kn d n )T k (t), k = 1, 2, . . . , n, (4.1) with initial conditions S (0) = φ 0 , T k (0) = φ k , k = 1, 2, . . . , n, where δ kj denotes the Kronecker delta such that δ kj = 1 if k = j and δ kj = 0 otherwise, and Γ = [0, ∞) n+1 . We assume that the conditions hold for the infection rates for k = 1, 2, . . . , n, i.e. we assume that the k-th strain infects those who are infected by a milder strain (including the non-infected) with the same rate. The notation d n stands for disease-induced death rate for the most infectious strain. In our previous work [2], we gave an iterative procedure (similar to the one introduced in Section 3 of the present paper) to calculate a sequence of reproduction numbers which completely determines the global dynamics of the system. In the general step of the procedure we consider the system and dT n−ℓ (t) dt = S(t)β n−ℓ,n−ℓ T n−ℓ (t) + T n−ℓ (t) T k (t), and (4.4) Again, (4.4) might be decoupled from the other equations (4.2). For R (n−ℓ) 0 ≤ 1, system (4.4) has only the trivial equilibrium (0, 0). But for R (n−ℓ) 0 > 1, system (4.4) has two equilibria: the trivial equilibrium (0, 0) and the non-trivial equilibrium Then, from (4.2), we obtain the systems if R (n−ℓ−1) 0 ≤ 1, which, again, are systems with the same structure. In the end, we arrive at the two-dimensional system dS(t) dt = B (n−1) − b (n−1) S(t) − β 11 S(t)T 1 (t) + θ 1 T 1 (t), which has the two equilibria B (n−1) b (n−1) , 0 and b (n−1) + θ 1 β 11 , with the latter one only existing if R (n) 0 := B (n−1) β 11 b (n−1) (b (n−1) + θ 1 ) > 1.
The dynamics of this system can be determined in a similar way as in the case of (4.4), and we obtain that the first equilibrium is globally asymptotically stable if R  Proof. Let us suppose that there exists a solution started with positive initial values whose limit is not the equilibrium E obtained at the end of the procedure described in Dénes, Muroya, Röst [2]. It follows from the procedure that the last coordinate tends to the last coordinate of E. There exists a maximal index k (0 ≤ k ≤ n − 1) such that the kth coordinate of the solution does not tend to the kth coordinate of E, while all coordinates with index larger than k do tend to the corresponding coordinate of E. Let us consider the kth equation in the original system: Introducing the notation T 0 (t) := S(t), let us defineÑ k (t) as with respect to the original system. Hence, we can write the equation for T k (t) as