1. Introduction

Many pathogenic microorganisms have several different genetic variants or subtypes which are called strains. Different strains competing for the same host may differ in their key epidemiological parameters, such as infectivity, length of infectious period or virulence. General multistrain models are typically difficult to analyse because of the large state space (see Kryazhimskiy et al. [8]), sometimes even showing chaotic dynamics, as in Bianco et al. [1]. Bichara et al. [2] considered multi-strain SIS and SIR models without superinfection and proved that under generic conditions a competitive exclusion principle holds. Martcheva [9] showed that in a periodic environment, the principle of competitive exclusion does not necessarily hold. In reality, a stronger strain might superinfect an individual already infected by another strain. As a consequence, different virus strains with different virulence may coexist even in a constant environment, as it has been illustrated by Nowak [10], who considered a basic model to provide an analytical understanding of the complexities introduced by superinfection.

In the present paper, we consider a multistrain SIS model with superinfection with $n$ infectious strains and show that it is possible to obtain a stable coexistence of any subgroup of the $n$ strains. We establish an iterative method for calculating a sequence of reproduction numbers which determine what strains are present in the globally asymptotically stable coexistence equilibrium.

The structure of the paper is the following: in Section 2, we introduce a multistrain SIS model with superinfection. We present the iterative method which can be used to determine which equilibrium of the system will be the stable coexistence steady state, and we prove that it is globally asymptotically stable. In the last section, we apply the method described in Section 2 in the case of three strains. For this case, we give a complete description of the global stability properties of the system depending on the different reproduction numbers.

---

2. Main result

We consider a heterogeneous virus population with different infectivities. We will assume that superinfection is possible, i.e. more infective strains outcompete the less infective ones in an infected individual. We assume that an infected individual is always infected by only one virus strain, i.e. after superinfection, the more infective strain completely takes over the host from the less infective one. Let $n$ denote the number of strains with different infectivity. The population is divided into $n+1$ compartments depending on the presence of any of the virus strains: the susceptible class is denoted by $S(t)$ and we have $n$ infected compartments $T_1,\dots,T_n$, where a larger index corresponds to a compartment of individuals infected by a strain with larger infectivity. Let $B$ denote birth rate and $b$ death rate. Let $\beta_{jk}$ denote the transmission rate by which the $k$-th strain infects those who are infected by the $j$-th strain. We refer to the transmission rates from susceptibles to strain $k$ by $\beta_{kk}$. The notation $\theta_k$ stands for recovery rate among those infected by the $k$-th strain. We allow the most infective strain to be lethal with disease-induced mortality rate $d_n$. Using these notations, we obtain the following model:

$\begin{aligned} \frac{\mathrm{d}S(t)}{\mathrm{d}t}={}&B-b S(t)- S(t) \sum\limits_{k=1}^n \beta_{kk} T_k(t)+\sum\limits_{k=1}^n \theta_k T_k(t), \\ \frac{\mathrm{d}T_k(t)}{\mathrm{d}t}={}&S(t) \beta_{kk} T_k(t)+T_k(t) \sum\limits_{j=1}^n (1-\delta_{kj}) \beta_{kj} T_j(t)\\ &-(b+\theta_k+\delta_{kn}d_n) T_k(t), k=1,2,\dots,n, \end{aligned}$

(1)

with initial condition

$\begin{aligned} & S\left( 0\right) =\phi_0, T_k\left( 0\right) =\phi_k, k=1,2,\dots, n, \\ &\left( \phi_0,\phi_1,\phi_2,\dots,\phi_n \right) \in \Gamma, \end{aligned}$

(2)

where $\delta_{kj}$ denotes the Kronecker delta such that $\delta_{kj}=1$ if $k=j$ and $\delta_{kj}=0$ otherwise, and $\Gamma=[0,\infty)^{n+1}$. We assume that the conditions

$\begin{aligned} \beta_{kj}&=\beta_{kk},&1 &\leq j \leq k, \\ \beta_{jk}&=-\beta_{kj}=-\beta_{jj},&k+1 &\leq j \leq n, \end{aligned}$

(3)

hold for the infection rates for $k=1,2,\ldots,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.

Note that for $n=2$, (1) is equivalent to the model by Dénes and Röst [4,5,6] describing the spread of ectoparasites and ectoparasite-borne diseases. In [4,5], a model with no disease-induced mortality is studied (this corresponds to $d_n=0$ in the present model), while [6] studies the impact of excess mortality ($d_n>0$).

Before we present the procedure for the global stability analysis of the system (1), we note that for a given $n$, it is enough to consider solutions started with initial values $\phi_1,\dots,\phi_n>0$, as any boundary subspace of $\Gamma$ (i.e. a subspace where one or more infected compartments are equal to 0) is invariant, so if any $m$ of the initial values $\phi_1,\dots,\phi_n$ is equal to $0$, we can reduce (1) to an $(n+1-m)$-dimensional system of the same structure. This also means that in the case of $n$ strains, solutions started from the boundary of $\Gamma$ can be studied in an analogous way as solutions started from the interior, but in smaller dimension. This observation is formalized in the following proposition.

Proposition 1. Let us consider the system (1) with $n$ different strains, and let $m<n$. On any $m+1$-dimensional boundary subspace of $\Gamma$ (i.e.when $n-m$ strains are not present), (1) can be reduced to an $(m+1)$-dimensional system with the same structure. The corresponding boundary subspace is invariant for the reduced equation.

Now we can start the description of the procedure for the global stability analysis of (1); from now on we only consider solutions started with initial values $\phi_1,\dots,\phi_n>0$. Let us introduce the new variable

$N_n(t)=S(t)+\sum\limits_{j=1}^n T_j(t)$

(4)

representing the total population. By (3), $\beta_{kj}=-\beta_{jk}\ \textrm{for} \ k \neq j$, and hence,

$\displaystyle{\sum\limits_{k=1}^n T_k(t) \sum\limits_{j=1}^n (1-\delta_{kj}) \beta_{kj} T_j(t)=0}$

holds. By (1), we have

$\displaystyle{\frac{\mathrm{d}N_n(t)}{\mathrm{d}t}=B-b N_n(t)-d_n T_n(t)},$

(5)

and (1) can be rewritten as follows:

