Many infectious diseases are comprised of multiple strains with examples including Influenza, tuberculosis, and Dengue virus.
The time evolution of such systems is linked to a complex landscape shaped by interactions between competing strains.  Possible long-term dynamics include the extinction of less competitive strains, convergence to multi-strain steady-states, or self-sustained oscillations.

This work considers a two-strain epidemic model in which the strains can interact indirectly via the immunity response generated following infections, and in which this immune response wanes with time.  In particular, we focus on scenarios where the rate of waning immunity is significantly faster than the rate of demographic turnover.   The first key result of this study is the explicit computation of the steady states of the nonlinear system of seven equations.   Following this result, we take advantage of the separation of time scales in the problem and use perturbation methods to analyze the stability of the fixed points.  In particular, we establish the conditions under which the system gives rise to the coexistence of the two strains and whether coexistence is attained via convergence to an endemic steady-state or via self-sustained oscillations.

Our study unveils two parameter regimes of distinct qualitative behavior of the system and characterizes the separatrix between them.  Within the first regime, the system gives rise to oscillatory coexistence for all feasible conditions.  In the second regime, the system's behavior is governed by a solution to a quadratic equation, potentially resulting in the convergence to a multi-strain endemic equilibrium or the persistence of oscillatory coexistence.

Keywords: Multiple strains, Multiple scales, Mathematical Epidemiology, Oscillations, Waning immunity, Cross-immunity, Coexistence

§ INTRODUCTION

Many infectious diseases such as seasonal influenza, are comprised of multiple strains with examples including seasonal influenza, Haemophilus influenza, Streptococcus pneumonia, human immunodeficiency virus (VIH), tuberculosis, and Dengue virus, see [13, 19] and references within.
These strains can interact indirectly via the immunity response generated following infections.  For example, in the case of cross-immunity, the immune response generated by one infection may provide some degree of protection against subsequent infections with closely related strains.  Alternatively, infection by one strain may enhance the ability of other strains to establish an infection leading to enhanced susceptibility.
The dynamics of multi-strain epidemic systems may lead to the gradual extinction of less competitive strains until only one dominant strain remains.  Indeed, following the exclusion principle, strains that are more transmissible or virulent will tend to outcompete other strains in a race of infecting susceptible individuals, leading to the extinction of less competitive strains over time.  However, multi-strain epidemic systems may also give rise to endemic steady-states in which multiple strains co-exist, or to self-sustained oscillations of multiple strains.  A common mechanism for self-sustained oscillations is that the interaction between the strains leads to a sequential wave of epidemics each dominated by a different strain.  Then, by the time the last epidemic wave of the sequence ends, there will have been sufficient births to increase susceptibility to the initial strain.  Examples of such systems have been presented and studied with as few as three interacting strains [1, 7, 4, 11].  Self-sustained oscillations or chaotic behavior had also been studied in systems with interacting  strains [3, 6, 8], and when additional effects such as quarantine or age structure are taken into account [15, 14, 16, 9].

Clearly, the persistence of a strain or of multiple strains in an epidemic system with recovery is possible only if the system includes a mechanism for replenishing the group of susceptible.  In all the works mentioned above, the sole mechanism considered was demographic turnover.
In the case of oscillations, one implication is that the oscillation period depends on the birth and death rates.  For example, the works [15, 14] showed that the period time is proportional to $1/\sqrt{\mu}$ where $\mu$ is the rate of both birth and death and presented simulations in which the period of oscillations was roughly 5 years.  While it is plausible that demographic turnover does impact epidemic dynamics over such a timescale, it is also of interest to explore other effects that may impact epidemic dynamics and induce oscillations with shorter periods, e.g., on the timescale of a season or a year.   Waning immunity is an example of such an alternative mechanism.  Immunity loss occurs, e.g., after infections by one of the strains of the human respiratory syncytial virus [18].  Furthermore, recent works show that waning immunity may better explain dengue dynamics [12].
The rate of waning immunity, $\tau$, can be 10-100 times faster than the rate of demographic turnover, i.e., $\tau\gg \mu$.   Therefore, whenever present, one can expect the effect of waning immunity to be dominant over the effect of birth and death.

In this work, we consider systems with two interacting strains for which

* We consider a general interaction between the strains.  Specifically, we assume that infection from one strain may provide either partial immunity (cross-immunity) or enhanced susceptibility to infection by the other strain.

* The mechanism for replenishing the group of susceptible is not demographic turnover, but rather waning immunity or loss of convalescent immunity over time.

We study this system with the aim of understanding its long-term dynamics.  Specifically, we ask under what conditions multiple strains persist and what is the nature of persistence. For example, whether the system converges to a coexistence steady-state solution or oscillates about it.

§.§ Paper outline and summary of results

In Section <ref> we extend the two-strain model [2, 13] to account for waning immunity and for a more general interaction between the strains. The key parameters of the model are the basic reproduction numbers (invasion numbers) $\mathcal{R}_1$, $\mathcal{R}_2$ for strains $1$ and $2$, respectively, the rate of waning immunity $\tau$, and the relative susceptibility $
\sigma_{ij}$ to strain $j$ for an individual previously infected with and recovered from strain $i$ ($i\ne j$).  Here $\sigma_{ij} = 0$ corresponds to total cross-immunity, $0 < \sigma_{ij} < 1$ corresponds to reduced susceptibility
(partial cross-immunity) and $\sigma_{ij} > 1$ corresponds to enhanced susceptibility.

In Section <ref> we compute the equilibrium points of the system.
The key result of this section is an explicit expression for the endemic multi-strain equilibrium point, aka the coexistence steady-state.   This is a surprising result since the coexistence steady-state is defined by a nonlinear algebraic system of seven equations. Hence, one can rarely expect to find an explicit solution for such a system.
This result opens the way to the analysis of the systems' behavior.  In particular, it enables studying the linear stability of the coexistence steady-state by obtaining an explicit expression for the characteristic polynomial of the relevant Jacobian matrix.  Yet, the characteristic polynomial gives rise to a septic equation that cannot be solved analytically in the general case.

In Section <ref>, we overcome the above difficulty by taking advantage of the separation of times scales in the problem and using perturbation methods to obtain an approximation of the characteristic polynomial roots when the rate of waning immunity is significantly smaller than the rate of recovery.  The arising Jacobian matrix, however, is asymmetric and suffers from a multiplicity of eigenvalues at leading orders.  Therefore, it does not lend itself easily to standard matrix perturbation methods [10].  We overcome these difficulties by application of perturbation theory to approximate the roots of the characteristic polynomial.
Using the above results, we map the systems' behavior as a function of its parameters, see Figure <ref>.  In particular, regions IV and V correspond to regions in which both strains coexist:  In Region IV the system tends toward an endemic multi-strain equilibrium, see, e.g., Figure <ref>$(A_3)$, while in region V the two strains are maintained via sustained oscillations, see, e.g., Figure <ref>$(A_2)$, or possibly by chaotic behavior.

Top graphs: Bifurcation diagram of system (<ref>) for the case $\tau=10^{-4}$, A) $\sigma_{12}=3$, $\sigma_{21}=2.5$ B) $\sigma_{12}=1$, $\sigma_{21}=1.5$ and C) $\sigma_{12}=0.5$, $\sigma_{21}=0.75$.  For each case, the $(\mathcal{R}_1,\mathcal{R}_2)$ plane is divided into as much as five subregions: Region I in which the system tends to a disease-free equilibrium, namely both strains become extinct.  Region II (III) in which only strain 1 (strain 2) will be maintained, see, e.g., bottom graph $(A_1)$.  Region IV in which the system dynamics give rise to an endemic multi-strain equilibrium, see, e.g., bottom graph $(A_3)$ and finally region V in which the two strains are maintained via sustained oscillations or chaotic behavior, see, e.g., bottom graph $(A_2)$. Bottom graphs: Solution $I_i(t)+J_i(t)$ of (<ref>) for $i=1,2$ corresponding to the points $(A_1), (A_2)$ and $(A_3)$ in graph A for which $\sigma_{12}=3$, $\sigma_{21}=2.5$ $\mathcal{R}_2=2$, and $A_1$) $\mathcal{R}_1=0.15$, $A_2$) $\mathcal{R}_1=1.5$ $A_3$) $\mathcal{R}_1=0.3$. In Figure ($A_1$), $I_1+J_1\equiv0$ (dashed-blue curve).

One can observe a qualitative difference between the properties of regions $V$ in Figures <ref>A and <ref>B.  In Section <ref>, we aim to understand the source of the above qualitative difference and to find the transition surface.  To do so, we focus on the behavior of the system when the basic reproduction numbers are close to one, aka near the organizing center $\mathcal{R}_1=\mathcal{R}_2=1$.  Our analysis reveals two distinct qualitative behaviors of the model separated by the transition surface $s=0$,

\begin{equation*}
\end{equation*}

For positive $s$, the system
either oscillates about the co-existence steady-state or develops chaos for any pair of basic reproduction numbers that are sufficiently close to one and for which both of the endemic single-strain steady-states are unstable, see, e.g., Figure <ref>A.
In contrast, for negative $s$, the parameter regime in which the system oscillates about the co-existence steady-state or develops chaos is determined by a solution of a quadratic equation.  This regime is generally only a sub-region of the feasible surrounding of the organizing center, see, e.g., Figure <ref>B.  In other cases, region V does not exist at all, see, e.g,. Figure <ref>C.

In Section <ref> we consider model generalizations and conduct a numerical investigation of their effect.  In particular, we consider the effect of birth and death, and relax the assumption that recovery rates and waning rates are equal for both strains.  When the waning rate is significantly larger than the rate of birth and death, the impact of demographic turnover is shown to be small compared to the impact of waning immunity.  We find, however, that the system is sensitive to other changes such as an asymmetric change in the recovery rates, providing an example of the rich mathematical structure of the system.
Concluding remarks are given in Section <ref>.

§ THE MATHEMATICAL MODEL

