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

DOI: 10.1002/mma.6636

R E S E A R C H A R T I C L E

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

Attila Dénes

Yoshiaki Muroya

Gergely Röst

1Bolyai Institute, University of Szeged, Szeged, Hungary

2Department of Mathematics, Waseda University, Tokyo, Japan

3Mathematical Institute, University of Oxford, Oxford, United Kingdom Correspondence

Attila Dénes, Bolyai Institute, University of Szeged, H-6720 Szeged, Hungary.

Email: denesa@math.u-szeged.hu Communicated by: R. Bravo de la Parra Funding information

EU-funded Hungarian, Grant/Award Number: EFOP-3.6.1-16-2016-00008;

János Bolyai Research Scholarship of the Hungarian Academy of Sciences; Ministry for Innovation and Technology,

Grant/Award Number:

TUDFO/47138-1/2019-ITM; National Research, Development and Innovation Fund of Hungary, Grant/Award Number:

FK 124016 and PD 128363

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.

K E Y WO R D S

global asymptotic stability, multigroup epidemic model, multistrain model, patch model M S C C L A S S I F I C AT I O N

92D30; 34D23

1 I N T RO D U CT I O N

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¹ considered a model to provide an analytical understanding of the complexities introduced by superinfection. In our earlier work,²we considered a multistrain SIS model with super infection withninfectious strains and showed that it is possible to obtain a stable coexistence of any subgroup of then 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.

Recently, there has been an increasing interest in the modelling of the spatial spread of infectious diseases (see, e.g., Arino and Portet,³Knipl,⁴Knipl and Röst,⁵Muroya, Kuniya and Enatsu,⁶Nakata and Röst⁷). There are several ways to model spatial spread: one might use partial differential equations (see, e.g., Peng and Zhao,⁸Allen et al.,⁹Ge et al.¹⁰) or one may apply ordinary or functional differential equations where individuals can travel between different patches (countries, regions, cities etc.).

Marvá et al.¹¹considered a spatially distributed periodic multistrain SIS epidemic model with patches of periodic migra- tion rates without superinfection. Considering global reproduction numbers in the nonspatialized aggregated system that

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© The Authors. Mathematical Methods in the Applied Sciences published by John Wiley & Sons Ltd

Math Meth Appl Sci. 2020;1–10. wileyonlinelibrary.com/journal/mma 1

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²to the general case ofppatches. In Section 2, we establish a multistrain SIS model with superinfection withninfectious 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.

2 T H E M O D E L

We consider a heterogeneous virus population withnvirus 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 indi- vidual taking over the host completely, that is, we assume that an infected individual is always infected by only one virus strain. Letndenote the number of strains with different virulences, whereaspstands for the number of patches. On each patch, the population is divided inton+1 compartments depending on the presence of any of the virus strains: the sus- ceptible class of patch𝓁is denoted byS^𝓁(t)and on each patch𝓁, there areninfected compartmentsT₁^𝓁,…,T_n^𝓁where a larger index corresponds to a compartment of individuals infected by a strain with larger virulence, so fori<j,Tjindi- viduals superinfectT_iindividuals. LetB^𝓁denote the birth rate andb^𝓁the death rate on the𝓁th patch. We denote by𝛽_k𝑗^𝓁 the transmission rate on patch𝓁by which thekth strain infects those who are infected by thejth strain. The transmission rates from susceptibles to strainkon patch𝓁will be denoted by𝛽_kk^𝓁. Recovery rate on patch𝓁among those infected by thekth strain will be denoted by𝜃^𝓁_k. Bym_𝓁i, we denote the travel rate from patchito𝓁, which, on a given patch is equal for all compartments on that patch. This assumption—which is natural in the case of mild diseases—will be important in the transformation of variables described in the next section. The parametersB^𝓁,b^𝓁,m_𝓁i,𝓁,i= 1,…,pare assumed to be nonnegative.

Using these notations, we consider the following multistrain SIS model with superinfection and patch structure:

dS^𝓁(t)

dt =B^𝓁−b^𝓁S^𝓁(t) −S^𝓁(t)

∑n

k=1

𝛽_kk^𝓁T_k^𝓁(t) +

∑n

k=1

𝜃^𝓁_kT^𝓁_k(t) +

∑p

i=1

(1−𝛿_𝓁i){

m_𝓁iSⁱ(t) −mi𝓁S^𝓁(t)} , dT^𝓁

k(t)

dt =S^𝓁(t)𝛽_kk^𝓁T_k^𝓁(t) +T_k^𝓁(t)

∑n 𝑗=1

(1−𝛿k𝑗)𝛽_k𝑗^𝓁T^𝓁_𝑗(t) −(

b^𝓁+𝜃_k^𝓁) T^𝓁_k(t) +