$\begin{aligned} \frac{\mathrm{d}S(t)}{\mathrm{d}t}={}&B-b S(t)- S(t) \sum\limits_{k=1}^n \beta_{kk} T_k(t)+\sum\limits_{k=1}^n \theta_k T_k(t), \\ \frac{\mathrm{d}T_k(t)}{\mathrm{d}t}={}&S(t) \beta_{kk} T_k(t)+T_k(t)\sum\limits_{j=1}^n (1-\delta_{kj}) \beta_{kj} T_j(t)\\ &-(b+\theta_k) T_k(t), k=1,2,\ldots,n-1, \\ \end{aligned}$

(6)

and

$\begin{aligned} \frac{\mathrm{d}T_n(t)}{\mathrm{d}t}={}&\biggl(N_n(t)-\sum\limits_{j=1}^n T_j(t)\biggr) \beta_{nn} T_n(t)+T_n(t)\sum\limits_{j=1}^n (1-\delta_{nj}) \beta_{nj} T_j(t)\\ &-(b+\theta_n+d_n) T_n(t) \\ ={}&T_n(t)\biggl(\beta_{nn}N_n(t)-\sum\limits_{j=1}^n \{\beta_{nn}-(1-\delta_{nj})\beta_{nj}\} T_j(t)-(b+\theta_n+d_n)\biggr), \\ \frac{\mathrm{d}N_n(t)}{\mathrm{d}t}={}&B-b N_n(t)-d_n T_n(t). \end{aligned}$

(7)

Then, by (3), introducing the new variable

$U_n(t)=B/b-N_n(t), t \geq 0,$

(7) is equivalent to the following system:

$\begin{aligned} \frac{\mathrm{d}T_n(t)}{\mathrm{d}t}={}&\beta_{nn}T_n(t)\biggl(\frac{B}{b}-\frac{b+\theta_n+d_n}{\beta_{nn}}-T_n(t)-U_n(t)\biggr), \\ \frac{\mathrm{d}U_n(t)}{\mathrm{d}t}={}&d_n T_n(t)-b \, U_n(t). \end{aligned}$

(8)

Let us define $\mathcal{R}_0^{(n)}$ as

$\mathcal{R}_{0}^{(n)}\ :\text{=}\frac{B{{\beta }_{nn}}}{b(b+{{d}_{n}}+{{\theta }_{n}})}.$

Since $B/b$ is the number of susceptibles at the disease free equilibrium of system (1), this formula can be interpreted as the basic reproduction number of the $n$-th strain.

The system (8) can be decoupled as a 2-dimensional Lotka-Volterra system with feedback controls (see, for example, Faria and Muroya [7]) and is independent from the system (6). We note that one can show global stability for (8) by applying the general result of Faria and Muroya [7,Theorem 3.8] for Lotka-Volterra systems with feedback controls. However, the application of that technique to our system (8) requires the condition $d_n<b$. Here, without the condition $d_n<b$, we give a proof of the global stability of (8). First, we show a simple theorem on the global dynamics of (8) to be applied later in determining the global dynamics of the full system.

Theorem 2.1. Let us consider system (8) on $\mathbb{R}^0_+\times\mathbb{R}$. For $\mathcal{R}_0^{(n)} \leq 1$, the system has only the trivial equilibrium $(0,0)$, which is globally asymptotically stable. For $\mathcal{R}_0^{(n)}>1$, system (8) has two equilibria: the trivial equilibrium $(0,0)$ and the globally asymptotically stable positive equilibrium

$(T_n^*,U_n^*)=\left(\frac{B\beta_{nn}-(b+d_n+\theta_n)b}{\beta_{nn}(b+d_n)},\frac{d_n\{B\beta_{nn}-(b+d_n+\theta_n)b\}}{b\beta_{nn}(b+d_n)}\right),$

which only exists if $\mathcal{R}_0^{(n)}>1$.

Proof. Obviously, if we start a solution with $T_n(0)=0$ then $U_n(t)\to 0$ as $t\to\infty$, so all such solutions tend to the trivial equilibrium.

We show that for $\mathcal{R}_0^{(n)}>1$, the limit set of all other solutions is the positive equilibrium. Let us suppose this is not true, i.e. there is a solution started from positive initial value in $T$ which tends to the trivial equilibrium. If this holds, then for any $\varepsilon>0$, there is a $T>0$ such that $T_n(t)<\varepsilon$ and $|U_n(t)|<\varepsilon$ holds for all $t>T$. Then for the derivative $T_n'(t)$, we can give the estimation

$\begin{align*} T_n'(t)&=\beta_{nn}T_n(t)\biggl(\frac{B}{b}-\frac{b+\theta_n+d_n}{\beta_{nn}}-T_n(t)-U_n(t)\biggr)\\ &>\beta_{nn}T_n(t)\left(\frac{B}{b}-\frac{b+\theta_n+d_n}{\beta_{nn}}-2\varepsilon\right), \end{align*}$

which is positive for $\varepsilon$ small enough as $\tfrac{B}{b}>\tfrac{b+\theta_n+d_n}{\beta_{nn}}$ follows from $\mathcal{R}_0^{(n)}>1$. This contradicts our assumption, thus for $\mathcal{R}_0^{(n)}>1$, the limit of any solution started from positive initial values is the positive equilibrium. Let us note that all solutions are bounded. We apply the Dulac function $D(T_n,U_n)=1/T_n$ to obtain

$\begin{align*} &\frac{\partial}{\partial T_n}\left[D(T_n,U_n)\biggl\{\beta_{nn}T_n\biggl(\frac{B}{b}-\frac{b+\theta_n+d_n}{\beta_{nn}}-T_n-U_n\biggr)\biggr\}\right]\\ &\mathrel{\phantom{=}}+\frac{\partial}{\partial U_n}\{D(T_n,U_n)\left(d_n T_n-b \, U_n\right)\}\\ & = -\beta_{nn}-\frac{b}{T_n}, \end{align*}$

which is negative for $T_n>0$. Hence, we obtain from the above and the Poincaré-Bendixson theorem that for $\mathcal{R}_0^{(n)} \leq 1$, the trivial equilibrium $(0,0)$ is globally attractive and for $\mathcal{R}_0^{(n)}>1$, there exists a unique positive equilibrium of (8) which is globally attractive on $\{\,(T_n,U_n)\in\mathbb{R}_+\times\mathbb{R}\,\}$. Stability follows from the fact that Dulac's criterion excludes homoclinic orbits too (cf. [3]).