We consider a model for the dynamics of two strains in which infection from one strain provides
1) full immunity to re-infection by the same strain for some time. 2) possibly reduced susceptibility  (cross-immunity) or enhanced susceptibility to infection by the other strain.

The mathematical model is based on the model presented in [2, 13], and extended to account for waning immunity.  The population is divided into eight different compartments: susceptibles (S), those infected with strain i ($I_i$, primary infection), those recovered from strain i ($R_i$, as a result of primary infection), those infected with strain i after they had recovered from strain $j\ne i$ ($J_i$, secondary infection) and those recovered from both strains ($R_3$). The population is assumed to mix randomly.  The model is given by

\begin{equation}   \begin{split}
\frac{\text{d}S}{\text{d}t} &=\mu-\sum_{i=1}^{2}\beta_i(I_i+J_i)S-\mu S+\sum_{k=1}^{3}\tau_k  R_k,\\[1ex]
\frac{\text{d}I_1}{\text{d}t}  &= \beta_1S(I_1+J_1)-(\mu+\gamma_1)I_1,\\[1ex]
\frac{\text{d}I_2}{\text{d}t}  &= \beta_2S(I_2+J_2)-(\mu+\gamma_2)I_2,\\[1ex]
\frac{\text{d}J_1}{\text{d}t}  &= \beta_1\sigma_{21}R_2(I_1+J_1)-(\mu+\gamma_1)J_1,\\[1ex]
\frac{\text{d}J_2}{\text{d}t}  &= \beta_2\sigma_{12}R_1(I_2+J_2)-(\mu+\gamma_2)J_2,\\[1ex]
\frac{\text{d}{R_1}}{\text{d}t} &= \gamma_1I_1-\beta_2\sigma_{12}(I_2+J_2)R_1-(\mu+\tau_1)R_1,\\[1ex]
\frac{\text{d}{R_2}}{\text{d}t} &= \gamma_2I_2-\beta_1\sigma_{21}(I_1+J_1)R_2-(\mu+\tau_2)R_2,\\[1ex]
\end{split}
\end{equation}

were $\mu$ is the rate at which individuals are born, as well as the mortality rate, $\beta_i$ denotes the transmission coefficient for strain $i$, $\gamma_i$ denotes the recovery rate from strain $i$, $\tau_i$ is the rate at which immunity to re-infection by strain $i$ wanes and $\tau_3$ is the rate at which immunity to re-infection by both strains wanes.  Finally, $\sigma_{ij}$ is the relative susceptibility to strain $j$ for an individual previously infected with and recovered from strain $i$ ($i\ne j$), so that $\sigma_{ij}=0$ corresponds to total cross-immunity, $0<\sigma_{ij}<1$ corresponds to reduced susceptibility (partial cross-immunity) and $\sigma_{ij}>1$ corresponds to enhanced susceptibility.

We consider initial conditions which satisfy

\begin{equation}
\begin{split}
&S(0)\ge0, I_i(0)\ge0, J_i(0)\ge0, R_i(0)\ge0,\quad i=1,2,\\
&S(0)+I_1(0)+I_2(0)+J_1(0)+J_2(0)+R_1(0)+R_2(0)\le 1,
\end{split}
\end{equation}

and, thus, ensure that for all $t>0$,

\begin{equation}\label{eq:boundOnVariables}
\begin{split}
&S(t)\ge 0, I_i(t)\ge0, J_i(t)\ge0, R_i(t)\ge0,\quad i=1,2,\\
&S(t)+I_1(t)+I_2(t)+J_1(t)+J_2(t)+R_1(t)+R_2(t)+R_3(t)\equiv 1.
\end{split}
\end{equation}

We set $\gamma_1=1$.  This is equivalent to scaling time so that a single time unit equals to the recovery period after infection from strain $1$.  With respect to this scaling, we consider the case $\mu\ll \tau_i \ll 1$ for $i=1,2,3$.  Namely, the case for which the rate of birth and death is much slower than the rates of waning immunity $\tau_i$, which themselves are much slower than the recovery rate.  For example, after a recovery period of a week, immunity is gained for a characteristic period of one year and demographic turnover occurs over a scale of 80 years.  Accordingly, in what follows we neglect the effect of demographic turnover $\mu=0$.  Furthermore, for simplicity, our analysis will consider the case of symmetric rates

\[
\gamma_1=\gamma_2=1,\quad \tau_1=\tau_2=\tau_3=\tau,
\]

so that the model of focus takes the form

\begin{align}
\frac{\text{d}S}{\text{d}t} &=-\sum_{i=1}^{2}\mathcal{R}_i(I_i+J_i)S+\tau\sum_{k=1}^{3}  R_k,\label{eq:eqS}\\[1ex]
\frac{\text{d}I_1}{\text{d}t}  &= \mathcal{R}_1(I_1+J_1)S-I_1,\quad  \frac{\text{d}I_2}{\text{d}t}  = \mathcal{R}_2(I_2+J_2)S-I_2,\label{eq:eqIi}\\[1ex]
\frac{\text{d}J_1}{\text{d}t}  &= \mathcal{R}_2\sigma_{21}(I_1+J_1)R_2-J_1,\quad          \frac{\text{d}J_2}{\text{d}t}  = \mathcal{R}_2\sigma_{12}(I_2+J_2)R_1-J_2\label{eq:eqJi},\\[1ex]
\frac{d {R_1}}{d t} &= I_1-\mathcal{R}_2\sigma_{12}(I_2+J_2)R_1-\tau R_1,\quad
\frac{d {R_2}}{d t} = I_2-\mathcal{R}_1\sigma_{21}(I_1+J_1)R_2-\tau R_2\label{eq:eqRi},\\[1ex]
\end{align}

where $\mathcal{R}_i=\beta_i$ is the basic reproduction number of strain $i$, see Section <ref>.

We note that in Section <ref>, we will consider the impact of the above assumptions by studying numerically (<ref>) and comparing it to (<ref>).

§ SYSTEM EQUILIBRIA POINTS -  EXPLICIT EXPRESSIONS.

In what follows we will show that the system (<ref>) can have up to four equilibrium point:

* The disease-free equilibrium (DFE), $\phi^{DFE}$ for which

\[
S^{DFE}=1,\quad R^{DFE}_3=0,\quad I^{DFE}_i=J^{DFE}_i=R^{DFE}_i=0,\qquad i=1,2,
\]

2,3. Two single-strain endemic equilibrium (EE) points, $\phi^{EE,i}$, $i=1,2$, in which strain $i$ persists while strain $j\ne i$ does not,

\[
I^{EE,i}_i>0,\quad I^{EE,i}_j=0,\quad i=1,2,\quad j\ne i,
\]

4. A multi-strain endemic equilibrium, $\phi^{CE}$, in which both strains co-exist (CE),

\[
I^{CE}_1+J^{CE}_1>0,\quad I^{CE}_2+J^{CE}_2>0,
\]

aka the coexistence steady-state.

A standard computation of the next generation matrix [13, 5, 17] yields that the basic reproduction number is given by

\begin{equation}\label{eq:R0}
\mathcal{R}_0=\max{\{\mathcal{R}_1,\mathcal{R}_2\}},
\end{equation}

By computation, (<ref>) determines the stability of $\phi^{DFE}$ [13],

The disease-free equilibrium, $\phi^{DFE}$, is linearly stable when $\mathcal{R}_0<1$ and is an unstable (saddle) whenever $\mathcal{R}_0>1$.

In what follows, we find explicit expressions for the additional three equilibrium points and consider their stability.

§.§ The single-strain endemic steady states and their stability

The endemic steady-states $\phi^{EE,i}$, $i=1,2$, are  given by

\begin{equation}\label{eq:phiEE}
\begin{split}
S^{EE,i}&=\frac1{\mathcal{R}_i},\quad I_i^{EE,i}=\frac{\tau}{1+\tau}\frac{\mathcal{R}_i-1}{\mathcal{R}_i},\quad
J_i^{EE,i}&=J_j^{EE,i}=I_j^{EE,i}=R_j^{EE,i}=R_3^{EE,i}=0,\quad j\ne i.
\end{split}
\end{equation}

The endemic steady state $\phi^{EE,i}$ exists only when $\mathcal{R}_i>1$.

The endemic steady state $\phi^{EE,i}$ exists when $I_i^{EE,i}>0$.  This implies that $\mathcal{R}_i>1$, see (<ref>).

The single-strain endemic steady state $\phi^{EE,i}$ (<ref>)  is linearly stable if

\begin{equation}\label{eq:condEE}
\mathcal{R}_j<\frac{(1+\tau)\mathcal{R}_i}{1+\tau+\sigma_{ij}(\mathcal{R}_i-1)},\quad j\ne i,
\end{equation}

See Appendix <ref>.

Regions II and III in Figure <ref> correspond to the stability regions of $\phi^{EE,2}$ and $\phi^{EE,1}$, respectively.  Lemma <ref> and Proposition <ref> explicitly define the boundaries of these regions.  The following Lemma uses these results to show that regions II and III do not overlap.

It is impossible for the single-species endemic states (<ref>), $\phi^{EE,1}$ and $\phi^{EE,2}$, to exist simultaneously and be stable.

Assume, without loss of generality that $\phi^{EE,1}$ exists and is stable.  Then, by Lemma <ref>, $\mathcal{R}_1>1$.  Thus,  condition (<ref>) implies that $\mathcal{R}_2<1$.    In this case, however, Lemma <ref> implies that $\phi^2_{EE}$ does not exist.  Similarly, if $\phi^{EE,2}$ exists and is stable, then $\phi^{EE,1}$ does not exist.

§.§ The coexistence steady state $\phi^{CE}$