∑p

i=1

(1−𝛿_𝓁i){

m_𝓁iT_kⁱ(t) −mi𝓁T_k^𝓁(t)} , k=1,2,…,n, 𝓁=1,2,…,p,

(1)

with initial conditions

S^𝓁(0) =𝜑^𝓁₀,T_k^𝓁(0) =𝜑^𝓁_k,k=1,2,…,n,𝓁=1,2,…,p, (𝜑¹₀, 𝜑¹₁, 𝜑¹₂,…, 𝜑¹n, 𝜑²₀, 𝜑²₁, 𝜑²₂,…, 𝜑²n,…, 𝜑^p₀, 𝜑^p₁, 𝜑^p₂,…, 𝜑^pn

)∈R^(n+1)p₊ ^{=∶ Γ}, (2)

where𝛿kjdenotes the Kronecker delta such that𝛿kj=1 ifk=jand𝛿kj=0 otherwise, and where 𝛽_k^𝓁_𝑗=𝛽_kk^𝓁, 1≤𝑗≤k, and

𝛽_k^𝓁_𝑗= −𝛽_𝑗𝑗^𝓁, k+1≤𝑗≤n, k=1,2,…,n, 𝓁=1,2,…,p. (3)

Note that forn=2 andp=1, (1) corresponds to the model by A. Dénes and G. Röst describing the spread of ectopar- asites and ectoparasite-borne diseases,^12,13 whereas forp= 1, it corresponds to the multistrain SIS model by A. Dénes, Y. Muroya and G. Röst.²

3 M A I N R E S U LT

Let us introduce the notation

N_n^𝓁(t) =S^𝓁(t) +

∑n 𝑗=1

T_𝑗^𝓁(t),𝓁=1,2,…,p. (4)

Then, by (3), we have𝛽_k^𝓁_𝑗= −𝛽_𝑗^𝓁_kfork≠jand hence,

∑n

k=1

T_k^𝓁(t)

∑n 𝑗=1

(1−𝛿k𝑗)𝛽_k^𝓁_𝑗T_𝑗^𝓁(t) =0,𝓁=1,2,…,p.

Thus, (1) is equivalent to dT^𝓁_k(t)

dt = (

N_n^𝓁(t) −

∑n 𝑗=1

T_𝑗^𝓁(t) )

𝛽_kk^𝓁T_k^𝓁(t) +T_k^𝓁(t)

∑n 𝑗=1

(1−𝛿k𝑗)𝛽_k^𝓁_𝑗T^𝓁_𝑗(t) −(

b^𝓁+𝜃_k^𝓁) T_k^𝓁(t)

+

∑p

i=1

(1−𝛿_𝓁i){

m_𝓁iT_kⁱ(t) −mi𝓁T_k^𝓁(t)}

, k=1,2,…,n−1,

(5a)

dT_n^𝓁(t) dt =

( N_n^𝓁(t) −

∑n 𝑗=1

T^𝓁_𝑗(t) )

𝛽nn^𝓁 T_n^𝓁(t) +T_n^𝓁(t)

∑n 𝑗=1

(1−𝛿n𝑗)𝛽_n𝑗^𝓁T^𝓁_𝑗(t)

−( b^𝓁+𝜃^𝓁n

)T_n^𝓁(t) +

∑p

i=1

(1−𝛿𝓁i){

m_𝓁iTⁱ_n(t) −m_i𝓁T^𝓁_n(t)}

=T^𝓁_n(t) (

𝛽nn^𝓁 N_n^𝓁(t) −

∑n 𝑗=1

{𝛽nn^𝓁 − (1−𝛿n𝑗)𝛽_n^𝓁_𝑗}

T_𝑗^𝓁(t) −( b^𝓁+𝜃n^𝓁

)) +

∑p

i=1

(1−𝛿𝓁i){

m_𝓁iT_nⁱ(t) −m_i𝓁T_n^𝓁(t)} ,

=T^𝓁_n(t)(

𝛽nn^𝓁 N_n^𝓁(t) −𝛽nn^𝓁T^𝓁_n(t) −( b^𝓁+𝜃n^𝓁

))+

∑p

i=1

(1−𝛿𝓁i){

m_𝓁iTⁱ_n(t) −m_i𝓁T^𝓁_n(t)} ,

(5b)

dN_n^𝓁(t)

dt =B^𝓁−b^𝓁N_n^𝓁(t) +

∑p

i=1

(1−𝛿𝓁i){

m_𝓁iN_nⁱ(t) −m_i𝓁N_n^𝓁(t)} , 𝓁=1,2,…,p.

(5c)

Equations (5b) and (5c) are clearly independent from the rest of the equations. In particular, Equation (5c) are also independent from Equation (5b). As the coefficient matrixAof the linear system of equations