To determine the global dynamics of (1), which is equivalent to the global dynamics of (6) and (8), first we substitute the equilibria calculated in Theorem 2.1 into the full model to obtain the following reduced system of (6):

$\begin{aligned} \frac{\mathrm{d}S(t)}{\mathrm{d}t}={}&B^{(1)}-b^{(1)} S(t)- S(t) \sum\limits_{k=1}^{n-1} \beta_{kk} T_k(t)+\sum\limits_{k=1}^{n-1} \theta_k T_k(t), \\ \frac{\mathrm{d}T_k(t)}{\mathrm{d}t}={}& S(t) \beta_{kk} T_k(t)+T_k(t) \sum\limits_{j=1}^{n-1} (1-\delta_{kj}) \beta_{kj} T_j(t)\\ & -\left(b^{(1)}+{\theta}_k \right) T_k(t), k=1,2,\dots,n-2, \\ \end{aligned}$

(9)

and

$\begin{aligned} \frac{\mathrm{d}T_{n-1}(t)}{\mathrm{d}t}={}&\biggl(N_{n-1}(t)-\sum\limits_{j=1}^{n-1} T_j(t)\biggr) \beta_{n-1,n-1} T_{n-1}(t) \\ &+T_{n-1}(t)\sum\limits_{j=1}^{n-1} (1-\delta_{n-1,j}) \beta_{n-1,j} T_j(t)-(b^{(1)}+\theta_{n-1}) T_{n-1}(t) \\ ={}&T_{n-1}(t)\biggl(\beta_{n-1,n-1}N_{n-1}(t) \\ &\phantom{T_{n-1}(t)\biggl(}-\sum\limits_{j=1}^{n-1} \{\beta_{n-1,n-1}-(1-\delta_{n-1,j})\beta_{n-1,j}\} T_j(t)\\ &\phantom{T_{n-1}(t)\biggl(}-(b^{(1)}+\theta_{n-1})\biggr), \\ \frac{\mathrm{d}N_{n-1}(t)}{\mathrm{d}t}={}&B^{(1)}-b^{(1)} N_{n-1}(t), \end{aligned}$

(10)

where $N_{n-1}$ is defined similarly as $N_n(t)$ in (4), as the sum of the susceptible compartment and the first $n-1$ (already modified) infected compartments:

$N_{n-1}(t)=S(t)+\sum\limits_{k=1}^{n-1}T_{k}(t)$

and the new coefficients are defined as

$B^{(1)}:= B+\theta_n T_n^*, b^{(1)}:= b+\beta_{nn}T_n^*,$

if $\mathcal{R}_0^{(n)}>1$ and

$B^{(1)}:= B, b^{(1)}:= b, \\ $

if $\mathcal{R}_0^{(n)} \leq 1$. We define the next reproduction number as

$\mathcal{R}_0^{(n-1)}:= \frac{B^{(1)} \beta_{n-1,n-1}}{b^{(1)}(b^{(1)}+\theta_{n-1})}.$

It is easy to see that the two systems (9) and (10) we obtained are of similar structure as (6) and (7), but with dimension $n-1$, the positivity of the new parameters follows from the conditions (3). This means that by repeating the above steps we can further reduce the dimension by substituting the values of the globally asymptotically equilibrium of the decoupled system (10) into the remaining equations.

In general, after performing the above steps $l$ times, we arrive at the system

$\begin{aligned} \frac{\mathrm{d}S(t)}{\mathrm{d}t}={}&B^{(l)}-b^{(l)} S(t)- S(t) \sum\limits_{k=1}^{n-l} \beta_{kk} T_k(t)+\sum\limits_{k=1}^{n-l} \theta_k T_k(t), \\ \frac{\mathrm{d}T_{k}(t)}{\mathrm{d}t}={}& S(t) \beta_{kk} T_k(t)+T_k(t) \sum\limits_{j=1}^{n-l} (1-\delta_{kj}) \beta_{kj} T_j(t)\\ & -\left(b^{(l)}+{\theta}_k\right) T_k(t), k=1,2,\dots,n-l-1, \\ \end{aligned}$

(11)

and

$\begin{aligned} \frac{\mathrm{d}T_{n-l}(t)}{\mathrm{d}t}={}&S(t) \beta_{n-l,n-l} T_{n-l}(t)+T_{n-l}(t)\sum\limits_{j=1}^{n-l} (1-\delta_{n-l,j}) \beta_{n-l,j} T_j(t)\\ &-\left(b^{(l)}+{\theta}_{n-l}\right) T_{n-l}(t), \\ \frac{\mathrm{d}N_{n-l}(t)}{\mathrm{d}t}={}&B^{(l)}-b^{(l)} N_{n-l}(t), \end{aligned}$

(12)

where

$N_{n-l}(t)=S(t)+\sum\limits_{k=1}^{n-l}T_k(t),$

$B^{(0)}=B, \ b^{(0)}=b$ and inductively for $l=1,2,\ldots,n-1$, we define

$\mathcal{R}_0^{(n-l)}:= \frac{B^{(l)} \beta_{n-l,n-l}}{b^{(l)}(b^{(l)}+\theta_{n-l})}$

and

$B^{(l)}:= B^{(l-1)}+\theta_{n-l+1} T_{n-l+1}^*, b^{(l)}:= b^{(l-1)}+\beta_{n-l+1,n-l+1}T_{n-l+1}^*,$

if $\mathcal{R}_0^{(n-l)}>1$ and

$B^{(l)}:= B^{(l-1)}, b^{(l)}:= b^{(l-1)},$

if $\mathcal{R}_0^{(n-l)} \leq 1$.

Now we introduce $U_{n-l}(t)=B^{(l)}/b^{(l)}-N_{n-l}(t)$, to rewrite the equation (12) as

$\begin{aligned} \frac{\mathrm{d}T_{n-l}(t)}{\mathrm{d}t}={}&\beta_{n-l,n-l}T_{n-l}(t)\left(\frac{B^{(l)}}{b^{(l)}}-\frac{b^{(l)}+{\theta}_{n-l}}{\beta_{n-l,n-l}}-T_ {n-l}(t)-U_{n-l}(t)\right), \\ \frac{\mathrm{d}U_{n-l}(t)}{\mathrm{d}t}={}&-b^{(l)} \ U_{n-l}(t). \end{aligned}$

(13)

Again, (13) might be decoupled from the other equations (11). For $\mathcal{R}_0^{(n-l)} \leq 1$, system (13) has only the trivial equilibrium $(0,0)$. But for $\mathcal{R}_0^{(n-l)}>1$, system (13) has two equilibria: the trivial equilibrium $(0,0)$ and the non-trivial equilibrium

$(T_{n-l}^*,U_{n-l}^*)=\left(\frac{B^{(l)}\beta_{n-l,n-l}-(b^{(l)}+\theta_{n-l})b^{(l)}}{\beta_{n-l,n-l}b^{(l)}},0 \right),$

which only exists if

$\mathcal{R}_0^{(n-l)}>1.$

Then, from (11), we obtain the systems

$\label{lsteps-3} \begin{aligned} \frac{\mathrm{d}S(t)}{\mathrm{d}t}={}&B^{(l+1)}-b^{(l+1)} S(t)- S(t) \sum\limits_{k=1}^{n-l-1} \beta_{kk} T_k(t)+\sum\limits_{k=1}^{n-l-1} \theta_k T_k(t), \\ \frac{\mathrm{d}T_{k}(t)}{\mathrm{d}t}={}& S(t) \beta_{kk} T_k(t)+T_k(t) \sum\limits_{j=1}^{n-l-1} (1-\delta_{kj}) \beta_{kj} T_j(t)\\ & -\left(b^{(l+1)}+{\theta}_k\right) T_k(t), k=1,2,\dots,n-l-2, \\ \end{aligned}$

and

$\label{lsteps-4} \begin{aligned} \frac{\mathrm{d}T_{n-l-1}(t)}{\mathrm{d}t}={}&S(t) \beta_{n-l-1,n-l-1} T_{n-l-1}(t)\\ &+T_{n-l-1}(t)\sum\limits_{j=1}^{n-l-1} (1-\delta_{n-l-1,j}) \beta_{n-l-1,j} T_j(t)\\ &-\left(b^{(l+1)}+{\theta}_{n-l-1}\right) T_{n-l-1}(t), \\ \frac{\mathrm{d}N_{n-l-1}(t)}{\mathrm{d}t}={}&B^{(l+1)}-b^{(l+1)} N_{n-l-1}(t), \end{aligned}$

where

$N_{n-l-1}(t)=S(t)+\sum\limits_{k=1}^{n-l-1}T_k(t),$

$B^{(0)}=B, \ b^{(0)}=b$ and inductively for $l=1,2,\ldots,n-1$, we define

$\mathcal{R}_0^{(n-l-1)}:= \frac{B^{(l+1)} \beta_{n-l-1,n-l-1}}{b^{(l+1)}(b^{(l+1)}+\theta_{n-l-1})}$

and

$B^{(l+1)}:= B^{(l)}+\theta_{n-l} T_{n-l}^*, b^{(l+1)}:= b^{(l)}+\beta_{n-l,n-l}T_{n-l}^*,$

if $\mathcal{R}_0^{(n-l-1)}>1$ and

$B^{(l+1)}:= B^{(l)}, b^{(l+1)}:= b^{(l)},$

if $\mathcal{R}_0^{(n-l-1)} \leq 1$, which, again, are systems with the same structure. In the end, we arrive at the two-dimensional system

$\begin{split} \frac{\mathrm{d}S(t)}{\mathrm{d}t}={}& B^{(n-1)}-b^{(n-1)}S(t)-\beta_{11}S(t)T_1(t)+\theta_1T_1(t),\\ \frac{\mathrm{d}T_1(t)}{\mathrm{d}t}={}&\beta_{11}S(t)T_1(t)-(b^{(n-1)}+\theta_1)T_1(t), \end{split}$

which has the two equilibria

$\left(\frac{B^{(n-1)}}{b^{(n-1)}},0\right)\textrm{and}\left(\frac{b^{(n-1)}+\theta_1}{\beta_{11}},\frac{B^{(n-1)}}{b^{(n-1)}}-\frac{b^{(n-1)}+\theta_1}{\beta_{11}}\right),$

with the latter one only existing if

$\mathcal{R}_0^{(n)}:=\frac{B^{(n-1)}\beta_{11}}{b^{(n-1)}(b^{(n-1)}+\theta_1)}>1.$

The dynamics of this system can be determined in a similar way as in the case of (8), and we obtain that the first equilibrium is globally asymptotically stable if $\mathcal{R}_0^{(n)}\leq 1$ and the second one is globally asymptotically stable if $\mathcal{R}_0^{(n)}> 1$.

Thus, by the above discussion, we can reach a conclusion by induction to the global dynamics of the model (1) and we claim the following.

Theorem 2.2. The multistrain SIS model (1) has a globally asymptotically stable equilibrium on the region $\Gamma_0$, where $\Gamma_0$ is the interior of $\Gamma$. The global dynamics is completely determined by the threshold parameters $\big(R_0^{(1)}, R_0^{(2)},\ldots,R_0^{(n)}\big)$ which can be obtained iteratively and determine which one of the equilibria is globally asymptotically stable.

Proof. The procedure described above can be used to prove the theorem. To show that at each step, after obtaining the globally asymptotically stable equilibrium of the two-dimensional system (13) we can really substitute the coordinates of this equilibrium into the remaining equations, we use the theory of asymptotically autonomous systems [11,Theorem 1.2]. According to this theorem, if $e$ is a locally asymptotically stable equilibrium of the system

$\dot{y}=g(y)$

(14)

which is the limit equation of the asymptotically autonomous equation

$\dot{x}=f(t,x),$

(15)

and $\omega$ is the $\omega$-limit set of a forward bounded solution $x$ of (15), then, if $\omega$ contains a point $y_0$ such that the solution of (14) through $(0,y_0)$ converges to $e$ for $t\to\infty$, then $\omega=\{e\}$, i.e. $x(t)\to e$ for $t\to\infty$.

Let us suppose that at the end of the procedure, we obtain the equilibrium $E=(\bar{S},\bar{T}_1,\dots,\bar{T}_n)$ where $\bar{T}_i=0$ or $\bar{T}_i=T_i^*$ depending on the reproduction numbers and let $E_k=(\bar{S},\bar{T}_1,\dots,\bar{T}_k)$ the equilibrium of the $(k+1)$-dimensional system obtained during the procedure, consisting of the first $k+1$ coordinates of $E$. Let $\mathcal{B}(M)$ denote the basin of attraction of a point $M$. Let us define the domains $\Gamma_0^k$ as $\Gamma_0^k :=\{\,(S,T_1,\dots,T_k,):S,T_1,\dots,T_k>0\,\}$ for $k=1,\dots,n-1$ and $\Gamma_0^n:=\Gamma_0$. By induction, we will show that $\Gamma_0\subset\mathcal{B}(E)$. Let us suppose that for a given $k$, $\Gamma_0^k\subset\mathcal{B}(E_k)$ holds and $E_k$ is stable. We will show that from this also $\Gamma_0^{k+1}\subset\mathcal{B}(E_{k+1})$ follows. We know from the procedure that on $\Gamma_0^{k+1}$, for all solutions of the $(k+2)$-dimensional system, the coordinate $T_{k+1}(t)$ tends to $\bar{T}_{k+1}$. Thus, for the limit set $\omega(x)$ of any $x\in\Gamma_0^{k+1}$ for the $(k+2)$-dimensional system, we have $\omega(x)\subset\Gamma_0^k\times\{\bar{T}_{k+1}\}$. Since $E_k$ was supposed to be stable, according to [11,Theorem 1.2], from this it follows that $x\in\mathcal{B}(E_{k+1})$, and, as $x$ was chosen arbitrarily, $\Gamma_0^{k+1}\subset\mathcal{B}(E_{k+1})$. Thus, we have shown global attractivity of the equilibrium.

To prove global asymptotic stability, and in order to be able to proceed with the induction using [11,Theorem 1.2], we still need to show that $E_{k+1}$ is a stable equilibrium of the $(k+1)$-dimensional reduced system in each step. Let us suppose that this does not hold, i.e. $E_{k+1}$ is unstable for some $k\leq n$. This means that there exists an $\varepsilon>0$ such that there exists a sequence $\{x_m\}\to 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},\varepsilon):= \{\,x\in\Gamma : |x-E_{k+1}|\leq \varepsilon\,\}$. Let us denote by $x_m^{\varepsilon}$ the first exit point from $B(E_{k+1},\varepsilon)$ of the solution started from $x_m$, reached at time $\tau_m$. There is a convergent subsequence of the sequence $x_m^{\varepsilon}$ (still denoted by $x_m^{\varepsilon}$) which tends to a point denoted by $x_{\varepsilon}^*\in S(E_{k+1},\varepsilon):= \{\,x\in\Gamma : |x-E_{k+1}|= \varepsilon\,\}$. We will show that the $\alpha$-limit set $\alpha(x_{\varepsilon}^*)$ is the singleton $\{E_{k+1}\,\}$. For this end, let us consider the set $S(E_{k+1},\frac{\varepsilon}{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},\frac{\varepsilon}{2})$. Let us denote the last exit point of each trajectory from this set before time $\tau_m$, respectively, by $x_m^{{\varepsilon/2}}$. Also this sequence has a convergent subsequence (still denoted the same way), let us denote its limit by $x_{{\varepsilon/2}}^*$. We will show that the trajectory started from this point will leave $B(E_{k+1},\varepsilon)$. If this is not the case, let us denote by $d>0$ the distance of this trajectory from $S(E_{k+1},\varepsilon)$. As $E_{k+1}$ is globally attractive, this trajectory will eventually enter $S(E_{k+1},\frac{\varepsilon}{4})$ at some time $T>0$. For continuity reasons, there exists $N\in\mathbb{N}$ such that if $m>N$ then $|x_{{\varepsilon/2}}^*t-x_m^{{\varepsilon/2}}t|<\max\{\frac{d}{2},\frac{\varepsilon}{8}\}$ for $0<t<T$. This means that for $m$ large enough, the trajectory started from $x_m^{{\varepsilon/2}}$ will enter again $S(E_{k+1},\frac{\varepsilon}{2})$ before exiting $S(E_{k+1},\varepsilon)$. This contradicts the choice of $x_m^{{\varepsilon/2}}$ as the last exit points from $S(E_{k+1},\frac{\varepsilon}{2})$. Again, continuity arguments show that the intersection point of the trajectory started from $x_{{\varepsilon/2}}^*$ and $S(E_{k+1},\varepsilon)$ is $x_{{\varepsilon}}^*$. Proceeding like this (taking neighbourhoods of radius $\varepsilon/4,\ \varepsilon/8$ etc.) we can show that the backward trajectory of $x_{{\varepsilon}}^*$ will enter any small neighbourhood of $E_{k+1}$ as $t\to-\infty$. From the above, it is also clear that $E_{k+1}$ is the only $\alpha$-limit point of this trajectory. Because of the global attractivity of $E_{k+1}$, the $\omega$-limit set of the trajectory is also $\{E_{k+1}\,\}$, thus the orbit is homoclinic.