The following proposition characterizes the steady state $\phi^{CE}$ of coexistence, namely a fixed point for which both $I_1>0$ and $I_2>0$.  The steady-state consists of seven variables $(S^*,I_1^*,I^*_2,J^*_1,J^*_2,R^*_1,R^*_2)$ that are a solution of a nonlinear algebraic system of seven equations.  One can rarely expect to find an explicit solution for such a system, unless in carefully selected cases.
Nevertheless, Proposition <ref> provides an explicit expression for $\phi^{CE}$, when exists, for any given set of parameters,

There exists a steady-state solution

\[
\phi^{CE}=(S^*,I_1^*,I^*_2,J^*_1,J^*_2,R^*_1,R^*_2)
\]

of (<ref>) in the parameter regime

\begin{equation}\label{eq:parm_regime_CE}
\mathcal{R}_1>\frac{(1+\tau)\mathcal{R}_2}{1+\tau+\sigma_{21}(\mathcal{R}_2-1)},\quad \mathcal{R}_2>\frac{(1+\tau)\mathcal{R}_1}{1+\tau+\sigma_{12}(\mathcal{R}_1-1)}
\end{equation}

which satisfies $I_1^*>0$ and $I_2^*>0$, and which is given explicitly by

\begin{equation}\label{eq:S*I*}
S^*=\frac{-b-\sqrt{b^2-4ac}}{2a},\quad I_1^*=a_1+b_1S^*,\quad I_2^*=a_2+b_2S^*,
\end{equation}

\begin{equation}\label{eq:S*I*_coeffs}
\begin{split}
a&=\mathcal{R}_1\mathcal{R}_2(\sigma_{21} - \sigma_{12})\left[(\mathcal{R}_1+ \mathcal{R}_2)(\sigma_{21}+ \sigma_{12})+\tau (\mathcal{R}_1-\mathcal{R}_2)\right],\\
b&=\mathcal{R}_1\mathcal{R}_2(\sigma_{12}+\sigma_{21})-\mathcal{R}_1\mathcal{R}_2\sigma_{12}\sigma_{21}(1+\mathcal{R}_1 + \mathcal{R}_2)-\tau(\sigma_{21}\mathcal{R}_1^2+\sigma_{12}\mathcal{R}_2^2),\\
\frac{(2 \tau + 1)\sigma_{12} \sigma_{21}a_1}{\mathcal{R}_1}&= \frac{ (\sigma_{21} - \sigma_{12})(\sigma_{12}^2 \sigma_{21} + \sigma_{12} \sigma_{21}^2 - \sigma_{12} \sigma_{21} + \sigma_{21}^2 \tau^2 - \sigma_{12}^2 \tau^2)}{(\mathcal{R}_1+\mathcal{R}_2) \sigma_{12} \sigma_{21} - \mathcal{R}_1 \sigma_{12} - \mathcal{R}_2 \sigma_{21} - (\mathcal{R}_1-\mathcal{R}_2) (\sigma_{12}-\sigma_{21}) \tau},\\
\frac{b_1}{2\mathcal{R}_1\mathcal{R}_2}&=\frac{{\mathcal{R}_2\sigma_{12}\tau - \mathcal{R}_2\sigma_{21}\tau + \mathcal{R}_1\sigma_{21} - \mathcal{R}_1\sigma_{12} + \mathcal{R}_1\sigma_{12}^2\tau - \mathcal{R}_1\sigma_{21}^2\tau - \mathcal{R}_2\sigma_{12}^2 + \mathcal{R}_2\sigma_{21}^2}}
{{\mathcal{R}_1\sigma_{12}\sigma_{21}  + \mathcal{R}_2\sigma_{12}\sigma_{21}  - \mathcal{R}_1\sigma_{12} - \mathcal{R}_2\sigma_{21}}-\tau(\mathcal{R}_1-\mathcal{R}_2) (\sigma_{12}-\sigma_{21})},\\
\frac{(2 \tau + 1)\sigma_{12} \sigma_{21}a_2}{\mathcal{R}_2} &=\frac{ (\sigma_{12} - \sigma_{21})(\sigma_{21}^2 \sigma_{12} + \sigma_{21} \sigma_{12}^2 - \sigma_{12} \sigma_{21} + \sigma_{12}^2 \tau^2 - \sigma_{21}^2 \tau^2)}{(\mathcal{R}_1+\mathcal{R}_2) \sigma_{12} \sigma_{21} - \mathcal{R}_1 \sigma_{12} - \mathcal{R}_2 \sigma_{21} - (\mathcal{R}_2-\mathcal{R}_1) (\sigma_{21}-\sigma_{12}) \tau},\\
\frac{b_2}{2\mathcal{R}_1\mathcal{R}_2}&=\frac{{\mathcal{R}_1\sigma_{21}\tau - \mathcal{R}_1\sigma_{12}\tau + \mathcal{R}_2\sigma_{12} - \mathcal{R}_2\sigma_{21} + \mathcal{R}_2\sigma_{21}^2\tau - \mathcal{R}_2\sigma_{12}^2\tau - \mathcal{R}_1\sigma_{21}^2 + \mathcal{R}_1\sigma_{12}^2}}
{{\mathcal{R}_1\sigma_{12}\sigma_{21}  + \mathcal{R}_2\sigma_{12}\sigma_{21}  - \mathcal{R}_1\sigma_{12} - \mathcal{R}_2\sigma_{21}}-\tau(\mathcal{R}_1-\mathcal{R}_2) (\sigma_{12}-\sigma_{21})},
\end{split}
\end{equation}

\begin{equation}\label{eq:Ri}
R_2^*=\frac{1-\mathcal{R}_1S^*}{\mathcal{R}_1\sigma_{21}},\quad R^*_1=\frac{1-\mathcal{R}_2S^*}{\mathcal{R}_2\sigma_{12}},
\end{equation}

\begin{equation}\label{eq:Ji}
J_1^*=I_2^*-\tau R^*_2,\quad J^*_2=I_1^*-\tau R_1^*.
\end{equation}

The steady-state solution satisfies $J_1^\prime(t)+R_2^\prime(t)=0=J_2^\prime(t)+R_1^\prime(t)$.  Thus, Equations (<ref>) and (<ref>) imply (<ref>).

We next derive relation (<ref>): Equations (<ref>) imply

\[
\frac{I_1^*}{I^*_1+J^*_1}=\mathcal{R}_1S^*,\quad \frac{I^*_2}{I^*_2+J^*_2}=\mathcal{R}_2S^*.
\]

\[
\frac{J_1^*}{I_1^*+J_1^*}=1-\mathcal{R}_1S^*,\quad \frac{J^*_2}{I^*_2+J^*_2}=1-\mathcal{R}_2S^*.
\]

Substituting the above expression in (<ref>) defines $R^*_i$ as a linear function of $S^*$

\[
\mathcal{R}_1\sigma_{21}R_2^*=\frac{J_1^*}{I_1^*+J^*_1}=1-\mathcal{R}_1S^*,\quad \mathcal{R}_2\sigma_{12}R^*_1=\frac{J_2^*}{I_2^*+J_2^*}=1-\mathcal{R}_2S^*.
\]

Substituting (<ref>) and (<ref>) into (<ref>) yields a nonlinear system of three equations $S^\prime(t)=I^\prime_1(t)=I^\prime_2(t)=0$ for three unknowns: $S^*$, $I_1^*$ and $I_2^*$. Solving this system gives rise to (<ref>) and (<ref>).

Finally, the requirement that $I_1^*>0$ and $I_2^*>0$ leads to the inequalities (<ref>).

Proposition <ref> explicitly defines the coexistence steady-state solution $\phi^{CE}$.  We note that the regime in which $\phi^{CE}$ is defined coincides with the regime in which both single-species endemic steady states $\{I_{EE,i}\}_{i=1,2}$ are unstable, see (<ref>) and (<ref>).  This implies that

\[
\mathcal{R}_1>1,\quad \mathcal{R}_2>\frac{(1+\tau)\mathcal{R}_1}{1+\tau+\sigma_{12}(\mathcal{R}_1-1)},\quad\mbox{or}\quad \mathcal{R}_2>1,\quad \mathcal{R}_1>\frac{(1+\tau)\mathcal{R}_2}{1+\tau+\sigma_{21}(\mathcal{R}_2-1)},
\]

then the solutions of the system (<ref>) may either converge to $\phi^{CE}$, converge to a limit cycle that oscillates about it, or develops chaos.  In all cases, the two species will persist.

To go beyond Lemma <ref> and determine whether the solutions of the system (<ref>) will converge to $\phi^{CE}$ or not, we next study the stability of $\phi^{CE}$.

The explicit expression of $\phi^{CE}$ as defined in Proposition <ref> suggests that one can determine the stability of $\phi^{CE}$ by considering the Jacobian matrix at $\phi^{CE}=(S^*,I_1^*,I_2^*,J_1^*,J_2^*,R_1^*,R_2^*)$, and studying the behaviour of the corresponding characteristic polynomial.  Although it is straightforward to obtain an explicit expression for the corresponding characteristic polynomial in this way, computing the eigenvalues requires solving a septic equation.  It is impossible to do so analytically in a general case.  To overcome this difficulty, in the next section, we take advantage of the separation of times scales in the problem and find an approximation of the eigenvalues when $\tau\ll1$.

§ COEXISTENCE - A PERTURBATIVE APPROACH

In what follows, we utilize the separation of time scales, $\tau\ll1$, and use a perturbative approach to study the stability of $\phi^{CE}$.
Using Proposition <ref>, we compute the Jacobian of (<ref>) at $\phi^{CE}$ and its corresponding characteristic polynomial.  Collecting orders of $\tau$ gives rise to the following form of characteristic polynomial
\begin{equation}\label{eq:P(lambda)}
P(\lambda)=\lambda^5(\lambda+1)^2+\tau \lambda^3 P_1(\lambda)+\tau^2 \lambda P_2(\lambda)+\tau^3 P_3(\lambda)+O(\tau^4),
\end{equation}

where $\qquad P_i(0)\ne 0,\quad i=1,2,3.$

Let $0<\tau\ll1$.
Then, the roots $\{\lambda_i\}_{i=1}^7$ of the polynomial (<ref>) are given by