(B¹

⋮ B^p

)

=

⎛⎜

⎜⎜

⎝

b¹+∑p

i=1(1−𝛿1i)m_i1 −m₁₂ … −m_1p

−m21 b²+∑p

i=1(1−𝛿2i)mi2 … −m2p

⋮ ⋮ ⋱ ⋮

−mp1 −mp2 … b^p+∑p

i=1(1−𝛿pi)mip

⎞⎟

⎟⎟

⎠ (N_n¹

⋮ N_n^p

)

is a strictly diagonally dominantZ-matrix, it is nonsingular, and its inverse is nonnegative (because of the nonnegativity of the parameters), hence, this algebraic system has a unique, positive solution

(N_n^1∗

⋮ N_n^p∗

)

=A⁻¹ (B¹

⋮ B^p

) .

Let us defineP_𝓁(t) ∶=N_n^𝓁(t) −N_n^𝓁∗, 𝓁=1,…,p, then forP_𝓁^′(t), we have the equation d

dt (P1(t)

⋮ P_p(t)

)

= −A (P1(t)

⋮ P_p(t)

)

. (6)

From the properties of the matrix−A, applying the Gershgorin circle theorem, we obtain thatP_𝓁(t)→0 exponentially ast→∞,𝓁=1,…,p. Hence, for Equation (5c), there exist positive constantsN_n^𝓁∗, 𝓁=1,2,…,psuch that

t→+∞lim N_n^𝓁(t) =N_n^𝓁∗, 𝓁=1,2,…,p, (7)

exponentially and (5b) has the following limit system:

dT^𝓁_n(t)

dt =T_n^𝓁(t)(

𝛽nn^𝓁 N_n^𝓁∗−( b^𝓁+𝜃^𝓁n

)−𝛽nn^𝓁 T^𝓁_n(t)) +

∑p

i=1

(1−𝛿𝓁i){

m_𝓁iT_nⁱ(t) −m_i𝓁T_n^𝓁(t)}

, 𝓁=1,2,…,p, (8)

which is ap-dimensional Lotka–Volterra system with patch structure, in the form as Equation (2.1) in Takeuchi et al.¹⁴ We introduce the notation

̃ m_ii=

∑p 𝓁=1

(1−𝛿i𝓁)m_i𝓁, i=1,2,…,p,

and define the connectivity matrix

M=

⎡⎢

⎢⎢

⎣

−m̃₁₁ m₁₂ … m_1p m21 −m̃22 … m2p

⋮ ⋮ ⋱ ⋮

mp1 mp2 … −m̃pp

⎤⎥

⎥⎥

⎦ .

Now, we define

c^𝓁_n=𝛽nn^𝓁 N_n^𝓁∗− (b^𝓁+𝜃n^𝓁), 𝓁=1,2,…,p,

and

Mn=

⎡⎢

⎢⎢

⎣

c¹_n−m̃11 m12 … m1p

m21 c²_n−m̃22 … m2p

⋮ ⋮ ⋱ ⋮

mp1 mp2 … c^p_n−m̃pp

⎤⎥

⎥⎥

⎦ .

Let us denote bys(L)the stability modulus of ap×pmatrixL, defined bys(L) ∶=max{Re𝜆∶𝜆is an eigenvalue ofL}. If Lhas nonnegative off-diagonal elements and is irreducible, thens(L)is a simple eigenvalue ofLwith a (componentwise) positive eigenvector (see, e.g., Theorem A.5 in Smith¹⁵).

Proposition 1 (see Theorem 2.1 in Takeuchi et al.¹⁴).Suppose that M_nis irreducible. Then, Equation (8) has a positive equilibrium which is globally asymptotically stable if s(Mn)>0. If s(Mn) ≤ 0, then 0 is a globally asymptotically stable equilibrium, and the populations go extinct in every patch.

Note that we may take that the populations go extinct in every patch not only ifs(Mn) < 0 but also ifs(Mn) = 0 (see Theorem 2.2 of Faria¹⁶).

Let E^∗_n = (T_n^1∗,T^2∗_n ,…,T_n^p∗) be the unique equilibrium of (8) which is globally asymptotically stable. Then,E_n^∗ = (0,0,…,0)ifs(M_n) ≤ 0, andE^∗_n = (T_n^1∗,T^2∗_n ,…,T_n^p∗)satisfiesT_n^𝓁∗ >0, 𝓁 =1,2,…,p, ifs(M_n) >0. Therefore, in the first case, the unique equilibrium of (8), is globally asymptotically stable on{(T_n¹,T_n²,…,T_n^p) ∈R^p+}, whereas in the sec- ond case, the unique positive equilibriumE_n^∗= (T_n^1∗,T_n^2∗,…,T^p∗_n )withT_n^𝓁∗>0, 𝓁=1,2,…,pis globally asymptotically stable with respect to{(T_n¹,T_n²,…,T_n^p) ∈R^p+}⧵{(0,0,…,0)}. Let us introduce the notations