As the equilibrium of the two-dimensional system (13) for $n-l=k+1$ is globally asymptotically stable, for any $\hat{\varepsilon}>0$, there exists an $\hat N\in\mathbb{N}$ such that if $m>\hat N$, then the $T_{k+1}$ coordinate of the solution started from $x_m$ is closer to $\bar{T}_{k+1}$ than $\hat \varepsilon$. Thus, the "limit trajectory" obtained above will entirely lie in the hyperplane $T_{k+1}=\bar{T}_{k+1}$. This means that we have found a homoclinic orbit in this hyperplane. However, on this hyperplane, our current $(k+2)$-dimensional system coincides with the $(k+1)$-dimensional system, for which global asymptotic stability of the equilibrium follows from the induction assumption. This excludes the presence of a homoclinic orbit and from this contradiction we obtain the global asymptotic stability of the equilibrium of the $(k+2)$-dimensional system.

Trivially, for $k=1$, the assertion holds, so repeating the inductive step leads to $\Gamma_0\in\mathcal{B}(E)$.

Example 2.1. As an example, let us assume that after performing the procedure described above, we obtain a sequence of reproduction numbers for which the relations $R_0^{(1)}\leq 1, \ R_0^{(2)}>1, \ R_0^{(3)} \leq 1, \ldots, \ R_0^{(n)}>1$ hold. This means that an equilibrium of the form $(S^*,0,T_2^*,0,\dots,T_n^*)$ is globally asymptotically stable on $\Gamma_0$. Following the above procedure we can obtain the equilibria which attract solutions started from $\Gamma\setminus\Gamma_0$ in each different case.