\begin{equation}\label{eq:lambda12}
\lambda_{1,2}=-1+O(\tau).
\end{equation}

\begin{equation}\label{eq:lambda3}
\lambda_3=-\tau\frac{P_3(0)}{P_2(0)}+O(\tau\sqrt{\tau}),
\end{equation}

and for $j=4,5,\cdots,7,$

λ_j=A_j√(τ)+B_j τ+O(τ√(τ)),     $|P_1^2(0)-4P_2(0)|\gg \tau$,

λ_j= a_j√(τ)+b_j τ^3/4+c_j τ+O(τ√(τ)),     $P_1^2(0)-4P_2(0)=O(\sqrt\tau)$,

where $\{A_j,B_j\}_{j=4}^7$ are defined by

\begin{equation}\label{eq:lambda47_alphabeta}
\begin{split}
A_j^2&=-\frac{P_{1}(0)\pm \sqrt{P_1^2(0) - 4P_2(0)}}2,\\
[P^\prime_2(0)-2P_2(0)]A_j^2 -[P^\prime_1(0)-2P_1(0)][A_j^2P_1(0)+P_2(0)] +P_3(0)}{ 2A_j^2P_1(0)+4P_2(0)},\quad
\end{split}
\end{equation}

and $\{a_j,b_j,c_j\}_{j=4}^7$ are defined by

\begin{equation}\label{eq:lambda47_ab}
\begin{split}
b_j&=\pm\frac12\sqrt{\frac{-P_1^3(0) + P_1^2(0)P^\prime_1(0) - 2P_1(0)P^\prime_2(0) + 2\sqrt{-2P_1(0)}\Delta P_2 + 4P_3(0)}{\sqrt{-2}P_1^{3/2}(0)}}
\end{split}
\end{equation}

\begin{equation}\label{eq:lambda47_c}
c_j=\frac{3(-2P_1(0))^{5/2} + 4(-2P_1(0))^{3/2}P^\prime_1(0) - 32b_i^2P_1(0) + 8\sqrt{-2P_1(0)}P^\prime_2(0) + 8\Delta P_2}{16\sqrt{-2P_1(0)}P_1(0)},
\end{equation}

\[
\Delta P_2=\lim_{\tau\to0^+}\frac{P_2(0)-\frac14P_1^2(0)}{\sqrt\tau}.
\]

For $\tau\ll1$, $P(\lambda)$ has two roots near $-1$ that satisfy (<ref>), and five roots near $\lambda=0$. When $\lambda\ll1$, dominant balance show that the dominant terms of $P(\lambda)$ are either $\lambda^5,\tau \lambda^3 P_1(\lambda)$ and $\tau^2 \lambda P_2(\lambda)$ or $\tau \lambda^3 P_1(\lambda)$ and $\tau^2 \lambda P_2(\lambda)$.  The latter case gives rise to a solution that satisfies (<ref>).  In the former case, $\lambda=O(\sqrt\tau)$.  Substituting (<ref>) in (<ref>) and equating orders yields (<ref>).  However, approximation (<ref>,<ref>) is not uniformly valid.  Indeed, when

\[
\]

then $B_j$ is undefined, see (<ref>), since it's denominator equals zero.
Approximation (<ref>) is valid only when

\[
B_j\tau\ll A_j \sqrt{\tau}.
\]

This implies the validity region $|P_1^2(0)-4P_2(0)|\gg\tau$.

To approximate the roots in the case $P_1^2(0)-4P_2(0)\ll1$, we consider

\[
P_2(0)=\frac14 P_1^2(0)+\sqrt\tau\Delta P_2.
\]

Substituting (<ref>) into (<ref>) and equating orders yields (<ref>,<ref>).

Proposition <ref> provides an approximation to the roots $\{\lambda_i\}_{i=1}^7$ of the relevant characteristic polynomial.
The next subsection is devoted to the implications of Proposition <ref>.
A systematic numerical verification of Proposition <ref> is presented subsequently in Section <ref>.

Proposition <ref> presents an approximation to the solution given the values of all parameters including the value of $\tau$.  The results are valid at the limit $\tau\to0^+$.   For example, Equation (<ref>) implies that

\[
\lim_{\tau\to0^+}\frac{\lambda_{1}+1}{\tau}= C,
\]

where $C$ is a (finite) constant.  In this case, the result does not depend on the behavior of the other parameters in the limit $\tau\to0^+$.  In contrast, the approximation of $\{\lambda_i\}_{i=4}^7$ at the limit of $\tau\to0^+$ does depend on the behavior of the other parameters as $\tau\to0^+$.  For example,  condition $P_1^2(0)-4P_2(0)=O(\sqrt\tau)$ in (<ref>) corresponds to a case $P_1^2(0)-4P_2(0)=0$ or alternatively to a case when some parameters concurrently change with $\tau$ so that $P_1^2(0;\mathcal{R}_1,\mathcal{R}_2,\sigma_{12},\sigma_{21})-4P_2(0;\mathcal{R}_1,\mathcal{R}_2,\sigma_{12},\sigma_{21})$ remains sufficiently small as $\tau\to0^+$.

§.§ Implications of Proposition <ref>

Proposition <ref> enables finding the stability region of the coexistence steady-state $\phi^{CE}$ by using the approximation to the roots $\{\lambda_i\}_{i=1}^7$ to test when $\mathrm{Re}(\lambda_i)<0$ for $i=1,2,\cdots,7$.  This approximation is explicit, where the explicit expressions for

\[
P_i(0)=P_i(0;\mathcal{R}_1,\mathcal{R}_2,\sigma_{12},\sigma_{21}),\quad P^\prime _j(0)=P_j^\prime(0;\mathcal{R}_1,\mathcal{R}_2,\sigma_{12},\sigma_{21}),
\]

where $i=1,2,3,\quad j=1,2$,
that appear in (<ref>) are presented in the supplementary material.
Nevertheless, these explicit expressions are long and cumbersome.  Therefore using Proposition <ref> to analytically study the stability of $\phi^{CE}$ is a notorious task.  In this subsection, we will use Proposition <ref> to study the stability of $\phi^{CE}$ numerically.

We first test the sign of $\{\lambda_i\}_{i=3}^7$ in a wide range of cases, including the cases considered in Figure <ref> and in all subsequent examples in the manuscript (data not shown).  In all cases we test numerically the expressions provided in the supplementary material
and find that

\[
P_2(0)>0,\quad P_3(0)>0, \quad P_1^2(0)-4P_2(0)>0.
\]

This implies that, to leading order, $\lambda_3<0$, see (<ref>).  Therefore, the stability of $\phi^{CE}$ is determined by $\{\lambda_i\}_{i=4}^7$.  The above findings also imply that $A_j$ and $a_j$ are pure imaginary number, and therefore that the sign of $\mathrm{Re}\lambda_i$ for $i=4,5,\cdots,7$ is determined by $B_j$ and $b_j$, respectively, see Equations (<ref>). This is the reason approximations (<ref>) go beyond the leading order.

We next use Proposition <ref> and Equations (<ref>) to compute the region in the $(\mathcal{R}_1, \mathcal{R}_2)$ plane in which the co-existence steady state, $\phi^{CE}$, of (<ref>) exists but is unstable.  In this case, all steady-state solutions of (<ref>) are unstable and hence the system either approaches a periodic solution or develops chaotic behavior.  We note that this region corresponds to region V in terms of Figure <ref>.

The region in the $(\mathcal{R}_1, \mathcal{R}_2)$ plane in which $\phi^{CE}$ is unstable is marked by blue in Figure <ref>.  In comparison, the region of instability approximated using Proposition <ref> is marked by a shaded region with dashed red boundaries.

Bifurcation diagram of $\phi^{CE}$ in the $(\mathcal{R}_1, \mathcal{R}_2)$ plane.  The coexistence steady-state $\phi^{CE}$ is a stable steady-state solution of (<ref>) in the green regions.  In the blue region, the coexistence steady-state exists but it is unstable.  The shaded region with dashed red boundaries is an approximation of the blue region computed using (<ref>).  All cases correspond to $\tau=10^{-3}$, $\sigma_{12}=2.5$ and A: $\sigma_{21}=0.5$,  B: $\sigma_{21}=1$, C: $\sigma_{21}=2$, and D: $\sigma_{21}=2$.

We observe that the approximated region well agrees with the numerically computed region of instability only for moderate values of $\mathcal{R}_1$ and $\mathcal{R}_2$, see Figure <ref>.  This can be expected, since Proposition <ref> relies on perturbation analysis that is valid for sufficiently small $\tau$.  For example, the derivation of (<ref>) relies on the identification of the relevant dominant terms in (<ref>), where, e.g., $\tau^3P_3(\lambda)$ is not one of the dominant terms since $\tau^3 P_3(\lambda)\ll\tau \lambda^3 P_1(\lambda)$. These condition implies that derivation (<ref>) is valid only when

\begin{equation}\label{eq:assumption_O1}
\sqrt{\tau}P_3(0)\ll P_1(0).
\end{equation}

In practice, however, the coefficients $P_i(0)$ and $P_i^\prime(0)$ involve high powers of $\mathcal{R}_1, \mathcal{R}_2$, $\sigma_{12}$ and $\sigma_{21}$, see supplementary material,
and can therefore be of varying magnitude.  For example, for $\mathcal{R}_1=\mathcal{R}_2=5$ and $\sigma_{12}=\sigma_{21}=2$,

\[
P_1(0)\approx 12.4,\quad P_3(0)\approx 642,
\]

so that for this choice of parameters (<ref>) is not satisfied for $\tau=10^{-3}$.  As expected, the approximation improves for smaller values of $\tau$, see, e.g., Figure <ref>.

Same as Figure <ref> for $\tau=5\cdot 10^{-4}$.