N_n−1^𝓁 (t) =S^𝓁(t) +

n−1∑

𝑗=1

T_𝑗^𝓁(t), 𝓁=1,2,…,p,

and

b^𝓁₍₁₎ =b^𝓁−𝛽_kn^𝓁T_n^𝓁∗ =b^𝓁+𝛽nn^𝓁 T_n^𝓁∗, k=1,2,…,n−1, 𝓁=1,2,…,p,

and

B^𝓁₍₁₎ =B^𝓁+𝜃n^𝓁T_n^𝓁∗, 𝓁=1,2,…,p,

where(T_n^1∗,…,T_n^1∗)is either equal to(0,…,0)(ifs(Mn) ≤ 0) or it is equal to the unique positive equilibrium of (8) (if s(Mn)>0). This way, substitutingT_n^i∗, 1 =1,…,pinto the place ofTⁱ_n(t)in (4) and (5), we may consider the following reduced system of (5) for the global stability of (1):

dT^𝓁

k(t) dt =

(

N_n−1^𝓁 (t) −

∑n−1 𝑗=1

T_𝑗^𝓁(t) )

𝛽_kk^𝓁T_k^𝓁(t) +T_k^𝓁(t)

n−1∑

𝑗=1

(1−𝛿k𝑗)𝛽_k𝑗^𝓁T_𝑗^𝓁(t) − (

b^𝓁₍₁₎+𝜃_k^𝓁) T_k^𝓁(t)

+

∑p

i=1

(1−𝛿𝓁i){

m_𝓁iT_kⁱ(t) −m_i𝓁T_k^𝓁(t)}

,k=1,2,…,n−2,

(9a)

dT_n−1^𝓁 (t)

dt =

(

N_n−1^𝓁 (t) −

∑n−1 𝑗=1

T^𝓁_𝑗(t) )

𝛽_n−1,n−1^𝓁 T^𝓁_n−1(t) +T_n−1^𝓁 (t)

∑n−1 𝑗=1

(1−𝛿n−1,𝑗)𝛽_n−1,𝑗^𝓁 T_𝑗^𝓁(t) − (

b^𝓁₍₁₎+𝜃_n−1^𝓁 ) T_n−1^𝓁 (t)

+

∑p

i=1

(1−𝛿𝓁i){

m_𝓁iTⁱ_n−1(t) −m_i𝓁T_n−1^𝓁 (t)}

=T_n−1^𝓁 (t)

(𝛽_n−1,n−1^𝓁 N_n−1^𝓁 (t) −𝛽_n−1,n−1^𝓁 T_n−1^𝓁 (t) − (

b^𝓁₍₁₎+𝜃_n−1^𝓁 )) +

∑p

i=1

(1−𝛿_𝓁i){

m_𝓁iTⁱ_n−1(t) −mi𝓁T_n−1^𝓁 (t)} ,

(9b) dN_n−1^𝓁 (t)

dt =B^𝓁₍₁₎−b^𝓁₍₁₎N_n−1^𝓁 (t) +

∑p

i=1

(1−𝛿_𝓁i){

m_𝓁iN_n−1ⁱ (t) −mi𝓁N_n−1^𝓁 (t)} , 𝓁=1,2,…,p.

(9c)