---

3. Application to three strains

To make the process for determining the global stability properties of the system (6) described in the previous section better visible and understandable, we consider the case $n=3$. For an even simpler case we refer the reader to [6] which corresponds to the case $n=2$.

Let us now turn to the case of three strains. To make our notations easier to follow and unambiguous, for the reproduction numbers and equilibria we will use the signs '$+$' and '$-$' in the upper indices, where adding a '$+$' sign denotes that the last reproduction number determined during the procedure was greater than 1 and adding a '$-$' sign means that the last reproduction number was less than or equal to 1. Also, to simplify the notations, we will use simple indices for the infection rates.

In the case of three strains, our system takes the form

$\begin{split} \frac{\mathrm{d}S(t)}{\mathrm{d}t}={}& B-bS(t)-S(t)\sum\limits_{k=1}^3\beta_k T_k(t)+\sum\limits_{k=1}^3\theta_k T_k(t),\\ \frac{\mathrm{d}T_1(t)}{\mathrm{d}t}={}&\beta_1S(t)T_1(t)-\beta_2T_1(t)T_2(t)-\beta_3T_1(t)T_3(t)-(b+\theta_1)T_1(t),\\ \frac{\mathrm{d}T_2(t)}{\mathrm{d}t}={}&\beta_2S(t)T_2(t)+\beta_2T_1(t)T_2(t)-\beta_3T_2(t)T_3(t)-(b+\theta_2)T_2(t),\\ \frac{\mathrm{d}T_3(t)}{\mathrm{d}t}={}&\beta_3S(t)T_3(t)+\beta_3T_1(t)T_3(t)+\beta_3T_2(t)T_3(t)-(b+d_3+\theta_3)T_3(t).\\ \end{split}$