The above examples demonstrate how Proposition <ref> and Equations (<ref>) can be used to map the behavior of system (<ref>) in different parameter regions.  These examples combine the asymptotic results of Proposition <ref> with numerical methods.  Yet, this approach is limited.  For example, while one can observe a qualitative difference between the stability regions in Figures <ref>A, <ref>B and Figures <ref>C, <ref>D, confirming the observed phenomenon occurs generically and finding a transition surface by analyzing the cumbersome explicit expressions of Proposition <ref> is a demanding task.
To overcome this problem, in Section <ref> we further restrict ourselves to the behavior of the system when $\mathcal{R}_1\approx \mathcal{R}_2\approx 1$.

§.§ Numerical verification of Proposition <ref>

Proposition <ref> is valid for sufficiently small $\tau>0$.  In what follows, we verify that the results of Proposition <ref> agree with numerical computations for small but finite values of $\tau$, and in particular for characteristic values of $\tau$ in an epidemic setting.

Figure <ref>A
presents the approximation errors

\begin{equation}\label{eq:Ei}
E_i=|\lambda_i-\lambda_i^{\rm approx}|,
\end{equation}

where $\lambda_i^{\rm approx}$ are is the approximation (<ref>) to the root $\lambda_i$ of (<ref>) for the case $\sigma_{12}=0.6$, $\sigma_{21}=3$, $\mathcal{R}_1=2$ and $\mathcal{R}_2=1.5$.  As expected, the error $E_{1,2}=\max_{i=1,2}E_i=O(\tau)$, see (<ref>), while $E_{3-7}=\max_{i=3}^7E_i=O(\tau\sqrt\tau)$, see (<ref>) and (<ref>).

Graph of approximation error $E_{1,2}=\max\{E_1,E_2\}$ ($\circ$) and $E_{3-7}=\max\{E_3,E_4,\cdots,E_7\}$ ($\square$) as a function of $\tau$ where $E_i$ is given by (<ref>) for the cases: A) $\mathcal{R}_1 = 2$, $\mathcal{R}_2 = 1.5$, $\sigma_{21}  = 0.6$ and $\sigma_{12}  = 3$.  Super-imposed are the curves $0.7\tau$ (dashes) and $2.3\,\tau\sqrt{\tau}$ (solid).  B) $\mathcal{R}_1 = 2$, $\mathcal{R}_2 = 1.8$, $\sigma_{21}  = 0.6$ and $\sigma_{12}  = 3$. Super-imposed are the curves $0.8\tau$ (dashes) and $3\,\tau\sqrt{\tau}$ (solid). C) $\mathcal{R}_1 = 3$, $\mathcal{R}_2 = 1.8$, $\sigma_{21}  = 2$ and $\sigma_{12}  = 3$. Super-imposed are the curves $2.1\tau$ (dashes) and $4.5\,\tau\sqrt{\tau}$ (solid).

Graph of $\mathrm{Re}(\lambda_7)$ (solid blue) for the case $\tau=10^{-3}$, $\mathcal{R}_1=2$, $\sigma_{12}=3$, and $\sigma_{21}=0.6$.
Super-imposed is approximation (<ref>) (dash-dotted black) and approximation (<ref>) (dashed red).

The proof of Proposition <ref> shows that approximation (<ref>) is not uniformly valid.  Figure <ref> demonstrates this point by presenting the eigenvalue approximation (<ref>) near a point $\mathcal{R}_2=\mathcal{R}_2^*\approx 1.78$ in which $P_1^2(0;\mathcal{R}_i,\sigma_{ij})=4P_2(0;\mathcal{R}_i,\sigma_{ij})$ and the denominator of $B_j$ equals zero, see (<ref>).  As expected, we observe that the eigenvalue
approximation (<ref>) for $\lambda_7$ (dash-dotted curve) well approximates the eigenvalue (solid blue curve) for $\mathcal{R}_2$ significantly far from $\mathcal{R}_2^*$.  Yet, for $\mathcal{R}_2\approx\mathcal{R}_2^*$, the approximation (<ref>) fails.  Rather, as expected, (<ref>) well approximates the eigenvalue for $\mathcal{R}_2\approx\mathcal{R}_2^*$.

Figure <ref> also demonstrates why is it important to approximate the eigenvalues near the point $\mathcal{R}_2^*$ where the approximation (<ref>) fails.  Indeed, we observe in Figure <ref> that $Re(\lambda_i)=0$, and hence the linear stability of $\phi^{CE}$ changes, near the point $\mathcal{R}_2^*$ in which $P_1^2(0;\mathcal{R}_i,\sigma_{ij})=4P_2(0;\mathcal{R}_i,\sigma_{ij})$.

§ COEXISTENCE AND OSCILLATIONS FOR REPRODUCTION NUMBERS ABOUT ONE

In this section, we focus our analysis on the behavior of the system (<ref>) at reproduction numbers about one, $\mathcal{R}_1\approx \mathcal{R}_2\approx 1$. In particular, we aim to understand the source of the qualitative difference between the stability regions in Figures <ref>A, <ref>B and Figures <ref>C, <ref>D, and to find the transition surface between them.  The first step of applying Proposition <ref> to reproduction numbers about one already guides to answers to these questions.

Let $0<\tau\ll1$, $\mathcal{R}_1=1+b_1\tau^2$ and $\mathcal{R}_2=1+b_2\tau^2$ in the domain (<ref>), and define

\begin{equation}\label{eq:s_def}
\end{equation}

Then, for $\tau\ll1$, the steady-state solution (<ref>)

\[
\phi^{CE}=(S^*,I_1^*,I^*_2,J^*_1,J^*_2,R^*_1,R^*_2)
\]

of (<ref>) takes the form

\begin{equation}\label{eq:S*_betanearOne}
S^*=\frac{2\sigma_{12}\sigma_{21}+s-\sqrt{s^2-2(\sigma_{12}+\sigma_{21})\tau s+c_3\tau^2}+(\sigma_{12}+\sigma_{21})\tau-c_2\tau^2}{2\sigma_{12}\sigma_{21}+2s+c_1\tau^2}+O(\tau^3),
\end{equation}