It is easy to see that (9) is of similar structure as (5), but with dimensionpn. The positivity of the new parameters fol- lows from the conditions (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 decoupledpdimensional Lotka–Volterra system into the remaining equations.

We proceed repeating the same steps for the newly arising reduced system, decreasing the dimension of our system in each round of the procedure byp. In each round, forqdecreasing fromn−1 to 1, we introduce the respective limits N_q^𝓁∗andT_q^𝓁∗, as well as the matricesM_qcorresponding to the reduced system in an analogous way as it was presented in the case of the original system. In the end, we arrive at apdimensional 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 forS¹(t)and (T₁¹(t),T¹₂(t),…,T_p¹(t)). Thus, by the above discussion, we can reach a conclusion by induction to the global dynamics of the model (1) and we formulate the following theorem.

Theorem 1. Assume that the connectivity matrix M is irreducible. Then the global dynamics of the multistrain, mul- tipatch SIS model (1) is completely determined by the threshold parameters(s(M1),s(M2),…,s(Mn))which can be obtained iteratively. There exists an equilibrium inΓwhich is globally asymptotically stable with respect to the regionΓ0, whereΓ0is the interior ofΓ.

Proof of Theorem 1. 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 limitsN_𝜅^𝓁∗andT_𝜅^𝓁∗, respectively, into the remaining 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.

1. We obtainN_n^𝓁∗ (𝓁=1,…,p)from the linear system (6).

2. We substitute the limitsN_n^𝓁∗ (𝓁=1,…,p)into Equation (5b) to obtain Equation (8).

3. We obtain the limitsT_n^𝓁∗(𝓁=1,…,p)of the Lotka–Volterra system (8).

4. We create the new variablesN_n−1^𝓁 (t), 𝓁=1,…,pand parametersb^𝓁₍₁₎,B^𝓁₍₁₎,𝓁=1,…,p.

5. We substitute the limitsT_n^𝓁∗ (𝓁 = 1,…,p)into Equation (5a) to obtain the reduced system (9) which has the same structure as the original one (5).

6. We repeat this cyclen−1 times, with the indices decreased by 1 every time.

For the validity of Step 3 in the qth cycle, we need to verify that M_n−q is irreducible. Because M_n−q = M + diag[c¹_n−q,…,c^p_n−q]and we assumed thatMis irreducible,Mn−qis 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.¹⁷To apply this theorem, we recall the quasimonotone condition¹⁷for a differential equationx^′(t) =f(t,x(t)): we say that the time-dependent vector field𝑓 ∶J×D→R^n(whereJ ⊂R^andD⊂ Rⁿ) satisfies the quasimonotone condition inDif for all(t,y),(t,z) ∈J×D, we have

𝑦≤zand𝑦i=z_iimplies𝑓i(t, 𝑦)≤𝑓i(t,z).

According to Theorem 4.1 of Hirsch and Smith,¹⁷ if𝑓,g ∶J×D→ Rⁿare continuous, Lipschitz on each compact subset ofD, at least one of them satisfies the quasimonotone condition, andf(t,y) ≤ g(t,y)for all(t,y) ∈J×D, then

𝑦,z∈Rⁿ, 𝑦≤zimpliesx(t;t0, 𝑦)≤x(t;t0,z)for allt>t0, wherex(t;t0,y)denotes the solution ofx^′(t) =f(t,x(t))started fromyatt=t0.

To show that the limitsT_𝜅^𝓁∗ obtained during the procedure by substituting the limits of (8) into (5a) are the same as the limit of the variablesT_𝜅^𝓁,𝜅 =1,…,n,𝓁=1,…,pin the original system, we will use an induction argument.

It is clear from the above that the claim is true for𝜅 =n. Let now 1 ≤ r ≤ n−1 and let us suppose that the claim holds for allT_𝜅^𝓁(t)forr< 𝜅 ≤ n. The limitsT_r^𝓁∗ are obtained by first substituting the limitsT_r+1^𝓁∗ into the equations forT_𝑗^𝓁(t), 1 ≤ j ≤ rand then substituting the limitsN_r^𝓁∗into the equations forT_r^𝓁(t), hence, we have to compare the limits of the two systems

dT_r^𝓁(t) dt =(

N_r+1^𝓁 (t) −2T_r+1^𝓁 (t) −T_r^𝓁(t))

𝛽rr^𝓁T_r^𝓁(t) − (

b^𝓁_(n−r+1)+𝜃r^𝓁

) T^𝓁_r(t) +

∑p

i=1

(1−𝛿_𝓁i){

m_𝓁iT_rⁱ(t) −mi𝓁T_r^𝓁(t)}

=(

N_r^𝓁(t) −T_r+1^𝓁 (t) −T_r^𝓁(t))

𝛽rr^𝓁T^𝓁_r(t) − (

b^𝓁_(n−r+1)+𝜃r^𝓁

) T_r^𝓁(t) +

∑p

i=1

(1−𝛿𝓁i){

m_𝓁iT_rⁱ(t) −m_i𝓁T^𝓁_r(t)} (10)

and

dT_r^𝓁(t) dt =(

N_r^𝓁∗−T_r^𝓁(t))

𝛽rr^𝓁T_r^𝓁(t) − (

b^𝓁_(n−r)+𝜃r^𝓁

) T^𝓁_r(t) +

∑p 𝓁=1

(1−𝛿𝓁i){

m_𝓁iT_rⁱ(t) −m_i𝓁T_r^𝓁(t)} , 𝓁=1,2,…,p.

(11)

We know thatN_r^𝓁(t)(𝓁=1,…,p) converge toN_r^𝓁∗(𝓁=1,…,p), whereas from the definition ofr, we have that T_r+1^𝓁 (t)(𝓁=1,…,p) converge toT^𝓁∗_r+1(𝓁=1,…,p). Then, for any𝜀 >0, there exists āt>0 such that|N_r^𝓁(t)−N_r^𝓁∗|<

𝜀and|T_r+1^𝓁 (t) −T_r+1^𝓁∗ |< 𝜀for allt > ̄t,𝓁 = 1,…,p. If we substituteT_r+1^1∗ +𝜀,…,T_r+1^p∗ +𝜀,N_r^1∗−𝜀,…,N_n−q^p∗ −𝜀, resp.T_r+1^1∗ −𝜀,…,T_r+1^p∗ −𝜀,N_r^1∗+𝜀,…,N_n−q^p∗ +𝜀into (10), we obtain two systems of the same structure as (11), and one of them is a lower, the other is an upper estimate of (10), and each has a globally asymptotically stable equilibrium (T¹_r(𝜀),…,T^p_r(𝜀))

and(T¹r(𝜀),…,T^pr(𝜀)), respctively, because of Proposition 1. It is easy to see that the original system (10), considered as a nonautonomous system with time-dependent coefficientsT_r+1¹ (t),…,T_r+1^p (t),N_r¹(t),…,N_r^p(t), satisfies the quasimonotone condition, as well as the systems obtained after the substitution. Hence, we can apply Theorem 4.1 of Hirsch and Smith¹⁷to obtain that for any solution(T_r¹(t),…,T_r^p(t))of (10),

T^𝓁_r(𝜀)≤lim inf

t→∞ T_r^𝓁(t)≤lim sup

t→∞

T^𝓁_r(t)≤T^𝓁r(𝜀), 𝓁=1,…,p. (12)

Solutions of limit Eq. (11) converge to a globally asymptotically stable equilibrium by Proposition 1, and by letting 𝜀→0, we find that this limit is the same as that of (10).

As we have assumed that for all larger indices, the limits of the compartments of the original system (5) are equal to the limits obtained during the procedure, using the equations forT¹_r(t),…,T_r^p(t)aftern−r+1 cycles of the procedure satisfy the quasimonotone condition and the comparison (12), the limits obtained for these have to coincide with those of the original system (forr=n, the statement follows directly).

To prove that not only attractivity but also global asymptotic stability holds, we will again use induction. LetE= (S̄¹,T¹1,…,T¹n,…, ̄S^p,T^p1,…,T^pn)denote the equilibrium obtained at the end of the procedure, whereT^𝑗i = 0 or T^𝑗i > 0 depending on the stability moduli(s(M1),s(M2),…,s(Mn))and letE_𝜅 = (S̄¹,T¹1,…,T¹_𝜅,…, ̄S^p,T^p1,…,T^p_𝜅) be the equilibrium of thep(𝜅+1)-dimensional system

dT_k^𝓁(t) dt =

( N_𝜅^𝓁(t) −

∑𝜅 𝑗=1

T^𝓁_𝑗(t) )

𝛽_kk^𝓁T^𝓁_k(t) +T^𝓁_k(t)

∑𝜅 𝑗=1

(1−𝛿k𝑗)𝛽_k^𝓁_𝑗T_𝑗^𝓁(t) − (

b^𝓁_(n−𝜅)+𝜃_k^𝓁) T_k^𝓁(t)

+

∑p

i=1

(1−𝛿_𝓁i){

m_𝓁iT_kⁱ(t) −mi𝓁T^𝓁_k(t)} , k=1,2,…, 𝜅−1, 𝓁=1,2,…,p,

(13a) dT_𝜅^𝓁(t)

dt =T_𝜅^𝓁(t)

(𝛽_𝜅,𝜅^𝓁 N_𝜅^𝓁(t) − (

b^𝓁_(n−𝜅)+𝜃^𝓁_𝜅)

−𝛽_𝜅,𝜅^𝓁 T_𝜅^𝓁(t) )

+

∑p

i=1

(1−𝛿𝓁i){

m_𝓁iT_𝜅ⁱ(t) −m_i𝓁T^𝓁_𝜅(t)} , 𝓁=1,2,…,p,

(13b)

dN_𝜅^𝓁(t)

dt =B^𝓁_(n−𝜅)−b^𝓁_(n−𝜅)N_𝜅^𝓁(t) +

∑p

i=1

(1−𝛿_𝓁i){

m_𝓁iN_𝜅ⁱ(t) −m_i𝓁N_𝜅^𝓁(t)} , 𝓁=1,2,…,p,

(13c)

obtained during the procedure withE_𝜅consisting of the firstp(𝜅+1)coordinates ofE. Let us suppose thatE_𝜅is a stable equilibrium of thep(𝜅+1)-dimensional reduced system for some𝜅 ≤ n. We will show that in each step,E_𝜅+1

is a stable equilibrium of thep(𝜅+2)-dimensional reduced system dT^𝓁

k(t) dt =

(

N_𝜅+1^𝓁 (t) −

𝜅+1∑

𝑗=1

T_𝑗^𝓁(t) )

𝛽_kk^𝓁T_k^𝓁(t) +T_k^𝓁(t)

𝜅+1∑

𝑗=1

(1−𝛿k𝑗)𝛽_k𝑗^𝓁T_𝑗^𝓁(t)

− (

b^𝓁_{(n−𝜅−1)}+𝜃_k^𝓁) T_k^𝓁(t) +

∑p

i=1

(1−𝛿𝓁i){

m_𝓁iT_kⁱ(t) −mi𝓁T^𝓁_k(t)} , k=1,2,…, 𝜅,

(14a)

dT_𝜅+1^𝓁 (t)

dt =

(

N_𝜅^𝓁₊₁(t) −

𝜅+1∑

𝑗=1

T_𝑗^𝓁(t) )

𝛽_𝜅^𝓁_+1,𝜅+1T^𝓁_𝜅+1(t) +T^𝓁_𝜅+1(t)

𝜅+1∑

𝑗=1

(1−𝛿𝜅+1,𝑗)𝛽_𝜅^𝓁_+1,𝑗T_𝑗^𝓁(t)

− (

b^𝓁_{(n−𝜅−1)}+𝜃_𝜅+1^𝓁 )

T_𝜅+1^𝓁 (t) +

∑p

i=1

(1−𝛿_𝓁i){

m_𝓁iT_𝜅+1ⁱ (t) −mi𝓁T^𝓁_𝜅+1(t)} ,

=T_𝜅^𝓁₊₁(t)

(𝛽_𝜅^𝓁_+1,𝜅+1N_𝜅^𝓁₊₁(t) −𝛽_𝜅^𝓁_+1,𝜅+1T_𝜅^𝓁₊₁(t) − (

b^𝓁_{(n−𝜅−1)}+𝜃_𝜅^𝓁₊₁)) +

∑p

i=1

(1−𝛿𝓁i){

m_𝓁iTⁱ_𝜅+1(t) −m_i𝓁T_𝜅^𝓁₊₁(t)} ,

(14b)

dN_𝜅+1^𝓁 (t)

dt =B^𝓁_{(n−𝜅−1)}−b^𝓁_{(n−𝜅−1)}N_𝜅+1^𝓁 (t) +

∑p

i=1

(1−𝛿_𝓁i){

m_𝓁iN_𝜅+1ⁱ (t) −mi𝓁N_𝜅+1^𝓁 (t)} , 𝓁=1,2,…,p.

(14c)

Suppose this does not hold, that is,E_𝜅+1 is unstable. In this case, there exist an 𝜀 > 0 and a sequence {x_m} → E_𝜅+1, |xm − E_𝜅+1| < 1∕m such that the orbits started from the points of the sequence leave B(E_𝜅+1, 𝜀) ∶=

{

x∈R^(𝜅+2)p₊ ^∶|x−E_𝜅+1|≤𝜀}

. By an exit point fromB(E_𝜅+1, 𝜀), we mean a pointxsuch that|E_𝜅+1−x| = 𝜀and for the trajectory throughx, there is an open intervalJ ∋ 0 such that for allt ∈ J,xt ∈ B(E_𝜅+1, 𝜀)ift ≤ 0 and xt∉B(E_𝜅+1, 𝜀)ift>0. Let us denote byx^𝜀_mthe first exit point fromB(E_𝜅+1, 𝜀)of the solution started fromx_m, reached at time𝜏m. There is a convergent subsequence of the sequencex^𝜀_m(still denoted byx_m^𝜀) which tends to a point denoted byx_𝜀^∗ ∈S(E_𝜅+1, 𝜀) ∶={

x∈R^(𝜅+2)p₊ ∶|x−E_𝜅+1|=𝜀}

. We will show thatE_𝜅+1 ∈𝛼(x^∗_𝜀). For this end, let us consider the setS

(

E_𝜅+1,^𝜀₂)

. Clearly, all solutions started from the pointsxm(we drop the first elements of the sequence, if necessary) will leave the setB

(

E_𝜅+1,^𝜀₂)

. We denote the last exit point of each trajectory from this set before time𝜏m, respectively, byx_m^𝜀∕2. Also this sequence has a convergent subsequence (still denoted the same way), let us denote its limit byx_𝜀∕2^∗ . We will show that the trajectory started fromx_𝜀∕2^∗ goes throughx_𝜀^∗. AsE_𝜅+1is globally attractive, this tra- jectory will eventually enterS

(

E_𝜅+1,₄^𝜀)

at some timeT>0. Let us suppose that the trajectory started fromx^∗_𝜀∕2does not go throughx_𝜀^∗and let us denote byd>0 the distance of this trajectory fromx^∗_𝜀. For continuity reasons, there is an N∈Nso that for anym>N,|x^∗_𝜀∕2t−x^𝜀∕2_m t|<max

{d 2,^𝜀₈}

for 0<t<T. This means that formlarge enough, the tra- jectory started fromx^𝜀∕2_m will enter againS

(

E_𝜅+1,^𝜀₂)

without getting close tox_𝜀^∗which contradicts eitherx^𝜀_mbeing the first exit point fromB(E_𝜅+1, 𝜀)orx^𝜀_m^∕2being the last exit point before𝜏mfromB

(

E_𝜅+1,^𝜀₂)

. Hence, we have shown that the trajectory started fromx_𝜀^∗_∕2goes throughx^∗_𝜀. Proceeding like this (taking neighbourhoods of radius𝜀∕4, 𝜀∕8 etc.), we obtain that the backward trajectory ofx^∗_𝜀enters any small neighbourhood ofE_𝜅+1. That is, there exists a decreas- ing sequencetn < 0 such that|xtn−E_𝜅+1| < _2n^𝜀. We have eithertn → t^* for somet^* < 0 ortn → −∞. In the first case,xt_n → xt^* = E_𝜅+1which contradicts the fact thatE_𝜅+1is an equilibrium. Hence,t_n → ∞, and we obtain that E_𝜅+1∈𝛼(x^∗_𝜀), while it follows from the global attractivity ofE_𝜅+1that the𝜔-limit set of the trajectory is{E_𝜅+1}. Let us denote this trajectory by𝛾(x^∗_𝜀)

We know that Equations (14b) for ^d

dtT_𝜅¹₊₁(t),…,_dt^dT_𝜅^p₊₁(t)and (14c) for ^d

dtN_𝜅¹₊₁(t),…,_dt^dN_𝜅^p₊₁(t)can be decoupled from the rest of the equations and using the exponential stability of the limits

t→+∞lim N_𝜅^𝓁₊₁(t) =N_𝜅^𝓁∗₊₁,𝓁=1,2,…,p,

and Proposition 1 we obtain thatT¹_𝜅+1,…,T^p_𝜅+1is a stable equilibrium of the system consisting of the system dT_𝜅+1^𝓁 (t)

dt =T_𝜅^𝓁₊₁(t)(

𝛽_𝜅^𝓁_+1,𝜅+1N_𝜅^𝓁∗₊₁−(

b^𝓁+𝜃_𝜅^𝓁₊₁)

−𝛽_𝜅^𝓁_+1,𝜅+1T_𝜅^𝓁₊₁(t)) +

∑p

i=1

(1−𝛿𝓁i){

m_𝓁iT_𝜅ⁱ₊₁(t) −m_i𝓁T^𝓁_𝜅+1(t)}

,𝓁=1,2,…,p.

Therefore, the equilibriumE_𝜅+1is stable in the coordinatesT_𝜅+1¹ ,…,T_𝜅+1^p in the sense that for anỹ𝜀 >0 there exists a ̃𝛿(̃𝜀) >0 such that for any initial valuexwith|x−E_𝜅+1| < ̃𝛿,|T_𝜅+1^𝓁 (t) −T^𝓁_𝜅+1| < ̃𝜀for allt > 0 and𝓁 = 1,…,p.

Thus, the trajectory𝛾(x_𝜀^∗)obtained above lies entirely in the subspace {

T_𝜅+1¹ =T¹_𝜅+1,…,T^p_𝜅+1=T^p_𝜅+1 }

. On the other hand, the currentp(𝜅+2)-dimensional system coincides with thep(𝜅+1)-dimensional system on this subspace.

For the latter system, stability of the equilibriumE_𝜅follows from the induction assumption. However, the existence of an orbit̃𝛾 (different from the equilibriumE_𝜅+1) whose𝜔-limit set is{E_𝜅+1}and whose𝛼-limit set containsE_𝜅+1 contradicts the stability of the equilibriumE_𝜅. Indeed, let us suppose thatE_𝜅is stable and there exists such an orbit

̃𝛾. The stability ofE_𝜅would implies that for any ̂𝜀 >0, there exists a ̂𝛿(̂𝜀)such that for any solution started from an initial valuePwith|P−E_𝜅|< ̂𝛿, we have|Pt−E_𝜅|< ̂𝜀for allt>0. Hence, this is also true for̂𝜀=|E_𝜅−P|∕2 for anỹ P̃ ∈ ̃𝛾, which is a contradiction, as a solution started from a point of̃𝛾clearly leavesS(E_𝜅, ̂𝜀). Hence, no such orbit ̂𝛾 can exist. This implies the global asymptotic stability of the equilibrium of thep(𝜅+2)-dimensional system.