(16)

Following the steps in the previous section, by introducing the notation $N_k(t)=S(t)+\sum_{j=1}^3T_j(t),\ k=1,2,3$, and further, the notation $U_n(t)=B/b-N_n(t)$, we may transcribe the above equation into the form

$\label{3dU} \begin{split} \frac{\mathrm{d}T_1(t)}{\mathrm{d}t}={}&\beta_1S(t)T_1(t)-\beta_2T_1(t)T_2(t)-\beta_3T_1(t)T_3(t)-(b+\theta_1)T_1(t),\\ \frac{\mathrm{d}T_2(t)}{\mathrm{d}t}={}&\beta_2S(t)T_2(t)+\beta_2T_1(t)T_2(t)-\beta_3T_2(t)T_3(t)-(b+\theta_2)T_2(t),\\ \frac{\mathrm{d}T_3(t)}{\mathrm{d}t}={}&\beta_3\left(\frac{B}{b}-U_3(t)\right)T_3(t)-(b+d_3+\theta_3)T_3(t),\\ \frac{\mathrm{d}U_3(t)}{\mathrm{d}t}={}&d_3 T_3(t)-b\, U_3(t). \end{split}$

The last two equations can be decoupled from the others and we obtain the two-dimensional system

$\label{3dred} \begin{split} \frac{\mathrm{d}T_3(t)}{\mathrm{d}t}={}&\beta_3\left(\frac{B}{b}-U_3(t)\right)T_3(t)-(b+d_3+\theta_3)T_3(t),\\ \frac{\mathrm{d}U_3(t)}{\mathrm{d}t}={}&d_3 T_3(t)-b\, U_3(t), \end{split}$

which has two equilibria: the trivial equilibrium $(0,0)$ and the positive equilibrium

$(T_3^*,U_3^*)=\left(\frac{B\beta_{3}-(b+d_3+\theta_3)b}{\beta_{3}(b+d_3)},\frac{d_3(B\beta_{3}-(b+d_3+\theta_3)b)}{b\beta_{3}(b+d_3)}\right),$

which only exists if

$\mathcal{R}_0 :=\frac{B \beta_{3}}{b(b+d_3+\theta_3)}>1.$

We know from the previous section that $(0,0)$ is a globally asymptotically stable equilibrium of (3) if ${\mathcal{R}_0}\leq 1$ and the positive equilibrium is globally asymptotically stable if $\mathcal{R}_0> 1$. Following the steps described in the previous section, we can reduce the system in both cases to 3 dimensions.

We start with the case $\mathcal{R}_0\leq 1$. In this case, we obtain the three-dimensional system

$\label{red3d_1} \begin{split} \frac{\mathrm{d}T_1(t)}{\mathrm{d}t}={}&\beta_1S(t)T_1(t)-\beta_2T_1(t)T_2(t)-(b+\theta_1)T_1(t),\\ \frac{\mathrm{d}T_2(t)}{\mathrm{d}t}={}&\beta_2S(t)T_2(t)+\beta_2T_1(t)T_2(t)-(b+\theta_2)T_2(t),\\ \frac{\mathrm{d}U_2(t)}{\mathrm{d}t}={}&-b\, U_2(t), \end{split}$

i.e. the last two equations have the form

$\begin{split} \frac{\mathrm{d}T_2(t)}{\mathrm{d}t}={}&\beta_2\frac{B}{b}T_2(t)-\beta_2U_2(t)T_2(t)-\beta_2(T_2(t))^2-(b+\theta_2)T_2(t),\\ \frac{\mathrm{d}U_2(t)}{\mathrm{d}t}={}&-b \,U_2(t). \end{split}$

This two-dimensional system has the two equilibria $(0,0)$ and

$(T_2^-,0)=\left(\frac{B}{b}-\frac{b+\theta_2}{\beta_2},0\right).$

The second equilibrium only exists if

$\mathcal{R}_0^-:=\frac{B \beta_2}{b(b+\theta_2)}>1,$

and $(0,0)$ is globally asymptotically stable if $\mathcal{R}_0^-\leq 1$, while $(T_2^-,0)$ is globally asymptotically stable if $\mathcal{R}_0^-> 1$.

Let us again proceed with the first case: if $\mathcal{R}_0^-\leq 1$ we obtain the two-dimensional system