\begin{equation*}
\begin{split}
c_3&=[(4s +2(\sigma_{12}+\sigma_{21})^2+2(\sigma_{12}-\sigma_{21}) ](b_1+b_2) + (\sigma_{21} + \sigma_{12})^2 - 4b_2\sigma_{12}\sigma_{21}(\sigma_{12} - \sigma_{21}),
\end{split}
\end{equation*}

\begin{equation}\label{eq:Ii_betanearOne}
\begin{split}
\end{split}
\end{equation}

\begin{equation}\label{eq:Ri_betanearOne}
\begin{split}
\end{split}
\end{equation}

\begin{equation}\label{eq:Ji_betanearOne}
J_1^*=I_2^*-\tau R^*_2,\quad J^*_2=I_1^*-\tau R_1^*.
\end{equation}

The results are obtained by computation the Taylor series expansion of (<ref>) about $\tau=0$, but without implicitly assuming that $\tau \ll |s|$.

Equations (<ref>-<ref>) of Proposition <ref> imply that the values $I_i^*$, $J_i^*$ and $R^*_i$ of the steady-state solution $\phi^{CE}$ are, to leading order, proportional to $1-S^*$.  Equation (<ref>) shows that, to leading order, $1-S^*$ depends on the result of $s-\sqrt{s^2+O(\tau)}$ or $s-|s|$.  In particular, $1-S^*$ changes its behavior when $s$ changes its sign.

Let $|s|\gg \tau$, $0<\tau\ll1$, where $s$ is defined by (<ref>).  Then, when $s$ is positive, (<ref>) satisfies for $s\gg\tau$

\begin{equation}\label{eq:S*_s_positive}
\frac{\sigma_{12}^2\sigma_{21}^2(b_1b_2\sigma_{21}\sigma_{12}-b_1\sigma_{21}-b_2\sigma_{12})}{s(\sigma_{12}\sigma_{21}+s)^2}\tau^2+O(\tau^3),
\end{equation}

and when $s$ is negative,

\begin{equation}\label{eq:S*_s_negative}
S^*=1+\frac{b_1\sigma_{12}+b_2\sigma_{12}}{s}\tau^2+O(\tau^3),\quad s\ll-\tau.
\end{equation}

Lemma <ref> show that, indeed, $1-S^*$ changes its qualitative behavior when $s$ changes its sign.
For positive $s$, $1-S^*=O(1)$ and is independent of the exact value of $\mathcal{R}_1$ and $\mathcal{R}_2$ as given by $b_1,b_2$, see (<ref>).
In contrast, for negative $s$, $1-S^*=O(\tau^2)$ and depends, to leading order, on $b_1,b_2$, see (<ref>).

Figure <ref>A presents $1-S^*$ as a function of $s=\sigma_{12}\sigma_{12}-\sigma_{12}-\sigma_{21}$ (solid blue curve).  As expected, $1-S^*$ changes its qualitative behavior as $s$ changes its sign.  The approximations (<ref>) and (<ref>) arising from Proposition <ref> and Lemma (<ref>), respectively, are super-imposed and agree well with the graph of $1-S^*$.

A: Graph of $S^*$ as a function of $s$ (solid blue curve) for the case $\sigma_{12}=2.5$,  $\mathcal{R}_1=1+\tau$, $\mathcal{R}_2=1+2\tau$ and $\tau=0.01$.  Note that $\sigma_{21}=(s+\sigma_{12})/(\sigma_{12}-1)$.  Super-imposed are the approximations given by (<ref>) (dashed black curve), and by (<ref>) (dash-dotted red curve).  B: Same data for $I_1^*$.  Note that (<ref>) is plotted outside the surrounding of $s=0$ to avoid the removable singularity.  The three curves are indistinguishable in each graph, except partially near $s=0$.

Figure <ref>B presents $I_1^*$ as a function of $s$.  As expected, see (<ref>), $I_1^*$ behaves qualitatively similar to $1-S^*$.

The above results focused on the value $\phi^{CE}$ as a function of $s$.  The following proposition goes further and studies the behavior of $\phi^{CE}$ as a function of $s$.  In particular, it considers the linear stability of $\phi^{CE}$ by approximating the eigenvalues of the Jacobian of (<ref>) at $\phi^{CE}$.

\[
\]

Then, for sufficiently small $\tau$, and for reproduction numbers sufficiently close to one

\[
\mathcal{R}_1=1+b_1\tau^2,\quad \mathcal{R}_2=1+b_2\tau^2,\qquad 0<\tau\ll1,
\]

the roots $\{\lambda_i\}_{i=1}^7$ of the characteristic polynomial $P(\lambda)$ of the Jacobian matrix of (<ref>) at $\phi^{CE}$
are given by (<ref>) for $s>0$ and by (<ref>) for $s<0$:

\begin{equation}\label{eq:lambda12_beta_around_one_positive}
\lambda_j= \pm\sqrt{s\tau}i\frac{\sqrt{(\sigma_{21} + \sigma_{12})(1+s)}}{\sigma_{12}\sigma_{21}+s}+O(\tau),\quad j=1,2,\quad s>0,
\end{equation}

\begin{equation}\label{eq:lambda35_beta_around_one_positive}
\lambda_j=\alpha_j\frac{1}{\sqrt[3]{(\sigma_{21} + \sigma_{12})\sigma_{12}\sigma_{21}}}(s\tau)^{\frac23}+O(\tau),\quad \alpha_j=\sqrt[3]{-1}=-1,\frac12\pm\frac{\sqrt{3}}2i,\quad j=3,4,5,
\end{equation}

\begin{equation}\label{eq:lambda67_beta_around_one_positive}
\lambda_{6,7}=-1+O(\tau).
\end{equation}

For $s<0$, three roots $\{\lambda_i\}_{i=1}^3$ of the characteristic polynomial $P(\lambda)$ satisfy

\begin{equation}\label{eq:lambda13_beta_around_one_minus}
\lambda_j= -\tau+O(\tau^2),\quad j=1,2,\cdots,5,\quad s<0,
\end{equation}
two additional roots $\{\lambda_i\}_{i=4}^5$ satisfy

\begin{equation}\label{eq:lambda45_beta_around_one_minus}
\lambda_j=x_j\tau^2 +O(\tau^3),
\end{equation}

where $\{x_j\}_{j=4}^5$ are the roots of the quadratic equation

\begin{equation}\label{eq:lambda45_beta_around_one_xj_minus}
x_j^2-\frac{b_1\sigma_{12}+b_2\sigma_{21}}{s}x_j -\frac{((\sigma_{12}-1)b_1+b_2)((\sigma_{21}-1)b_2+b_1)}{s}=0,
\end{equation}

and the last two roots satisfy

\begin{equation}\label{eq:lambda67_beta_around_one_minus}
\lambda_{6,7}=-1+O(\tau^3).
\end{equation}

When $s>0$, then for sufficiently small $\tau$ so that Lemma <ref> applies, the characteristic polynomial  takes the form

\begin{equation}\label{eq:P(lambda)_beta_around_one}
P(\lambda)=\lambda^5(\lambda+1)^2+\tau \lambda^3 P_1(\lambda)+\tau^2 \lambda^2 P_2(\lambda)+\tau^3 P_3(\lambda)+O(\tau^4),\quad P_i(0)\ne0,\quad i=1,2,3.
\end{equation}

$P(\lambda)$ has two roots near $-1$ that satisfy (<ref>), and five roots near $\lambda=0$. When $\lambda\ll1$, dominant balance show that the dominant terms of $P(\lambda)$ are either $\lambda^5$ and $\tau \lambda^3 P_1(\lambda)$ or only $\tau \lambda^3 P_1(\lambda)$ and $\tau^3 P_3(\lambda)$.  In the former case, $\lambda=O(\sqrt\tau)$ and satisfies to leading order a quadratic equation.  Further perturbation analysis yields that for $j=1,2$

\begin{equation*}
\lambda_j= \pm\sqrt{\tau}\sqrt{-P_1(0)}+\frac{2P_1^3(0) - P_1^2(0)P^\prime_1(0) + P_1(0)P_2(0) - P_3(0)}{2P_1^2(0)}\tau+O(\tau\sqrt\tau),
\end{equation*}

giving rise to (<ref>).

The latter case implies

\begin{equation*}
\lambda_j=\sqrt[3]{-\frac{P_3(0)}{P_1(0)}}\tau^{\frac23}-\frac{P_1(0)P_2(0) - P_3(0)}{3P_1^2(0)}\tau+O(\tau\sqrt\tau),\quad j=3,4,5,
\end{equation*}

yielding (<ref>).

When $s<0$, then for sufficiently small $\tau$ so that Lemma <ref> applies, the characteristic polynomial takes the form

\begin{equation}\label{eq:P(lambda)_beta_around_one_negative}
\begin{split}
\lambda^4 +3\tau^2 \lambda^3 +\tau^3 \lambda^2 P_3(\lambda)+\tau^4 \lambda^2 P_4(\lambda)+\tau^5\lambda P_5(\lambda)\right]+\\&\tau^6\lambda P_6(\lambda)+\tau^7 P_7(\lambda),
\end{split}
\end{equation}

where $P_i(0)\ne 0$ for $i=1,2,\dots,7.$
$P(\lambda)$ has two roots near $-1$ that satisfy (<ref>), and five roots near $\lambda=0$.
When $\lambda\ll1$, dominant balance shows that either $\lambda=O(\tau)$ or $\lambda=O(\tau^2)$.  In the former case,
three roots satisfy

\[
\lambda_j=x_j\tau +O(\tau^2),\quad j=1,2,3,
\]

Where $x_j$ are the roots of the cubic polynomial

\[
x_j^3+x_j^2P_1(0) +x_jP_2(0) +P_3(0)=(x_j+1)^3=0.
\]

This polynomial has a triple root $x_j=-1$, yielding (<ref>).
In the latter case, the two additional roots satisfy $\lambda_j=x_j\tau^2 +O(\tau^3)$ where

\[
x_j^2 P_3(0)+x_jP_5(0) +P_7(0),\quad j=4,5,
\]

yielding (<ref>) and (<ref>).

In the following, we first discuss the implications of Proposition <ref>.  Systematic numerical verification of Proposition <ref> is presented subsequently in Section <ref>, and its validity region is considered in <ref>.

Proposition <ref> gives rise to the following picture:

* For positive $s$ and sufficiently small $\tau$, the coexistence steady-state $\phi^{CE}$ is linearly unstable since, for example, ${\rm Re}\lambda_4>0$.  This result in independent of $b_1$ and $b_2$, thus $\phi^{CE}$ is unstable in the whole region of the parameter plane $(\mathcal{R}_1,\mathcal{R}_2)$ which is close to $(1,1)$ and within (<ref>).  In this region, the system (<ref>) either oscillates about the co-existence steady-state or develops chaos.  This is in agreement with the observed behavior in Figures <ref>C and <ref>D corresponding to $s=0.5$ and $s=1.5$, respectively, as well as to the behavior in Figure <ref>A ($s=2$).

* For negative $s$, the stability of $\phi^{CE}$ is determined by the sign of the real part of the roots of the quadratic equation (<ref>).  Equation (<ref>) depends on $b_1,b_2$ and $s$, and therefore, one can expect that in general $\phi^{CE}$ will be unstable (if at all) only in a sub-region of the parameter plane $(\mathcal{R}_1,\mathcal{R}_2)$ which is close to $(1,1)$. See, e.g, Figures <ref>A and <ref>B corresponding to $s=-1.75$ and $s=-1$, respectively, as well as Figures <ref>B ($s=-1$) and <ref>C ($s=-0.875$).

* For negative $s$ and when the interaction between the strains is characterized solely by cross-immunity, $0<\sigma_{12},\sigma_{21}<1$, we will show subsequently in Lemma <ref> that $\phi^{CE}$ is linearly stable in a surrounding of $(1,1)$ in the parameter plane $(\mathcal{R}_1,\mathcal{R}_2)$.  Hence, in this case, the system (<ref>) will not give rise to oscillatory behavior.
See, for example, Figure <ref>C with corresponding $s=-0.875$.

In particular, Proposition <ref> reveals the source of the qualitative difference between the stability regions in Figures <ref>A, <ref>B and Figures <ref>C, <ref>D, and to find the transition surface between them.

Under the conditions of Proposition <ref>, one can also deduce that when the interaction between the strains is characterized solely by cross-immunity, $0<\sigma_{12},\sigma_{21}<1$, the system (<ref>) will not give rise to oscillatory behavior.

Let $0<\sigma_{12},\sigma_{21}<1$.  Then, for sufficiently small $\tau$ and for reproduction numbers sufficiently close to one

\begin{equation}\label{eq:Rnearone}
\mathcal{R}_1=1+b_1\tau^2,\quad \mathcal{R}_2=1+b_2\tau^2,\quad 0<\tau\ll1,
\end{equation}

the steady-state $\phi^{CE}$ of (<ref>), if exists, is linearly stable.

Proposition <ref> implies that $\phi^{CE}$ exists in the parameter regime (<ref>) or in terms of (<ref>),

\begin{equation}\label{eq:parm_regime_CE_bi}
\end{equation}

The case $0<\sigma_{12},\sigma_{21}<1$ corresponds to $s<0$.
By Proposition <ref>, the roots of the characteristic polynomial of the Jacobian matrix of (<ref>) at $\phi^{CE}$ are approximated by (<ref>).

The only roots that may have a positive real part are $\{\lambda_j\}_{j=4}^5$, given by (<ref>) and the quadratic equation (<ref>).  Equation (<ref>) has one branch of negative solutions for any $b_1,b_2>0$.  The second branch of solutions
vanishes in the boundaries of (<ref>), and negative in the regime (<ref>).  Namely, for any $0<\sigma_{12},\sigma_{21}<1$, $x_j=0$ when