$\label{red2d_1} \begin{split} \frac{\mathrm{d}T_1(t)}{\mathrm{d}t}={}&\beta_1\frac{B}{b}T_1(t)-\beta_1U(t)T_1(t)-\beta_1(T_1(t))^2-(b+\theta_1)T_1(t),\\ \frac{\mathrm{d}U_1(t)}{\mathrm{d}t}={}&-b\, U_1(t), \end{split}$

which again has two equilibria: $(0,0)$ and

$(T_1^{--},0)=\left(\frac{B}{b}-\frac{b+\theta_1}{\beta_1},0\right).$

The second equilibrium only exists if

$\mathcal{R}_0^{--}:=\frac{B \beta_1}{b(b+\theta_1)}>1,$

and the trivial equilibrium is globally asymptotically stable if $\mathcal{R}_0^{--}\leq 1$, while $(T_1^{--},0)$ is globally asymptotically stable if $\mathcal{R}_0^{--}> 1$.

If $\mathcal{R}_0^->1$, then the equilibrium $(T_2^-,0)$ is globally asymptotically stable. Thus, we obtain the system

$\begin{split} \frac{\mathrm{d}T_1(t)}{\mathrm{d}t}={}&\beta_1\frac{B^{-+}}{b^{-+}}T_1(t)-\beta_1U_1(t)T_1(t)-\beta_1(T_1(t))^2-(b^{-+}+\theta_1)T_1(t),\\ \frac{\mathrm{d}U_1(t)}{\mathrm{d}t}={}&-b^{-+} U_1(t), \end{split}$

where $B^{-+}:= B+\theta_2T_2^-$ and $b^{-+}:= b+\beta_2T_2^-$. This system has the two equilibria $(0,0)$ and

$(T_1^{-+},0):=\left(\frac{B^{-+}}{b^{-+}}-\frac{b^{-+}+\theta_1}{\beta_1},0\right).$

The latter equilibrium only exists if

$\mathcal{R}_0^{-+}:=\frac{B^{-+}\beta_1}{b^{-+}(b^{-+}+\theta_1)}>1.$

If $\mathcal{R}_0^{-+}\leq 1$, then $(0,0)$ is globally asymptotically stable, while if $\mathcal{R}_0^{-+}>1$, then $(T_1^{-+},0)$ is globally asymptotically stable.

Now we turn to the case $\mathcal{R}_0>1$. In this case, all solutions tend to the positive equilibrium $(T_3^*,U_3^*)$. By substituting these values into the rest of the equations, introducing the notations $B^+:= B+\theta_3T_3^*$ and $b^+:= b+\beta_3T_3^*$, we obtain the three-dimensional system

$\begin{split} \frac{\mathrm{d}T_1(t)}{\mathrm{d}t}={}&\beta_1S(t)T_1(t)-\beta_2T_1(t)T_2(t)-(b^++\theta_1)T_1(t),\\ \frac{\mathrm{d}T_2(t)}{\mathrm{d}t}={}&\beta_2S(t)T_2(t)+\beta_2T_1(t)T_2(t)-(b^++\theta_2)T_2(t),\\ \frac{\mathrm{d}U_2(t)}{\mathrm{d}t}={}&-b^+ U_2(t), \end{split}$

(17)

and further, by decoupling the last two equations and rewriting them,

$\begin{split} \frac{\mathrm{d}T_2(t)}{\mathrm{d}t}={}&\beta_2\frac{B^+}{b^+}T_2(t)-\beta_2U_2(t)T_2(t)-\beta_2(T_2(t))^2-(b^++\theta_2)T_2(t),\\ \frac{\mathrm{d}U_2(t)}{\mathrm{d}t}={}&-b^+ \,U_2(t). \end{split}$

This two-dimensional system has two equilibria, the trivial equilibrium $(0,0)$ and the positive equilibrium

$(T_2^+,0)=\left(\frac{B^+}{b^+}-\frac{(b^++\theta_2)}{\beta_{2}},0\right),$

which only exists if

$\mathcal{R}_0^+ :=\frac{B^+ \beta_{2}}{b^+(b^++\theta_2)}>1.$

The trivial equilibrium is globally asymptotically stable if $\mathcal{R}_0^+\leq 1$ and the equilibrium $(T_2^+,U_2^+)$ is globally asymptotically stable if $\mathcal{R}_0^+>1$.

Let us proceed with the case $\mathcal{R}_0^+\leq 1$. In this case, (17) can be reduced to the system

$\label{red3d_2} \begin{split} \frac{\mathrm{d}T_1(t)}{\mathrm{d}t}={}&\beta_1\frac{B^+}{b^+}T_1(t)-\beta_1U(t)T_1(t)-\beta_1(T_1(t))^2-(b^++\theta_1)T_1(t),\\ \frac{\mathrm{d}U_1(t)}{\mathrm{d}t}={}&-b^+ U_1(t). \end{split}$

This system has the two equilibria $(0,0)$ and

$(T_1^{+-},0)=\left(\frac{B^+}{b^+}-\frac{b^++\theta_1}{\beta_1}\right),$

with the second equilibrium only existing if

$\mathcal{R}_0^{+-}:= \frac{B^+ \beta_1}{b^+(b^++\theta_1)}>1.$

The equilibrium $(0,0)$ is globally asymptotically stable if $\mathcal{R}_0^{+-}\leq 1$, and the positive equilibrium is globally asymptotically stable for $\mathcal{R}_0^{+-}>1$.

In the case $\mathcal{R}_0^+>1$, we can reduce the system to

$\label{red2d_2} \begin{split} \frac{\mathrm{d}T_1(t)}{\mathrm{d}t}={}&\beta_1\frac{B^{++}}{b^{++}}T_1(t)-\beta_1T_1(t)U_1(t)-\beta_1(T_1(t))^2-(b^{++}+\theta_1)T_1(t)T_1(t),\\ \frac{\mathrm{d}U_1(t)}{\mathrm{d}t}={}&-b^{++}U_1(t), \end{split}$

where $B^{++}:= B^++\theta_2T_2^+$ and $b^{++}:= b^++\beta_2T_2^+$. Again, this system has two equilibria, $(0,0)$ and

$(T_1^{++},U_1^{++}) =\left(\frac{B^{++}}{b^{++}}-\frac{(b^{++}+\theta_1)}{\beta_{1}},0\right),$

with the latter one only existing if

$\mathcal{R}_0^{++}:=\frac{B^{++} \beta_{1}}{b^{++}(b^{++}+\theta_1)}>1.$

If $\mathcal{R}_0^{++}\leq 1$, then $(0,0)$ is globally asymptotically stable, while if $\mathcal{R}_0^{++}>1$, then $(T_1^{++},U_1^{++})$ is globally asymptotically stable.

Based on the above calculations and Theorem 2.2, we can formulate the following theorem on the global dynamics of the three-strain model (16). Similarly to the general case, we use the notation

$\Gamma_0=\{\,(S,T_1,T_2,T_3):S> 0,\ T_1>0,\ T_2>0,\ T_3>0\,\}.$

Theorem 3.1. The following statements hold for the stability of the equilibria of (16).

(ⅰ) If $\mathcal{R}_0\leq 1$, $\mathcal{R}_0^-\leq 1$ and $\mathcal{R}_0^{--}\leq 1$, then the equilibrium $\big(\tfrac{B}{b},0,0,0\big)$ is globally asymptotically stable on $\Gamma_0$.

(ⅱ) If $\mathcal{R}_0\leq 1$, $\mathcal{R}_0^-\leq 1$ and $\mathcal{R}_0^{--}> 1$, then the equilibrium $\big(\tfrac{B}{b}-T_1^{--},T_1^{--},0,0\big)$ is globally asymptotically stable on $\Gamma_0$.

(ⅲ) If $\mathcal{R}_0\leq1$, $\mathcal{R}_0^{-}>1$ and $\mathcal{R}_0^{-+}\leq 1$, then the equilibrium $\big(\tfrac{B^{-+}}{b^{-+}}-T_2^-,0,T_2^-,0\big)$ is globally asymptotically stable on $\Gamma_0$.

(ⅳ) If $\mathcal{R}_0\leq1$, $\mathcal{R}_0^->1$ and $\mathcal{R}_0^{-+}> 1$, then the equilibrium $\big(\tfrac{B^{-+}}{b^{-+}}-T_1^{-+}-T_2^-,T_1^{-+},T_2^-,0\big)$ is globally asymptotically stable on $\Gamma_0$.

(ⅴ) If $\mathcal{R}_0>1$, $\mathcal{R}_0^+\leq 1$ and $\mathcal{R}_0^{+-}\leq 1$, then $\big(\tfrac{B^+}{b^+}-T_3^*,0,0,T_3^*\big)$ is globally asymptotically stable on $\Gamma_0$.

(ⅵ) If $\mathcal{R}_0>1$, $\mathcal{R}_0^+\leq 1$ and $\mathcal{R}_0^{+-}>1$, then the equilibrium $\big(\frac{B^+}{b^+}-T_1^{+-}-T_3^*,T_1^{+-},0,T_3^*\big)$ is globally asymptotically stable on $\Gamma_0$.

(ⅶ) If $\mathcal{R}_0>1$, $\mathcal{R}_0^+> 1$ and $\mathcal{R}_0^{++}\leq 1$ then the equilibrium $\big(\frac{B^{++}}{b^{++}}-T_2^+-T_3^*,0,T_2^+,T_3^*\big)$ is globally asymptotically stable on $\Gamma_0$.

(ⅷ) If $\mathcal{R}_0>1$, $\mathcal{R}_0^+>1$ and $\mathcal{R}_0^{++}>1$, then the equilibrium $\big(\frac{B^{++}}{b^{++}}-U_1^{++}-T_1^{++}-T_2^+-T_3^*,T_1^{++},T_2^+,T_3^*\big)$ is globally asymptotically stable on $\Gamma_0$.

Theorem 3.1 gives a complete description for solutions started from $\Gamma_0$. For solutions started with initial values 0 for any of the infected compartments, we can refer to Proposition 1. However, as an example, we show the case $T_2(t)\equiv 0$.

In this case, the system (16) takes the form

$\begin{split} \frac{\mathrm{d}S(t)}{\mathrm{d}t}={}& B-bS(t)-\beta_1S(t)T_1(t)-\beta_3S(t)T_3(t)+\theta_1 T_1(t)+\theta_3 T_3(t),\\ \frac{\mathrm{d}T_1(t)}{\mathrm{d}t}={}&\beta_1S(t)T_1(t)-\beta_3T_1(t)T_3(t)-(b+\theta_1)T_1(t),\\ \frac{\mathrm{d}T_3(t)}{\mathrm{d}t}={}&\beta_3S(t)T_3(t)+\beta_3T_1(t)T_3(t)-(b+d_3+\theta_3)T_3(t).\\ \end{split}$

(18)

This system is also of the form (1), for $n=2$. For a complete description of the global dynamics of this system, we refer the reader to [6].

---

4. Discussion

We established an SIS model for a disease with multiple strains where a more infective strain can superinfect an individual infected by another strain. We developed a method to determine the global stability properties of the system. The main idea is that after some transformation of the model, two equations are decoupled from the others and the global stability of this two-dimensional system is completely described using the Dulac-Bendixson criterion and the Poincaré-Bendixson theorem. The dimension of our system can consequently be decreased by 1 substituting the values of the limit points of the two-dimensional system into the remaining equations, applying the theory of asymptotically autonomous differential equations. At each step, we derive a reproduction number that selects from the two possible equilibria of the two-dimensional system the globally asymptotically stable one. At the end of this procedure, we find the globally asymptotically stable equilibrium of the full system, where some of the strains coexist, depending on the sequence of the reproduction numbers.

We note that in the present paper only the most infective strain could be lethal. Without this assumption, in the presence of multiple lethal strains, the transformations (6)-(7) cannot be performed, so our procedure cannot be applied. Therefore, it remains an open question to investigate a similar model where all virus strains may be lethal. We conjecture that a similar global asymptotic stability result holds in that more general setting, too.

---

Acknowledgments

A. Dénes was supported by Hungarian Scientific Research Fund OTKA PD 112463. Y. Muroya was supported by Scientific Research (c), No. 24540219 of Japan Society for the Promotion of Science. G. Röst was supported by ERC Starting Grant Nr. 259559 and by Hungarian Scientific Research Fund OTKA K109782.

---

Dedicated to the memory of Professor Yoshiaki Muroya who passed away in 2015 after the submission of this paper