\[
\frac{b_2}{b_1}=1-\sigma_{12},\quad\mbox{or}\quad \frac{b_1}{b_2}=1-\sigma_{21},
\]

and $x_j<0$ when (<ref>).  Thus, for sufficiently small $\tau$, $\mathrm{Re}\lambda_i<0$ for $i=1,2,\cdots,7$ implying that $\phi^{CE}$ is linearly stable.

§.§ Numerical verification of Proposition <ref>

A) Case $s<0$: Graph of approximation error $E_{1-3}=\max\{E_1,E_2,E_3\}$ ($\circ$) and $E_{4-7}=\max\{E_4,\cdots,E_7\}$ ($\square$) as a function of $\tau$ where $E_i$ is given by (<ref>) and  $\lambda_i^{\rm approx}$ is given by (<ref>) for the case $\mathcal{R}_1 = 1+1.5\tau^2$, $\mathcal{R}_2 = 1+0.5\tau^2$, $\sigma_{21}  = 1$ and $\sigma_{12}  = 3$ corresponding to $s=-1$. Super-imposed are the curves $3.5\tau^2$ (dashes) and $4.15\tau^3$ (solid).  B) Case $s>0$: Graph of approximation error $E_{1-7}=\max\{E_1,E_2,\cdots,E_7\}$ ($\circ$) where $E_i$ is given by (<ref>) and  $\lambda_i^{\rm approx}$ is given by (<ref>)
for the case $\mathcal{R}_1 = 1+1.5\tau^2$, $\mathcal{R}_2 = 1+0.5\tau^2$, $\sigma_{21}  = 2$ and $\sigma_{12} = 3$ corresponding to $s=1$.   Super-imposed is the curve $0.9\tau$ (dashes).

Figure <ref>
presents the approximation error of (<ref>) and (<ref>) as a function of $\tau$ for cases with positive and negative values of $s$.  As expected, in both cases, the numerical error is in full agreement with the results of Proposition <ref>.

§.§ Validity region of approximations (<ref>) and (<ref>)

Proposition <ref> shows that given a value of $s$, there exists $0<\tau\ll1$ such for which approximations (<ref>) and (<ref>) are valid.    Section <ref> shows that the results of Proposition <ref> are in full agreement with the numerical results.  We now consider a different question: Given a value of $0<\tau\ll1$, what is the range of values of $s$ for which approximations (<ref>) and (<ref>) are valid.

Lemma <ref> relies on the condition $|s|\gg \tau$ since, only in this regime, the Taylor series expansion of the expression for $S^*$ about $\tau=0$ is valid.
It is not clear, however, whether the series expansions used to derive approximation (<ref>) are also valid in the regime $|s|\gg \tau$.
The correction term (<ref>) of reads as

\begin{equation*}
\begin{split}
\lambda_j&=\pm\sqrt{\tau}\sqrt{-P_1(0)}+\frac{2P_1^3(0) - P_1^2(0)P^\prime_1(0) + P_1(0)P_2(0) - P_3(0)}{2P_1^2(0)}\tau+O(\tau\sqrt{\tau})\\
&=\pm \sqrt{s\tau}\frac{\sqrt{(\sigma_{21}+\sigma_{21})(1+s)}}{s+\sigma_{12}\sigma_{21}}-\left[\frac{1}2+O(s)\right]\tau+O(\tau\sqrt{\tau}).
\end{split}
\end{equation*}

This approximation is valid only when the correction term is significantly smaller than the leading order term, $\sqrt{s\tau}\gg \tau$, implying the condition $s\gg\tau$.  Yet, this argument does not ensure that the next order correction term is also significantly smaller than the correction term when $s\gg\tau$.  Indeed, our analysis does not imply any restriction on the behavior of the series coefficient as $s\to0$.  Therefore, more restrictive conditions may arise when considering higher-order terms or other expansions.  For example, the correction term (<ref>) of reads as

\begin{equation*}
\begin{split}
\lambda_3&=\sqrt[3]{-\frac{P_3(0)}{P_1(0)}}\tau^{\frac23}-\frac{P_1(0)P_2(0) - P_3(0)}{3P_1^2(0)}\tau+O(\tau\sqrt\tau)\\&=\sqrt[3]{\frac{s\tau}{(\sigma_{12}+\sigma_{21})\sigma_{12}\sigma{21}}}+\left(2+\frac{(\sigma_{21}-1)^2}{\sigma_{21}^2}s+O(s^2)\right)\frac{\tau}3+O(\tau\sqrt{\tau}).
\end{split}
\end{equation*}

Therefore, Approximation (<ref>) is valid only when  $(s\tau)^{2/3}\gg \tau$, implying the more restrictive condition $s\gg\sqrt{\tau}$.  Inspection of the next order term may in turn further restrict the region of validity.

We further attempt to determine numerically the validity region of approximations (<ref>) and (<ref>), e.g., for $j=4,5$.  To do so, we first plot the approximation error as a function of $s$, see Figure <ref>A.  As expected, we observe that the error increases sharply near $s=0$, and the interval of large error around $s=0$ decreases with $\tau$.  In Figure <ref>B we plot the location of the peak of the error as a function of $\tau$, and observe that the peak location behaves as $\sqrt{\tau}$.  Therefore, we conclude that the relevant approximations in the case considered are valid in the region $|s|\gg\sqrt{\tau}$.

A) Graph of approximation error $E_{4,5}$ as a function of $s$ where $E_i$ is given by (<ref>) and  $\lambda_i^{\rm approx}$ is given by (<ref>) for $s\le0$ and (<ref>) for $s>0$. For the cases $\tau=1e-4$
$\mathcal{R}_1 = 1+1.5\tau^2$, $\mathcal{R}_2 = 1+0.5\tau^2$, $\sigma_{12}  = 3$ corresponding to $s=-1$. Super-imposed are the curves $3.5\tau^2$ (dashes) and $4.15\tau^3$ (solid).  B) Case $s>0$: Graph of approximation error $E_{1-7}=\max\{E_1,E_2,\cdots,E_7\}$ ($\circ$) where $E_i$ is given by (<ref>) and  $\lambda_i^{\rm approx}$ is given by (<ref>)
for the case $\mathcal{R}_1 = 1+1.5\tau^2$, $\mathcal{R}_2 = 1+0.5\tau^2$, $\sigma_{21}  = 2$ and $\sigma_{12} = 3$ corresponding to $s=1$.   Super-imposed is the curve $0.9\tau$ (dashes).

§ NUMERICAL STUDY OF THE GENERAL MODEL (<REF>)

In Sections <ref>-<ref>, we focused on the model (<ref>).  This model is derived from the more general model (<ref>) by neglecting demographic turnover, $\mu=0$, and assuming that recovery rates and waning rates are symmetric $\gamma_1=\gamma_2$, and $\tau_1=\tau_2=\tau_3=\tau$.

In this section, we study numerically the general model (<ref>) in aim of testing the effect of relaxing the assumptions that give rise to the reduced model (<ref>).

We first consider the effect of demographic turnover, $\mu>0$.
The rate of waning immunity $\tau$ can be as much as 10-100 times faster than the rate of demographic turnover, e.g., for a waning time of one year versus a lifespan of 80 years.  Thus, plausible values of the ratio $\mu/\tau$ are in the range of $1/100$ to $1/10$.
In Figure <ref>A, we rerun the simulation presented in Figure <ref>$(A_2)$ with $\mu=\tau/10$, i.e., a value in the upper range of plausible values.  As expected, we observe that the effect of slow demographic turnover is perturbative.  In Figures <ref>B and <ref>C, we consider faster rates of demographic turnover with values of $\mu=\tau/5$ and $\mu=\tau/2$, respectively.  These relatively high ratios of $\mu/\tau$ lead to a change in amplitude and period of the oscillations but do not change the qualitative behavior of the solution.

Effect of $\mu$:
Solution of (<ref>) (solid blue) for the same parameters as in Figure <ref>$(A_2)$ expect A: $\mu=\tau/10$, B: $\mu=\tau/5$, C: $\mu=\tau/2$.  The solution of Figure <ref>$(A_2)$ is super-imposed in all graphs  (dotted black).

Effect of $\tau$:
Solution of (<ref>) (solid blue) for the same parameters as in Figure <ref>$(A_2)$ expect $tau_1=(1+\delta)\tau$ and $\tau_2=(1-\delta)\tau$ where A: $\delta=0.05$, B: $\delta=0.2$, C: $\delta=0.3$.  The solution of Figure <ref>$(A_2)$ is super-imposed in all graphs (dotted black).

In Figure <ref>, we relax the assumption $\tau_1=\tau_2=\tau_3$ and study the effect of an asymmetry in the waning rates.  To do so, we rerun the simulation presented in Figure <ref>$A_2$ with $\tau_1=(1-\delta)\tau$, $\tau_2=(1+\delta)\tau,$ and $\tau_3=\tau$ where $\delta=0.05,0.2$ and $0.3$. We observe that varying these parameters by  $5-30\%$ leads to a change in the amplitude and period of the oscillations but does not change the qualitative behavior of the solution.

Finally, In Figure <ref>, we relax the assumption $\gamma_2=\gamma_1=1$ and study the effect of an asymmetry in the recovery rate.  To do so, we rerun the simulation presented in Figure <ref>$A_2$ with $\gamma_2=1.01, 1.025$ and $1.05$, while maintaining $\gamma_1=1$.  We observe the results are sensitive to changes in the recovery rates.  Indeed, with already $5\%$ change, the oscillations undergo a period-doubling bifurcation.  At slightly larger values of $\gamma_2$, the oscillatory pattern breaks, and the system converges to a multi-strain endemic equilibrium (data not show).

Effect of $\gamma$: Solution of (<ref>) (solid blue) for the same parameters as in Figure <ref>$(A_2)$ expect where A: $\gamma_2=1.01$, B: $\gamma_2=1.025$, C: $\gamma_2=1.05$.  The solution of Figure <ref>$(A_2)$ is super-imposed in all graphs (dotted black).

§ CONCLUDING REMARKS

We have presented a study of a two-strain epidemic model in which the strains can interact indirectly via the immunity response generated following infections, and in which this immune response wanes with time.
In particular, we have provided explicit expression for all equilibrium states of the underlying system, including the endemic multi-strain equilibrium or co-existence steady state, see Section <ref>.  We have also provided an approximation to the parameter regime in which the co-existence steady state is stable, see Section <ref>.  Our results reveal two parameter regimes of distinct qualitative behavior of the system, and we characterize the transition surface between them, see Section <ref>.  In particular, our work sheds light on the question 'When does an epidemic model give rise to self-sustained oscillations?'.   Indeed, we show that in the model considered self-sustained oscillations arise when $s=\sigma_{12}\sigma_{21}-\sigma_{12}-\sigma_{21}>0$, and arise conditionally when $s$ is negative.  We also show that the model does not give rise to self-sustained oscillations when $0<\sigma_{21},\sigma_{12}<1$.

Our analysis focused on the case in which waning immunity, $\tau_i$, is dominant over demographic turnover, $\mu\ll\tau_i$.  We are aware that some of our results do not extend in the absence of waning immunity.  For example, while we showed
that the model does not give rise to self-sustained oscillations when $0<\sigma_{21},\sigma_{12}<1$, the work [3] considers the model (<ref>) in the absence of waning immunity, $\tau_i=0$, but for $\mu>0$, and provides an example of periodic solutions when $0<\sigma_{21},\sigma_{12}<1$, yet $\sigma_{12}\approx1$.

From a methodological point of view, the explicit computation of the coexistence steady-state presented in Proposition <ref> was a key result which opened the way to the analysis of the problem.  Such a result is unexpected since the coexistence steady-state is defined by a nonlinear algebraic system of seven equations, and to the best of our knowledge does not appear elsewhere in the literature considering similar systems.  The same approach can be used to obtain an explicit solution to the more general system (<ref>) with $\mu>0$.  Yet, for $\mu>0$, we find that the solution for $S^*$ is expressed in terms of a solution to a cubic equation, rather than a quadratic equation, see (<ref>), and the overall expressions are cumbersome.  Thus, in retrospect, neglecting demographic turnover by setting $\mu=0$ was key in the explicit computation of the coexistence steady-state and in applying this result in further analysis.
The generalization of Proposition <ref> to (<ref>) will be presented elsewhere.

Determining the stability of the coexistence steady-state requires finding the eigenvalues of the Jacobian matrix.  Although Proposition <ref> enables computing this matrix explicitly, one cannot compute analytically the eigenvalues of the underlying 7 by 7 matrix in the general case.  We overcome the above difficulty by taking advantage of the separation of times scales in the problem and using perturbation methods.  The Jacobian matrix, however, is asymmetric and suffers from a multiplicity of eigenvalues at leading orders.  Therefore, it does not lend itself easily to standard matrix perturbation methods [10].  This is evident, for example, from the $O(\tau)$ and $O(\sqrt{\tau})$ scale of the different eigenvalues of the matrix when $\mathcal{R}_i=O(1)$ (Section <ref>), and how this scale varies near $\mathcal{R}_i\approx1$ (Section <ref>).  To overcome this difficulty, in Propositions <ref> and <ref> we analyze directly the characteristic polynomial roots.

Our perturbative results, e.g., Proposition <ref>, are valid when $\mathcal{R}_i=O(1)$ as $\tau\to0$.  Such a region of validity can be expected when applying perturbation methods, and is of most relevance for studying common infectious diseases, see [13, 19] and references within.  Yet, there is also theoretical interest in understanding the behavior of the system in an ultimate limit in which one of the strains has a vast advantage over the other, $\mathcal{R}_1\ll \mathcal{R}_2$.  An asymptotic study in such a regime of large $\mathcal{R}_i$ will be presented elsewhere.

Our numerical results suggest that the generalized model (<ref>) with asymmetric rates between the strains and demographic turnover exhibits a rich mathematical behavior including additional Hopf bifurcations and period-doubling bifurcations.  Additional effects such as the impact of additional strains, quarantine, age structure, sub-populations, vaccination, and treatment are also expected to effect the systems' behavior.
A study including these generalizations will be presented in a subsequent work.

§ ACKNOWLEDGMENTS

We are grateful to Prof. Guy Katriel for his most useful comments.

[1]
V. Andreasen, J. Lin, and S. A. Levin, The dynamics of cocirculating
influenza strains conferring partial cross-immunity, Journal of mathematical
biology, 35 (1997), pp. 825–842.

[2]
C. Castillo-Chavez, H. W. Hethcote, V. Andreasen, S. A. Levin, and W. M.
Liu, Epidemiological models with age structure, proportionate mixing,
and cross-immunity, Journal of mathematical biology, 27 (1989),
pp. 233–258.

[3]
K. Chung and R. Lui, Dynamics of two-strain influenza model with
cross-immunity and no quarantine class, Journal of mathematical biology, 73
(2016), pp. 1467–1489.

[4]
J. Dawes and J. Gog, The onset of oscillatory dynamics in models of
multiple disease strains, Journal of mathematical biology, 45 (2002),
pp. 471–510.

[5]
O. Diekmann, J. A. P. Heesterbeek, and J. A. Metz, On the definition
and the computation of the basic reproduction ratio R0 in models for
infectious diseases in heterogeneous populations, Journal of mathematical
biology, 28 (1990), pp. 365–382.

[6]
N. Ferguson, R. Anderson, and S. Gupta, The effect of
antibody-dependent enhancement on the transmission dynamics and persistence
of multiple-strain pathogens, Proceedings of the National Academy of
Sciences, 96 (1999), pp. 790–794.

[7]
M. G. M. Gomes and G. F. Medley, Dynamics of multiple strains of
infectious agents coupled by cross-immunity: a comparison of models, in
Mathematical approaches for emerging and reemerging infectious diseases:
models, methods, and theory, Springer, 2002, pp. 171–191.

[8]
S. Gupta, N. Ferguson, and R. Anderson, Chaos, persistence, and
evolution of strain structure in antigenically diverse infectious agents,
Science, 280 (1998), pp. 912–915.

[9]
M. A. Kuddus, E. S. McBryde, A. I. Adekunle, and M. T. Meehan, Analysis and simulation of a two-strain disease model with nonlinear
incidence, Chaos, Solitons & Fractals, 155 (2022), p. 111637.

[10]
R.-C. Li, Matrix perturbation theory, in Handbook of linear
algebra, Chapman and Hall/CRC, 2006, pp. 15–1.

[11]
J. Lin, V. Andreasen, and S. A. Levin, Dynamics of influenza a
drift: the linear three-strain model, Mathematical biosciences, 162 (1999),
pp. 33–51.

[12]
L. López, R. E. Paul, V.-M. Cao-Lormeau, and X. Rodó, Considering waning immunity to better explain dengue dynamics, Epidemics, 41
(2022), p. 100630.

[13]
M. Martcheva, An introduction to mathematical epidemiology,
vol. 61, Springer, 2015.

[14]
M. Nuno, C. Castillo-Chavez, Z. Feng, and M. Martcheva, Mathematical
models of influenza: the role of cross-immunity, quarantine and
age-structure, Mathematical epidemiology,  (2008), pp. 349–364.

[15]
M. Nuño, Z. Feng, M. Martcheva, and C. Castillo-Chavez, Dynamics
of two-strain influenza with isolation and partial cross-immunity, SIAM
Journal on Applied Mathematics, 65 (2005), pp. 964–982.

[16]
H. R. Thieme, Pathogen competition and coexistence and the evolution
of virulence, Mathematics for Life Science and Medicine,  (2007),
pp. 123–153.

[17]
P. Van den Driessche and J. Watmough, Reproduction numbers and
sub-threshold endemic equilibria for compartmental models of disease
transmission, Mathematical biosciences, 180 (2002), pp. 29–48.

[18]
L. White, M. Waris, P. Cane, D. J. Nokes, and G. Medley, The
transmission dynamics of groups a and b human respiratory syncytial virus
(hrsv) in england & wales and finland: seasonality and cross-protection,
Epidemiology & Infection, 133 (2005), pp. 279–289.

[19]
G. P. Wormser and B. Pourbohloul, Modeling infectious diseases in
humans and animals by matthew james keeling and pejman rohani princeton, nj:
Princeton university press, 2008. 408 pp., illustrated (hardcover), 2008.

§ THE ENDEMIC STEADY STATES AND THEIR STABILITY

The endemic steady-states $\phi^{EE,i}$, $i=1,2$, are  given by (<ref>).

Using (<ref>) gives rise to an explicit expression for the Jacobian matrix of (<ref>) at $\phi^{EE,i}$, $i=1,2$ whose eigenvalues are given by

\begin{equation*}
\begin{split}
\lambda_1&=-\tau,\quad  \lambda_{2,3}=-1,\quad \lambda_4=-\tau-\frac{\tau}{\tau+1}(\mathcal{R}_i-1)\sigma_{ji},
\\\lambda_5&=\frac{[\sigma_{ij}(\mathcal{R}_i-1)+\tau+1]\mathcal{R}_j}{(\tau+1)\mathcal{R}_i}-1,
\\\lambda_{6,7}&=-\tau\frac{\mathcal{R}_i+\tau \pm\sqrt{\tau^2(\mathcal{R}_i+\tau)^2-4\tau(1+\tau)^2(\mathcal{R}_i-1)}}{2(\tau+1)}.
\end{split}
\end{equation*}

Non-negativeness of each of the variable values of $\phi^{EE,i}$ implies that the endemic steady-state exists only when $\mathcal{R}_i>1$, see (<ref>).  Therefore, all eigenvalues, except possibly $\lambda_5$, have a negative real part.  The condition $\lambda_5>0$ gives rise to (<ref>).