Revisiting the exclusion principle in epidemiology at its ultimate limit

Nir Gavish¹¹1Faculty of Mathematics, Technion, Haifa, 32000, Israel; [email protected]

Abstract

The competitive exclusion principle in epidemiology implies that when competing strains of a pathogen provide complete protection for each other, the strain with the largest reproduction number outcompetes the other strains and drives them to extinction. The introduction of various trade-off mechanisms may facilitate the coexistence of competing strains, especially when their respective basic reproduction numbers are close so that the competition between the strains is weak. Yet, one may expect that a substantial competitive advantage of one of the strains will eventually outbalance trade-off mechanisms driving less competitive strains to extinction. The literature, however, lacks a rigorous validation of this statement.
In this work, we challenge the validity of the exclusion principle at an ultimate limit in which one strain has a vast competitive advantage over the other strains. We show that when one strain is significantly more transmissible than the others, and under broad conditions, an epidemic system with two strains has a stable endemic equilibrium in which both strains coexist with comparable prevalence. Thus, the competitive exclusion principle does not unconditionally hold beyond the established case of complete immunity.

1 Introduction

One of the important principles of theoretical biology is the competitive exclusion principle which states that no two species can indefinitely occupy the same ecological niche [17, 14, 10, 28]. Many epidemiological systems, e.g., Dengue, Influenza or Malaria, consist of multiple strains of a pathogen which all compete for susceptibles, see [19, 15] and references within. The competitive exclusion principle in epidemiology implies that when competing strains of a pathogen provide complete protection for each other, the strain with the largest reproduction number outcompetes the other strains and drives them to extinction [2]. The result of [2] applies to directly transmitted infections with cross-immunity in a homogeneous population. Subsequent works showed that the exclusion principle extends to a broad family of epidemic systems including those described by vector-borne disease models [9, 32, 3], age-structured models [21], network structures [29], and also to models incorporating distributed delays and environmental transmission [3, 7].

In various cases, e.g., Dengue or Malaria, multiple strains of pathogens of infectious diseases stably coexist for long periods despite competing with each other, see [19, 15] and references within. Extensive research efforts have been dedicated to understanding how such competing species coexist rather than following the competitive exclusion principle and driving less competitive strains to extinction. These studies showed coexistence is possible when strain interactions and model dynamics extends beyond complete cross-immunity. Particularly, these studies revealed numerous trade-off mechanisms that may enable the coexistence of pathogen variants [20, 26].

Partial cross-immunity, for example, is a trade-off mechanism generating an immune response after infection that provides partial protection against subsequent infections with related strains, but not complete protection. As a result, competition is reduced, facilitating coexistence between strains [6]. Similarly, super-infection and co-infection reduce competition and facilitate coexistence by enabling an infected individual to be infected by other strains before recovery [22, 4]. Other examples of trade-off mechanisms involve adaptation by mutation, quarantine or age structure [25, 24, 26, 16].

The above studies of trade-off mechanisms, however, consider cases in which the basic reproduction numbers of the two strains are comparable, and, therefore, the competition between the two strains is relatively weak. For example, the study of a multi-strain system in [30] shows that the coexistence of strains is possible when their respective basic reproduction numbers are close and cross-immunity is weak. Similarly, a series of works shows how strain structure and partial cross-immunity between strains with close respective basic reproduction numbers can lead to long-period oscillatory dynamics without external forcing, see [15, Sec. 4.1.4.2.] and references within. The arising picture is that weaker competition implies weaker competitive exclusion, allowing various trade-off mechanisms to dominate dynamics.

Intuitively, the competitive exclusion principle should be most valid and robust in the limit of fierce competition between the strains. Thus, if one of the strains has a significantly larger competitive advantage over the other, one may expect it will outbalance trade-off mechanisms and become the dominant and eventually the sole strain. However, to the best of our knowledge, the literature lacks a rigorous validation of this statement. This is the focus of this study.

In this work, we test the validity of the exclusion principle in an ultimate regime in which one strain has a vast competitive advantage over the other strains, in the sense that one of the strains is significantly more transmissible than the other strains. We address the following questions:
1) Does the competitive exclusion principle always hold when one strain has a vast competitive advantage over the other strains? Namely, beyond the established case of complete cross-immunity, would a sufficiently large competitive advantage of one of the strain outbalance trade-off mechanisms?

Our findings are counter-intuitive and show that in a very general setting, the competitive exclusion principle does not hold, unconditionally, beyond the case of complete cross-immunity. Rather, we show that such systems may have a stable endemic equilibrium in which both species coexist. This raises further questions,

2) What trade-off mechanisms contribute to coexistence, even when one strain has a vast competitive advantage over the others?

Our results reveal the key role of partial cross-immunity in enabling coexistence.

2 Methods

2.1 Transmission model

Refer to caption — Figure 1: Schematic diagram of disease dynamics when the host is exposed to two co-circulating strains. Susceptibles (S) may be infected with strain $i=1$ or $i=2$ ( $I_{i}$ , primary infection). Those recovered from strain $i$ ( $R_{i}$ , as a result of primary infection) are temporarily immune to reinfections by strain $i$ but are possibly susceptible to infections by strain $j\neq i$ ( $J_{j}$ , secondary infection). Here, $\sigma_{ij}$ is the relative susceptibility to strain $j$ for an individual previously infected with and recovered from strain $i$ ( $i\neq j$ ), so that $\sigma_{ij}=0$ corresponds to total cross-immunity, $0<\sigma_{ij}<1$ corresponds to partial cross-immunity and $\sigma_{ij}>1$ corresponds to enhanced susceptibility.

We consider a two-strain epidemic model in which the strains can interact indirectly via the immune response generated following infections, and in which this immune response possibly wanes with time, see diagram in Figure 1. The mathematical model is based on the model presented in [5, 19], and extended to account for waning immunity and possible enhanced or reduced infectivity following an infection. Each individual resides in one of eight compartments: susceptibles (S), those infected with strain $i$ ( $I_{i}$ , primary infection), those infected with strain $i$ ( $J_{i}$ , secondary infection), after having been previously infected (and recovered) by strain $j\neq i$ , those recovered from strain $i$ ( $R_{i}$ , as a result of primary infection), and those recovered from both strains ( $R_{3}$ ). The population is assumed to mix randomly. The model is given by


$\displaystyle\frac{\text{d}S}{\text{d}t}$	$\displaystyle=\mu-\sum_{i=1}^{2}\beta_{i}(I_{i}+\eta_{i}J_{i})S-\mu S+\sum_{k=1}^{3}\delta_{k}R_{k},$	(1a)
$\displaystyle\frac{\text{d}I_{1}}{\text{d}t}$	$\displaystyle=\beta_{1}S(I_{1}+\eta_{1}J_{1})-(\mu+\gamma_{1})I_{1},$	(1b)
$\displaystyle\frac{\text{d}I_{2}}{\text{d}t}$	$\displaystyle=\beta_{2}S(I_{2}+\eta_{2}J_{2})-(\mu+\gamma_{2})I_{2},$	(1c)
$\displaystyle\frac{\text{d}J_{1}}{\text{d}t}$	$\displaystyle=\beta_{1}\sigma_{21}R_{2}(I_{1}+\eta_{1}J_{1})-(\mu+\gamma_{1})J_{1},$	(1d)
$\displaystyle\frac{\text{d}J_{2}}{\text{d}t}$	$\displaystyle=\beta_{2}\sigma_{12}R_{1}(I_{2}+\eta_{2}J_{2})-(\mu+\gamma_{2})J_{2},$	(1e)
$\displaystyle\frac{\text{d}{R_{1}}}{\text{d}t}$	$\displaystyle=\gamma_{1}I_{1}-\beta_{2}\sigma_{12}(I_{2}+\eta_{2}J_{2})R_{1}-(\mu+\delta_{1})R_{1},$	(1f)
$\displaystyle\frac{\text{d}{R_{2}}}{\text{d}t}$	$\displaystyle=\gamma_{2}I_{2}-\beta_{1}\sigma_{21}(I_{1}+\eta_{1}J_{1})R_{2}-(\mu+\delta_{2})R_{2},$	(1g)
$\displaystyle\frac{\text{d}R_{3}}{\text{d}t}$	$\displaystyle=\gamma_{1}J_{1}+\gamma_{2}J_{2}-(\mu+\delta_{3})R_{3},$	(1h)

were $\mu$ is the rate at which individuals are born, as well as the mortality rate, $\beta_{i}>0$ denotes the transmission coefficient for strain $i$ , $\gamma_{i}$ denotes the recovery rate from strain $i$ , $\delta_{i}$ is the rate at which immunity to re-infection by strain $i$ wanes and $\delta_{3}$ is the rate at which immunity to re-infection following infection by both strains wanes. The relative infectivity of those infected with strain $i$ after previous infection with strain $j\neq i$ is given by $\eta_{i}$ . Finally, $\sigma_{ij}$ is the relative susceptibility to strain $j$ for an individual previously infected with and recovered from strain $i$ ( $i\neq j$ ), so that $\sigma_{ij}=0$ corresponds to total cross-immunity, $0<\sigma_{ij}<1$ corresponds to partial cross-immunity and $\sigma_{ij}>1$ corresponds to enhanced susceptibility.

2.1.1 Assumptions on model parameters

The feasible range of the problem parameters is

	$\beta_{i}>0,\quad\gamma_{i}>0,\quad\mu\geq 0,\quad\delta_{k}\geq 0,\quad\eta_{i}>0,\quad\sigma_{ij}\geq 0,\qquad i=1,2,\quad j\neq i\quad k=1,2,3.$	(2a)
In what follows, we assume that the susceptible group is replenished by demographic turnover ( $\mu>0$ ) and/or by waning of the immune response generated following infections ( $\delta_{k}>0$ ),
	$\max\{\mu,\delta_{1},\delta_{2},\delta_{3}\}>0.$	(2b)

This condition enables the system to converge to an endemic equilibrium, rather than gradually exhausting the susceptible pool and converging to a disease-free state.

2.1.2 Basic reproduction number

A standard computation of the next generation matrix [19, 8, 27] yields that the basic reproduction number equals

\mathcal{R}_{0}=\max\{\mathcal{R}_{1},\mathcal{R}_{2}\}

where $\mathcal{R}_{i}$ is the basic reproduction number of strain $i$

\mathcal{R}_{i}=\frac{\beta_{i}}{\mu+\gamma_{i}}.

(3)

See Appendix A for details.

2.2 Analysis in the case of a highly transmissible strain

We focus on the case in which strain one has a vast competitive advantage over strain two, and ask whether, in such a case, the system can give rise to long-term dynamics in which both strains coexist. In particular, we consider the case in which strain one is significantly more transmissible than strain two, $\beta_{1}\gg 1$ , so that in terms of its basic reproduction number,

\mathcal{R}_{1}\gg 1.

(4)

The following Theorem characterizes single-strain endemic equilibrium points of (1) in the asymptotic regime (4).

Theorem 2.1.

Let $\mu$ , $\{\gamma_{i}\}_{i=1}^{2}$ , $\{\delta_{k}\}_{k=1}^{3}$ , $\{\eta_{i}\}_{i=1}^{2}$ , $\sigma_{12}$ and $\sigma_{21}$ be parameters that satisfy (2). Then, for sufficiently large $\beta_{1}>\gamma_{1}$ ,

The system (1) has a unique single-strain endemic equilibrium $\phi^{EE,1}$ with $I_{1}>0$ and $I_{2}=0$ that satisfies

\begin{split}&S^{EE,1}=\frac{1}{\mathcal{R}_{1}},\quad I_{1}^{EE,1}=\frac{\delta_{1}+\mu}{1+\delta_{1}+\mu}\frac{\mathcal{R}_{1}-1}{\mathcal{R}_{1}},\quad R_{1}^{EE,1}=\frac{1}{\delta_{1}+\mu}I_{1}^{EE,1},\quad I_{2}^{EE,1}=J_{1}^{EE,1}=J_{2}^{EE,1}=R_{2}^{EE,1}=0,\end{split}

where $\mathcal{R}_{i}$ is given by (3).

The solution $\phi^{EE,1}$ is locally unstable when

\hat{\mathcal{R}}^{1}_{2}:=\frac{\sigma_{12}\gamma_{1}\eta_{2}}{\gamma_{1}+\delta_{1}+\mu}\mathcal{R}_{2}+\left[1-\frac{\sigma_{12}\gamma_{1}\eta_{2}}{\gamma_{1}+\delta_{1}+\mu}\right]\frac{\mathcal{R}_{2}}{\mathcal{R}_{1}}>1,

(5)

and locally stable when $\hat{\mathcal{R}}^{1}_{2}<1$ .

If $\mathcal{R}_{2}>1$ , the system (1) has an additional unique single-strain endemic equilibrium $\phi^{EE,2}$ with $I_{2}>0$ and $I_{1}=0$ which satisfies

\begin{split}&S^{EE,2}=\frac{1}{\mathcal{R}_{2}},\quad I_{2}^{EE,2}=\frac{\delta_{2}+\mu}{1+\delta_{2}+\mu}\frac{\mathcal{R}_{2}-1}{\mathcal{R}_{2}},\quad R_{2}^{EE,2}=\frac{1}{\delta_{2}+\mu}I_{1}^{EE,2},\quad I_{1}^{EE,2}=J_{1}^{EE,2}=J_{2}^{EE,2}=R_{1}^{EE,2}=0.\end{split}

This solution is unstable.

Proof.

See Appendix B. ∎

Theorem 2.1 implies that when $\mathcal{R}_{1}\gg 1$ and the invasion number, $\hat{\mathcal{R}}^{1}_{2}$ , of strain two is bigger than one, then both single-strain endemic equilibrium points, if exist, are unstable. The disease-free equilibrium is also unstable in the regime $\mathcal{R}_{1}\gg 1$ , see Section 2.1.2. Thus, the system may potentially give rise to coexistence in the regime $\mathcal{R}_{1}\gg 1$ and $\hat{\mathcal{R}}^{1}_{2}>1$ . However, instability of the disease-free equilibrium and the single-strain endemic equilibria does not assure the persistence of the two strains [18]. Moreover, if system (1) does give rise to coexisting strains in the long run, Theorem 2.1 does not provide any information on the nature of coexistence, e.g., does the system oscillate between infection waves of different strains, converge to an endemic equilibrium or develop chaos.

The following theorem complements the results of Theorem 2.1 by characterizing cases in which the system of study has a stable endemic equilibrium in which both strains coexist, despite the vast competitive advantage of strain one:

Theorem 2.2.

Let $\beta_{2}$ , $\mu$ , $\{\gamma_{i}\}_{i=1}^{2}$ , $\{\delta_{k}\}_{k=1}^{3}$ , $\{\eta_{i}\}_{i=1}^{2}$ , $\sigma_{12}$ and $\sigma_{21}$ be parameters that satisfy (2). If, $\sigma_{12}>0$ and

\mathcal{R}_{2}>\frac{\gamma_{1}+\delta_{1}+\mu}{\gamma_{1}\eta_{2}\sigma_{12}},

(6)

where $\mathcal{R}_{i}$ is given by (3), then for sufficiently large $\beta_{1}>\gamma_{1}$
1) There exists a unique steady-state solution

\phi^{*}=(S^{*},I_{1}^{*},I_{2}^{*},J_{1}^{*},J_{2}^{*},R_{1}^{*},R_{2}^{*},R_{3}^{*})

of the system (1) with $I_{1}^{*}>0$ and $I_{2}^{*}>0$ that satisfies

	$S^{}=\frac{1}{\mathcal{R}_{1}}-\frac{s_{2}}{\mathcal{R}_{1}^{2}}+\mathcal{O}\left(\frac{1}{\mathcal{R}_{1}^{3}}\right),\quad\quad I_{1}^{}=a_{1}-\frac{b_{1}}{\mathcal{R}_{1}}+\mathcal{O}\left(\frac{1}{\mathcal{R}_{1}^{2}}\right),\quad I_{2}^{*}=\frac{b_{2}}{\mathcal{R}_{1}}+\mathcal{O}\left(\frac{1}{\mathcal{R}_{1}^{2}}\right),$	(7a)
	$J_{2}^{}=\frac{\gamma_{1}}{\gamma_{2}+\mu}I_{1}^{}-\frac{\mu+\delta_{1}}{\gamma_{2}+\mu}R_{1}^{}=\frac{\gamma_{1}a_{1}}{\gamma_{2}+\mu}-\frac{\delta_{1}+\mu}{\eta_{2}(\gamma_{2}+\mu)\sigma_{12}\mathcal{R}_{2}}+\mathcal{O}\left(\frac{1}{\mathcal{R}_{1}}\right),\qquad R^{}_{1}=\frac{1-\mathcal{R}_{2}S^{*}}{\eta_{2}\mathcal{R}_{2}\sigma_{12}}=\frac{1}{\eta_{2}\sigma_{12}\mathcal{R}_{2}}+\mathcal{O}\left(\frac{1}{\mathcal{R}_{1}}\right)$	(7b)
	$b_{2}=\frac{\gamma_{1}\eta_{2}(\delta_{3}+\mu)}{\gamma_{1}\gamma_{2}+(\gamma_{1}+\gamma_{2}+\mu)(\delta_{3}+\mu)}\left[\mathcal{R}_{2}-\frac{\gamma_{1}+\delta_{1}+\mu}{\gamma_{1}\eta_{2}\sigma_{12}}\right],\quad a_{1}=\frac{(\gamma_{2}+\mu)\sigma_{12}b_{2}+\delta_{1}+\mu}{\gamma_{1}\eta_{2}\sigma_{12}\mathcal{R}_{2}},$	(7c)
where for $\sigma_{21}>0$
	$J_{1}^{}=\frac{\gamma_{2}}{\gamma_{1}+\mu}I_{2}^{}-\frac{\mu+\delta_{2}}{\gamma_{1}+\mu}R_{2}^{}=\frac{\gamma_{2}b_{2}}{\gamma_{1}+\mu}\frac{1}{\mathcal{R}_{1}}+\mathcal{O}\left(\frac{1}{\mathcal{R}_{1}^{2}}\right),\qquad R_{2}^{}=\frac{1-\mathcal{R}_{1}S^{*}}{\eta_{1}\mathcal{R}_{1}\sigma_{21}}=\frac{\gamma_{2}b_{2}}{a_{1}\sigma_{21}(\gamma_{1}+\mu)}\frac{1}{\mathcal{R}_{1}^{2}}+\mathcal{O}\left(\frac{1}{\mathcal{R}_{1}^{3}}\right),$	(7d)
	$s_{2}=\frac{\gamma_{2}\eta_{1}b_{2}}{(\mu+\gamma_{1})a_{1}},\quad b_{1}=\frac{\gamma_{2}+\mu}{\gamma_{1}+\mu}b_{2}+\frac{\sigma_{12}(\gamma_{2}+\mu)(\delta_{3}+\mu)+(\gamma_{2}+\delta_{3}+\mu)\delta_{1}-\gamma_{2}\delta_{3}}{(\gamma_{1}+\gamma_{2})\delta_{3}+\gamma_{1}\gamma_{2}+(\gamma_{1}+\gamma_{2}+\delta_{3})\mu+\mu^{3}}\frac{1}{\sigma_{12}},$	(7e)
and for $\sigma_{21}=0$ ,
	$J_{1}^{}=0,\qquad R_{2}^{}=\frac{\gamma_{2}}{\delta_{2}+\mu}I_{2}^{*}=\frac{\gamma_{2}b_{2}}{\delta_{2}+\mu}\frac{1}{\mathcal{R}_{1}}+\mathcal{O}\left(\frac{1}{\mathcal{R}_{1}^{2}}\right),$	(7f)
	$s_{2}=0,\quad b_{1}=\frac{(\gamma_{2}+\mu)(\delta_{3}+\mu)(\gamma_{2}+\delta_{2}+\mu)\eta_{2}\sigma_{12}b_{2}+(\mu+\delta_{2})\left[(\gamma_{2}+\mu)(\delta_{3}+\mu)\eta_{2}\sigma_{12}+\gamma_{2}(\delta_{1}-\delta_{3})+\delta_{1}(\delta_{3}+\mu)\right]}{[(\gamma_{1}+\gamma_{2})\delta_{3}+\gamma_{1}\gamma_{2}+(\gamma_{1}+\gamma_{2}+\delta_{3})\mu+\mu^{3}](\delta_{2}+\mu)\eta_{2}\sigma_{12}}.$	(7g)

2) The solution $\phi^{*}$ is locally stable.

Proof.

See appendix C ∎

Theorem 2.2 considers cases when strain one is highly transmissible, and strain two is transmissible enough to satisfy (6). In these cases, the theorem proves the existence, uniqueness and stability of an endemic equilibrium in which both strains coexist. Condition (6) identifies with condition (5) in the limit $\beta_{1}\to\infty$ and can be read as

\lim_{\beta_{1}\to\infty}\hat{\mathcal{R}}^{1}_{2}>1,

i.e., the invasion number of strain two is bigger than one when strain one is highly transmissible. We note that condition (6) is stronger than condition (5). Namely, the coexistence endemic equilibrium $\phi^{*}$ does not exist when the single-strain endemic equilibrium $I_{EE}^{1}$ is stable. Furthermore, the theorem provides an approximation to the variable values at the endemic equilibrium, see (7). Condition (6), for example, originates from approximation (7). Indeed, the requirement $I_{2}^{*}>0$ implies $b_{2}>0$ , which in turn leads to condition (6), see (7a) and (7c). Approximation (7) also implies that despite the vast competitive advantage of strain one, at the equilibrium of coexistence, the prevalence of both strains is comparable

\frac{I_{2}^{*}+J_{2}^{*}}{I_{1}^{*}+J_{1}^{*}}=\frac{\gamma_{1}\sigma_{12}b_{2}}{(\gamma_{2}+\mu)\sigma_{12}b_{2}+\delta_{1}+\mu}+\mathcal{O}\left(\frac{1}{\mathcal{R}_{1}}\right),

and the prevalence of strain two may exceed that of strain one.

3 Results

Theorem 2.2 implies that two strains can coexist despite fierce competition between the strains. In the system of study, a necessary condition for coexistence is that infection by the dominant strain does not provide complete protection from infections by the other strain. Indeed, in the limit of complete cross-immunity, $\sigma_{12}=0$ , condition (6) cannot be satisfied. Therefore, Theorem 2.2 implies that the competitive exclusion principle does not unconditionally hold beyond the established case of complete cross-immunity [2], even if one strain has a vast competitive advantage over the others.

In what follows we present two numerical examples. In both cases, we fix all parameters except the transmission parameter of strain one and study the system’s behavior as strain one becomes more and more transmissible. We choose the timescales so that they roughly correspond to Influenza or COVID-19, where one unit of time is a week: A typical infectious period of one week ( $\gamma_{i}\approx 1$ ), characteristic convalescent immunity period of a year ( $\delta_{i}\approx 1/50$ ) [31, 13], and lifespan of 80 years ( $\mu\approx 1/4000$ ). In both examples, the parameters are chosen so that condition (6) of Theorem 2.2 is satisfied. Hence, we expect the systems to converge to a multi-strain endemic equilibrium when strain one is sufficiently transmissible.

We first consider an example in which the strain properties, excluding transmissibility, are comparable with only slight differences in their values. We chose the parameters so that they correspond to partial cross-immunity $0<\sigma_{ij}<1$ and neutral or reduced relative infectivity $0<\eta_{i}\leq 1$ in secondary infections. Since these values do not impose any advantage for secondary infections, one may expect that they would not facilitate coexistence in the case of fierce competition between the strains. We chose the initial condition $I_{1}(0)=10^{-2}$ and $I_{2}(0)=10^{-8}$ so that the initial prevalence of strain two is significantly smaller than the initial prevalence of strain one. This choice sets the initial values far from the equilibrium $\phi^{*}$ of coexistence and tests for a possible invasion of strain two. Figure 2A presents the prevalence of the strains, $(I_{1}+J_{1})(t)$ and $(I_{2}+J_{2})(t)$ , in the case both strains have the same basic reproduction number $\mathcal{R}_{1}=\mathcal{R}_{2}=1.5$ . We observe that strain two invades and the system reaches an endemic equilibrium in which both strains coexist with a comparable prevalence of the two strains, $I_{1}^{*}+J_{1}^{*}\approx 4\cdot 10^{-3}$ and $I_{2}^{*}+J_{2}^{*}\approx 7.2\cdot 10^{-3}$ . Next, we consider the same system when strain one is more transmissible. As can be expected from Theorem 2.2, we observe that strain two invades and co-exists with strain one when $\mathcal{R}_{1}=5$ and $\mathcal{R}_{1}=10$ , see Figures 2B and 2C, respectively. Furthermore, the prevalence of the strains is comparable with $I_{1}^{*}+J_{1}^{*}\approx 1.5\cdot 10^{-2}$ and $I_{2}^{*}+J_{2}^{*}\approx 0.5\cdot 10^{-2}$ .

To demonstrate that partial cross-immunity is necessary for coexistence, as follows from Theorem 2.2, we solve (1) in the case of complete cross-immunity, $\sigma_{12}=0$ , and otherwise with the same parameters as in Figure 2. As expected, we observe that strain two is driven to extinction for both moderate and large values of $\mathcal{R}_{1}$ , see Figure 3.

Theorem 2.2 proves the local stability of $\phi^{*}$ under appropriate conditions, and leaves open the question of whether $\phi^{*}$ is globally stable. In Figure 2, the initial condition is set far from $\phi^{*}$ . The less competitive strain successfully invades and the system converges to the coexistence equilibrium $\phi^{*}$ . These results, as well as results from additional simulations (data not shown), suggest that the equilibrium is globally stable.

In the next example, we consider a case in which the basic reproduction number of strain two is smaller than one, yet the relative infectivity in secondary infections with strain two, $\eta_{2}=3$ , is set relatively high so that condition (6) is satisfied. In this case, strain two cannot persist alone, yet Theorem 2.2 does imply that it will co-exist with strain one when strain one is sufficiently transmissible. As expected, for moderate values of $\mathcal{R}_{1}$ , strain two is driven to extinction, see Figure 4A. Yet, for $\mathcal{R}_{1}=5$ and $\mathcal{R}_{1}=10$ , the two strains coexist, see Figures 4B and 4C, respectively. Intuitively, strain two benefits from the co-circulation of strain one since the infectivity of those infected with strain two after being infected by strain one is higher than the infectivity of those with primary infections with strain two. These simulations also provide an example of a case in which the prevalence of the less dominant strain is higher in the long term than the prevalence of the dominant strain, see Figures 4B and 4C.

A systematic numerical verification of approximation (7) is presented in Appendix D.

4 Conclusions

Our results imply that the competitive exclusion principle does not extend beyond the established case of complete cross-immunity, even when one species has a vast competitive advantage over the other. Rather, we show that the system can converge to an endemic equilibrium in which both strains coexist, despite fierce competition between the strains. Particularly, we provide an example of coexistence when the basic reproduction number of the less dominant strain is less than one, and therefore this strain cannot persist without an interaction with the dominant strain. Moreover, we show that, in equilibrium, the prevalence of the less dominant strain is comparable in magnitude to that of the dominant strain, and even exceeds it in certain cases.

To understand why a vast competitive advantage of strain one does not drive strain two to extinction, it is useful to consider an epidemic outbreak in which both strains are introduced to a susceptible population. Due to the high basic reproduction number of strain one, in the epidemic’s early stages, infections by strain one are dominant and widespread. While those recovered from strain one remain immune to reinfections by strain one for some time, in the absence of complete cross-immunity, they remain susceptible, to some degree, to infections by strain two. Strain two can potentially thrive on this susceptible pool without suffering from competition with strain one. More concretely, strain two will be able to spread if its effective reproduction number $\mathcal{R}_{2}^{\rm eff}$ in a population that has all been infected with strain one is larger than one

\mathcal{R}_{2}^{\rm eff}=\sigma_{12}\mathcal{R}_{2}\eta_{2}>1.

Here, $\mathcal{R}_{2}^{\rm eff}$ equals the basic reproduction number $\mathcal{R}_{2}$ of strain two (for primary infections) while accounting for the relative changes in susceptibility, $\sigma_{12}$ , and infectivity, $\eta_{2}$ of those recovered from strain one. This condition is closely related to condition (6) which also takes the effect of waning immunity, $\delta_{1}$ , mortality, $\mu$ , and the infectious period, $\gamma_{1}$ , on the size of the population immune to infections by strain one, and can be written as,

\frac{\gamma_{1}}{\gamma_{1}+\delta_{1}+\mu}\mathcal{R}_{2}^{\rm eff}>1.

Following the explanation above, one may also view strain two as an invading strain after strain one widely spread in the population and reached equilibrium. Under this view, we note that condition (6) can also be read as the asymptotic invasion number $\hat{\mathcal{R}^{1}_{2}}$ of strain two at the limit of large transmissibility of strain one, see (5).

The analysis in this work combines two approaches. First, we follow a common approach and focus on the behavior of single-strain equilibriums while taking advantage of the fact that it is easier to compute the single-strain equilibrium values and determine their local stability, see Theorem 2.1. Yet, the instability of the single-strain and disease-free equilibrium points does not imply the persistence of multiple strains. Proving persistence may be demanding in complex epidemic systems with interacting strains [18], and attaining such a result does not reveal the nature of coexistence in the system, e.g., whether the system converges to an endemic equilibrium of coexistence or oscillates. Therefore, we complement our analysis by studying the possible equilibrium(s) of coexistence. Since this analysis involves problems that do not have an analytic solution or have a solution that is cumbersome to work with, we use asymptotic methods that enable us to tackle these problems at the accuracy required to attain the solution, e.g., determine the stability of the coexistence state. Using this approach we show the existence and stability of an endemic equilibrium in which both strains coexist. Therefore, the analysis not only shows coexistence but also characterizes its nature.

Our results show the local stability of an endemic equilibrium of coexistence. Numerical simulations show that the system converges to the above endemic equilibrium even when the initial conditions are far from equilibrium, suggesting that the equilibrium is globally stable. Proving global stability is left open as an interesting mathematical question for future research.

One implication of our results is that in the absence of complete cross-immunity, less transmissible, yet more dangerous, strains can invade and coexist with highly transmissible strains. This suggests that the evolution of pathogens and the emergence of new strains may not maximize $\mathcal{R}_{0}$ with time, but rather other quantities. Such a phenomenon was demonstrated in an epidemic system with superinfection [23]. The study of evolution in a multi-strain system with partial cross-immunity will be presented elsewhere.

Our approach challenged the validity of the exclusion principle in an ultimate limit where, intuitively, it should be most valid. From a mathematical point of view, studying such a limit enables using approximation methods that allow analysis of the system despite its complexity. Indeed, many details become negligible at the limit, making the analysis tractable. Concurrently, other quantities arise naturally at the limit. Thus, from a biological point of view, such an approach provides valuable insights as it highlights quantities of interest while naturally sorting out details. The unexpected results revealed using this approach refine our understanding of the exclusion principle and multi-strain dynamics.

This study was motivated by our experience during the SARS-CoV-2 pandemic. Indeed, the emergence of the Omicron variant created an unprecedented scenario of an epidemic driven by a highly transmissible variant. Our data-driven research during that time [12] showed that a variant with a basic reproduction number as high as ten can defy conventional theory in certain circumstances. One example is counter-intuitive optimal vaccine allocations at elevated reproduction numbers [11]. The results presented in this work provide an additional example to an unexpected result that arises when considering highly transmissible strains.

The counter-intuitive results presented in this study touch upon a fundamental principle in epidemiology and shed light on the mechanisms that enable competing strains to coexist. While we have focused on the competitive exclusion principle in epidemiology, it is plausible that similar phenomena will also occur in other ecological systems. For example, microbial populations that evolve near plant roots, and in particular accumulation of host-specific pathogens, may generate negative feedback on the plant and facilitate diversity or coexistence of plant species [1], plausibly, even if one plant has a vast competitive advantage over the other. Extensions of our study beyond epidemiological systems will be presented elsewhere.

Acknowledgements

This research was supported by the Israel Science Foundation (grant no. 3730/20 to N.G.) within the KillCorona-Curbing Coronavirus Research Program, and by the Israeli Science Foundation (ISF) grant 1596/23.

We thank Guy Katriel for the most helpful comments and discussions on this work.

Appendix A Computation of $\mathcal{R}_{0}$

Following a standard computation of the next generation matrix [19, 8, 27] we split the right-hand side in the equations for the infected compartments $(I_{1},I_{2},J_{1},J_{2})$ of (1) in the following way:

\begin{split}\frac{\text{d}I_{1}}{\text{d}t}&=\underbrace{\beta_{1}S(I_{1}+\eta_{1}J_{1})}_{\mathcal{F}_{1}}-\underbrace{(\mu+\gamma_{1})I_{1}}_{\mathcal{V}_{1}},\\ \frac{\text{d}I_{2}}{\text{d}t}&=\underbrace{\beta_{2}S(I_{2}+\eta_{2}J_{2})}_{\mathcal{F}_{2}}-\underbrace{(\mu+\gamma_{2})I_{2}}_{\mathcal{V}_{2}},\\ \frac{\text{d}J_{1}}{\text{d}t}&=\underbrace{\beta_{1}\sigma_{21}R_{2}(I_{1}+\eta_{1}J_{1})}_{\mathcal{F}_{3}}-\underbrace{(\mu+\gamma_{1})J_{1}}_{\mathcal{V}_{3}},\\[4.30554pt] \frac{\text{d}J_{2}}{\text{d}t}&=\underbrace{\beta_{2}\sigma_{12}R_{1}(I_{2}+\eta_{2}J_{2})}_{\mathcal{F}_{4}}-\underbrace{(\mu+\gamma_{2})J_{2}}_{\mathcal{V}_{4}},\end{split}

(8)

where $\mathcal{F}_{i}$ is the appearance rate of new infections in the relevant compartments, and $\mathcal{V}_{i}$ incorporates the remaining transitional terms. Linearization of the system of infected compartments, (8), about the disease-free equilibrium in which all the population is susceptible, $S=1$ , yields

x^{\prime}(t)=(F-V)x,\quad x=(I_{1},I_{2},J_{1},J_{2})^{T},

where the matrices $F$ and $V$ originated from the terms $\mathcal{F}_{i}$ and $\mathcal{V}_{i}$ , respectively,

F=\left[\begin{array}[]{cccc}\beta_{1}&&\beta_{1}\eta_{1}&\\ &\beta_{2}&&\beta_{2}\eta_{2}\\ 0&0&0&0\\ 0&0&0&0\end{array}\right],\quad V=\left[\begin{array}[]{cccc}\gamma_{1}+\mu&&&\\ &\gamma_{2}+\mu&&\\ &&\gamma_{1}+\mu&\\ &&&\gamma_{2}+\mu\end{array}\right].

The basic reproduction number is the maximal eigenvalue of the next-generation matrix

\mathcal{R}_{0}=\rho(FV^{-1})=\max\{\mathcal{R}_{1},\mathcal{R}_{2}\},\quad\mathcal{R}_{i}=\frac{\beta_{i}}{\gamma_{i}+\mu}.

Appendix B Proof of Theorem 2.1

We first prove part (i) of Theorem 2.1 by seeking a steady-state solution of (1) with $I_{1}>0$ and $I_{2}=0$ . Substituting $I_{2}=0$ in $J_{1}^{\prime}(t)+R_{2}^{\prime}(t)=0$ , see (1d) and (1g), implies that $R_{2}=0$ . Thus, $J_{1}^{\prime}(t)=0$ implies $J_{1}=0$ , see (1d). Finally, $I_{2}^{\prime}=0$ implies $S=0$ or $J_{2}=0$ , see (1c). In the former case, however, (1b) implies that $I_{1}=0$ at steady-state. Thus, in the case of interest, $I_{1}>0$ , it follows that $S\neq 0$ and $J_{2}=0$ . Overall, the above computation shows that the model (1) has a unique single-strain endemic equilibrium point $\phi^{EE,1}$ with $I_{1}>0$ and $I_{2}=0$ that satisfies

I_{2}^{EE,1}=J_{1}^{EE,1}=J_{2}^{EE,1}=R_{2}^{EE,1}=0.

Substituting the above in (1) and solving the resulting algebraic system yields that the other variables of $\phi^{EE,1}$ uniquely satisfy

S^{EE,1}=\frac{1}{\mathcal{R}_{1}},\quad I_{1}^{EE,1}=\frac{\delta_{1}+\mu}{1+\delta_{1}+\mu}\frac{\mathcal{R}_{1}-1}{\mathcal{R}_{1}},\quad R_{1}^{EE,1}=\frac{1}{\delta_{1}+\mu}I_{1}^{EE,1}.

Note that the condition $I_{1}^{EE,1}>0$ implies $\mathcal{R}_{1}>1$ . This condition is satisfied for sufficiently large $\beta_{1}>\gamma_{1}$ . Similarly, when $\mathcal{R}_{2}>1$ , model (1) has a unique single-strain endemic steady-state $\phi^{EE,2}$ with $I_{2}>0$ and $I_{1}=0$ that satisfies

\begin{split}&S^{EE,2}=\frac{1}{\mathcal{R}_{2}},\quad I_{2}^{EE,2}=\frac{\delta_{2}+\mu}{\gamma_{2}+\delta_{2}+\mu}\frac{\mathcal{R}_{2}-1}{\mathcal{R}_{2}},\quad R^{EE,2}_{2}=\frac{\gamma_{2}}{\delta_{2}+\mu}I^{EE,2}_{2},\\ &I^{EE,2}_{1}=J^{EE,2}_{1}=J_{2}^{EE,2}=R_{1}^{EE,2}=0.\end{split}

To study the stability of $\phi_{EE}^{1}$ , we compute the invasion number $\hat{\mathcal{R}}^{1}_{2}$ of strain two at the equilibrium $\phi_{EE}^{1}$ of strain one. To do so, we compute the maximal eigenvalue of the next-generation matrix of strain two: Similar to the computation in A, we linearize the relevant parts in the system of infected compartments, (8), about $\phi^{EE,1}$ to obtain the next-generation matrix and then compute its’ maximal eigenvalue. The resulting invasion number is given by (5), see Theorem 2.1(ii).

Similarly, the invasion number $\hat{\mathcal{R}}^{2}_{1}$ of strain one at the equilibrium $\phi_{EE}^{2}$ of strain two is given by

\hat{\mathcal{R}}^{2}_{1}:=\left[1+\frac{\sigma_{21}\gamma_{2}\eta_{1}}{\gamma_{2}+\delta_{2}+\mu}(\mathcal{R}_{2}-1)\right]\frac{\mathcal{R}_{1}}{\mathcal{R}_{2}}.

For sufficiently large $\mathcal{R}_{1}$ , $\phi_{EE}^{2}$ is unstable since $\hat{\mathcal{R}}^{2}_{1}>1$ . Thus proving part iii.

Appendix C Proof of Theorem 2.2

The following Lemma reduces the computation of the equilibrium point $\phi^{*}$ from a solution to a nonlinear system of seven equations, determined by the steady-state of (1), to a solution of a nonlinear system of three equations for three unknown $(S^{*},I_{1}^{*},I_{2}^{*})$ :

Lemma C.1.

Let $\sigma_{12}>0$ , and let $\phi^{*}=(S^{*},I_{1}^{*},I_{2}^{*},J_{1}^{*},J_{2}^{*},R_{1}^{*},R_{2}^{*})$ be a steady-state solution of the system (1) with $I_{1}^{*}>0$ and $I_{2}^{*}>0$ . Then, $\{S^{*},I_{1}^{*},I_{2}^{*}\}$ satisfy


	$\displaystyle\mu-\sum_{i=1}^{2}\beta_{i}[I_{i}^{}+\eta_{i}J_{i}^{}(S^{},I_{1}^{},I_{2}^{})]S^{}-\mu S^{}+\sum_{k=1}^{3}\delta_{k}R_{k}^{}(S^{},I_{1}^{},I_{2}^{*})=0,$	(9a)
	$\displaystyle\beta_{1}S^{}[I_{1}^{}+\eta_{1}J_{1}^{}(S^{},I_{2}^{})]-(\mu+\gamma_{1})I_{1}^{}=0,$	(9b)
	$\displaystyle\beta_{2}S^{}[I_{2}^{}+\eta_{2}J_{2}^{}(S^{},I_{1}^{})]-(\mu+\gamma_{2})I_{2}^{}=0,$	(9c)
where
	$J_{1}^{}(S^{},I_{2}^{})=\begin{cases}\frac{\gamma_{2}}{\gamma_{1}+\mu}I_{2}^{}-\frac{\mu+\delta_{2}}{\gamma_{1}+\mu}R_{2}^{}(S^{})&\sigma_{21}>0\\ 0&\sigma_{21}=0\end{cases},\quad J_{2}^{}(S^{},I_{1}^{})=\frac{\gamma_{1}}{\gamma_{2}+\mu}I_{1}^{}-\frac{\mu+\delta_{1}}{\gamma_{2}+\mu}R_{1}^{}(S^{}).$	(9d)
and
	$R^{}_{1}(S^{})=\frac{1-\mathcal{R}_{2}S^{}}{\eta_{2}\mathcal{R}_{2}\sigma_{12}},\qquad R_{2}^{}(S^{})=\begin{cases}\frac{1-\mathcal{R}_{1}S^{}}{\eta_{1}\mathcal{R}_{1}\sigma_{21}},&\sigma_{21}>0\\ \frac{\gamma_{2}}{\delta_{2}+\mu}I_{2}^{*},&\sigma_{21}=0,\end{cases}.$	(9e)

Proof.

When $I_{1}^{*}>0$ and $I_{2}^{*}>0$ , Equations (1b) and (1c) imply that the steady-state solution satisfies

\frac{\eta_{1}J_{1}^{*}}{I_{1}^{*}+\eta_{1}J_{1}^{*}}=1-\mathcal{R}_{1}S^{*},\quad\frac{\eta_{2}J^{*}_{2}}{I^{*}_{2}+\eta_{2}J^{*}_{2}}=1-\mathcal{R}_{2}S^{*}.

When $\sigma_{21}>0$ , substituting the above expressions in (1d) and (1e) defines $R^{*}_{i}$ as a linear function of $S^{*}$

\eta_{1}\mathcal{R}_{1}\sigma_{21}R_{2}^{*}(S^{*})=\frac{\eta_{1}J_{1}^{*}}{I_{1}^{*}+\eta_{1}J^{*}_{1}}=1-\mathcal{R}_{1}S^{*},\quad\eta_{2}\mathcal{R}_{2}\sigma_{12}R^{*}_{1}(S^{*})=\frac{\eta_{2}J_{2}^{*}}{I_{2}^{*}+\eta_{2}J_{2}^{*}}=1-\mathcal{R}_{2}S^{*},

yielding (9e) for $\sigma_{21}>0$ . Furthermore, the steady-state solution satisfies $J_{i}^{\prime}(t)=R_{i}^{\prime}(t)=0$ for $i=1,2$ , hence $J_{1}^{\prime}(t)+R_{2}^{\prime}(t)=0=J_{2}^{\prime}(t)+R_{1}^{\prime}(t)$ , implying (9d). Substituting (9e) and (9d) in (1a-1c), and setting $S^{\prime}(t)=I_{1}^{\prime}(t)=I_{2}^{\prime}(t)=0$ gives rise to (9a-9c).

When $\sigma_{21}=0$ , the above computation remains valid for $J_{2}^{*}$ and $R_{1}^{*}$ , while equations (1d) and (1g) directly imply $J_{1}^{*}=0$ and $(\delta_{2}+\mu)R_{2}^{*}=\gamma_{2}I_{2}^{*}$ . Therefore, giving rise to (9d) and (9e) for $\sigma_{21}=0$ . ∎

We now rely on Lemma C.1 to prove the first part of Theorem 2.2. For brevity, we consider separately the cases $\sigma_{21}>0$ and $\sigma_{21}=0$ :

Proposition C.1 (Existence, uniqueness and approximation of $\phi^{*}$ for $\sigma_{21}>0$ ).

Let $\sigma_{21}>0$ . Under the conditions of Theorem 2.2, if $\sigma_{12}>0$ , then there exists a unique steady-state solution

\phi^{*}=(S^{*},I_{1}^{*},I_{2}^{*},J_{1}^{*},J_{2}^{*},R_{1}^{*},R_{2}^{*})

of the system (1) with $I_{1}^{*}>0$ and $I_{2}^{*}>0$ . This solution satisfies (7).

Proof.

By Lemma C.1, $\phi^{*}$ satisfies (9). Denote $\varepsilon=1/\mathcal{R}_{1}$ and consider an expansion of the form

S^{*}=s_{0}+\varepsilon s_{1}+\varepsilon^{2}s_{2}+\mathcal{O}(\varepsilon^{3}),\quad I^{*}_{i}=a_{i}+\varepsilon b_{i}+\mathcal{O}(\varepsilon^{2}).

(10)

Substituting (10) into (9a) and equating the $\mathcal{O}(1/\varepsilon)$ terms in the arising expansion yields

[((\gamma_{1}+\mu)a_{1}+\eta_{1}\gamma_{2}a_{2})\sigma_{21}+(\mu+\delta_{2})s_{0}]s_{0}=0,

Since $a_{1},a_{2},s_{0}\geq 0$ , and due to (2), it follows that $s_{0}=0$ . Further substituting (10) with $s_{0}=0$ in (9c), yields

(\gamma_{2}+\mu)a_{2}=\mathcal{O}(\varepsilon),

hence $a_{2}=0$ . Similarly, substituting (10) in (9b) and in (9c) yields

(\gamma_{1}+\mu)(s_{1}-1)a_{1}=\mathcal{O}(\varepsilon),\quad\left[a_{1}\sigma_{12}\mathcal{R}_{2}\eta_{2}\gamma_{1}-(\delta_{1}+\mu)\right]s_{1}-(\gamma_{2}+\mu)\sigma_{12}b_{2}=\mathcal{O}(\varepsilon),

respectively. Since $s_{1},a_{1},b_{2}\geq 0$ , the only feasible solution of the above system is

s_{1}=1,\quad a_{1}=\frac{(\gamma_{2}+\mu)\sigma_{12}b_{2}+(\delta_{1}+\mu)}{\mathcal{R}_{2}\eta_{2}\gamma_{1}}.

At the next order, the expansion of (9b) yields

s_{2}=-\frac{\gamma_{2}\eta_{1}b_{2}}{(\gamma_{1}+\mu)a_{1}}.

Finally, substituting (10) into (9a) yields (7e).

Requiring $I_{2}^{*}>0$ implies $b_{2}>0$ , giving rise to condition (6), see (7c). This condition also ensures that, for sufficiently large $\beta_{1}>\gamma_{1}$ , $I_{1}^{*}>0$ since $a_{1}>0$ when $b_{2}>0$ , see (7c).

The above computation uniquely determined the values of the coefficients in (7). We now use the implicit function theorem to show that there exists a unique solution to (9a), which can be expressed as an expansion of the form (10) up to arbitrary order: Let ${\bf F}(\tilde{S}^{*}(\varepsilon),I_{1}^{*}(\varepsilon),I_{2}^{*}(\varepsilon);\varepsilon)$ be the algebraic system (9) where $S^{*}=\frac{\gamma_{1}+\mu}{\beta_{1}}(1+\tilde{S}^{*})$ . Then, as follows from the analysis above, ${\bf F}(\tilde{S}^{*},I_{1}^{*},I_{2}^{*};\varepsilon=0)$ has a unique solution $\tilde{S}^{*}(0)=0$ , $I_{1}^{*}(0)=a_{1}$ , $I_{2}^{*}(0)=0$ . The Jacobian matrix ${\bf J}$ of (9a) at $\varepsilon=0$ equals

\left.{\bf J}\right|_{(\tilde{S}^{*}(0),I_{1}^{*}(0),I_{2}^{*}(0))}=\left[\begin{array}[]{ccc}-a_{1}(\gamma_{1}+\mu)&-(\gamma_{1}+\mu)-\delta_{3}\frac{\gamma_{1}+\gamma_{2}+\mu}{\gamma_{2}+\mu}&-\eta_{1}\gamma_{2}-\delta_{3}\frac{\gamma_{1}+\gamma_{2}+\mu}{\gamma_{2}+\mu}\\ a_{1}(\gamma_{1}+\mu)&0&\eta_{1}\gamma_{2}\\ 0&0&-(\gamma_{2}+\mu)\end{array}\right].

The determinant of ${\bf J}$ equals

\mbox{Det}({\bf J})=-(\gamma_{1}+\mu)(\mu^{2}+(\delta_{3}+\gamma_{1}+\gamma_{2})\mu+(\delta_{3}+\gamma_{2})\gamma_{1}+\gamma_{2}\delta_{3})a_{1}.

Since $a_{1}>0$ and due to (2), $\mbox{Det}({\bf J})<0$ , and hence the Jacobian matrix is invertible. Therefore, by the implicit function theorem, there exists a constant $\varepsilon_{0}$ and unique functions $\tilde{S}^{*}(\varepsilon),I_{1}^{*}(\varepsilon),I_{2}^{*}(\varepsilon)$ defined for $0<\varepsilon<\varepsilon_{0}$ such that ${\bf F}(\tilde{S}^{*}(\varepsilon),I_{1}^{*}(\varepsilon),I_{2}^{*}(\varepsilon);\varepsilon)=0$ . Furthermore, these functions can be represented by a series expansion of the form (10).

∎

Next, we continue to focus on the case $\sigma_{21}>0$ and prove that $\phi^{*}$ is stable,

Proposition C.2 (Stability of $\phi^{*}$ for $\sigma_{21}>0$ ).

Let $\sigma_{21}>0$ , and let $\phi^{*}=(S^{*},I_{1}^{*},I_{2}^{*},J_{1}^{*},J_{2}^{*},R_{1}^{*},R_{2}^{*})$ be the steady-state solution of the system (1) with $I_{1}^{*}>0$ and $I_{2}^{*}>0$ . Then, for sufficiently large $\beta_{1}>\gamma_{1}$ , $\phi^{*}$ is locally stable.

Proof.

By Proposition C.1, $\phi^{*}$ satisfies (7). Substituting approximation (7) of the multi-strain endemic equilibrium point, $\phi^{*}$ , into the Jacobian matrix of (1) and computing the corresponding characteristic polynomial $P(\lambda)$ yields

	$\varepsilon^{2}P(\lambda)=a_{1}^{2}\sigma_{21}P_{0}(\lambda)+\varepsilon P_{1}(\lambda)+\varepsilon^{2}P_{2}(\lambda)+\varepsilon^{3}P_{3}(\lambda)+\mathcal{O}(\varepsilon^{4}),\qquad\varepsilon=\frac{1}{\mathcal{R}_{1}},\quad 0<\varepsilon\ll 1,$	(11a)
where
	$\begin{split}P_{0}(\lambda)&=(\lambda+\gamma_{1}+\mu)(\lambda+\gamma_{2}+\mu)Q(\lambda),\quad Q(\lambda)=\lambda^{3}+\alpha_{2}\lambda^{2}+\alpha_{1}\lambda+\alpha_{0},\\ P_{1}(\lambda)&=(\sigma_{21}+1)a_{1}\lambda^{6}+c_{1,5}\lambda^{6}+\cdots+c_{1,0},\\ P_{2}(\lambda)&=\lambda^{7}+c_{2,6}\lambda^{6}+\cdots+c_{2,0},\\ P_{3}(\lambda)&=c_{3,6}\lambda^{6}+\cdots+c_{3,0},\quad c_{3,6}\neq 0,\end{split}$	(11b)
and
	$\begin{split}\alpha_{2}&=b_{2}(\mu+\gamma_{2})\sigma_{12}+2\mu+\delta_{1}+\delta_{3}+\gamma_{1},\\ \alpha_{1}&=(2b_{2}\sigma_{12}+1)\mu^{2}+(b_{2}(\delta_{3}+\gamma_{1}+3\gamma_{2})\sigma_{12}+\delta_{1}+\delta_{3}+\gamma_{1})\mu+b_{2}\gamma_{2}(\delta_{3}+\gamma_{1}+\gamma_{2})\sigma_{12}+\delta_{3}(\delta_{1}+\gamma_{1})=\\ &b_{2}\sigma_{12}(2\mu^{2}+(\delta_{3}+\gamma_{1}+3\gamma_{2})\gamma_{2}+\mu(\delta_{3}+\gamma_{1}+\gamma_{2}))+\mu^{2}+(\delta_{1}+\delta_{3}+\gamma_{1})\mu+\delta_{3}(\delta_{1}+\gamma_{1})\\ \alpha_{0}&=b_{2}\sigma_{12}[\mu^{2}+(\delta_{3}+\gamma_{1}+\gamma_{2})\mu+(\delta_{3}+\gamma_{1})\gamma_{2}+\gamma_{1}\delta_{3})(\mu+\gamma_{2})],\\ c_{j,0}&\neq 0\quad j=1,2,3.\end{split}$	(11c)

Since $b_{2}>0$ , see (7c), and since $\delta_{i}>0$ or $\mu>0$ , it follows that the coefficients, $\alpha_{j}>0$ , $j=0,1,2$ , are all strictly positive.

The characteristic polynomial $P(\lambda)$ has overall seven roots. Five of the roots $\lambda_{i}$ of $P(\lambda)$ identify, to leading order, with the roots of $P_{0}(\lambda)$ . These roots satisfy $\lambda_{i}=\lambda_{i}^{\rm approx}+\mathcal{O}(\varepsilon)$ where

	$\lambda_{i}^{\rm approx}=-\mu-\gamma_{j}<0,\quad i=1,2,$	(12a)
and $\{\lambda_{i}^{\rm approx}\}_{i=3}^{5}$ are the roots of
	$Q(x)=x^{3}+\alpha_{2}x^{2}+\alpha_{1}x+\alpha_{0}=0,$	(12b)
see (11). The coefficients of $Q(x)$ are all strictly positive, see (11c), and direct computation shows that $\alpha_{2}\alpha 1>\alpha_{0}$ . Thus, the Routh–Hurwitz criterion implies that all roots of $Q(x)$ have a negative real part.

Dominant balance shows that $P(\lambda)$ has two additional roots that satisfy

\lambda_{i}^{\rm approx}=\frac{x_{i}}{\varepsilon},\quad\frac{|\lambda_{i}^{\rm approx}-\lambda_{i}|}{|\lambda_{i}|}=\mathcal{O(\varepsilon)},\qquad i=6,7,

(12c)

where $x_{i}$ are the roots of

x^{2}+(\sigma_{21}+1)a_{1}x+a_{1}^{2}\sigma_{21}=(x+\sigma_{21}a_{1})(x+a_{1}).

Hence,

x_{6}=-\sigma_{21}a_{1}<0,\quad x_{7}=-a_{1}<0.

Overall, all roots have a negative real part for sufficiently large $\beta_{1}>\gamma_{1}$ . Therefore, $\phi^{*}$ is linearly stable. ∎

Finally, we consider the case $\sigma_{21}=0$ . The proofs of the results, in this case, follow closely the proofs of Propositions C.1 and C.2.

Proposition C.3 (Existence, uniqueness and approximation of $\phi^{*}$ for $\sigma_{21}=0$ ).

Let $\sigma_{21}=0$ . Under the conditions of Theorem 2.2, if $\sigma_{12}>0$ , then there exists a unique steady-state solution

\phi^{*}=(S^{*},I_{1}^{*},I_{2}^{*},J_{1}^{*},J_{2}^{*},R_{1}^{*},R_{2}^{*})

of the system (1) with $I_{1}^{*}>0$ and $I_{2}^{*}>0$ . This solution satisfies (7).

Proof.

By Lemma C.1, $\phi^{*}$ satisfies (9). In the case, $\sigma_{21}=0$ , it follows that $J_{1}^{*}=0$ see (1d) and thus $S^{*}=1/\mathcal{R}_{1}$ when $I_{1}^{*}>0$ , see (1b). Denote $\varepsilon=1/\mathcal{R}_{1}$ and consider an expansion of the form

S^{*}=\varepsilon,\quad I^{*}_{i}=a_{i}+\varepsilon\tilde{b}_{i}+\mathcal{O}(\varepsilon^{2}),\qquad\sigma_{21}=0.

(13)

Substituting (13) in (9c), yields at leading order $a_{2}=0$ , and at the next order

a_{1}=\frac{(\gamma_{2}+\mu)\sigma_{12}b_{2}+(\delta_{1}+\mu)}{\mathcal{R}_{2}\eta_{2}\gamma_{1}},

as in the case of $\sigma_{21}>0$ . Further substituting (13) into (9a) yields (7g). As in the case of $\sigma_{21}>0$ , condition (6) ensures $I_{1}^{*},I_{2}^{*}>0$ .

The above computation uniquely determines the values of the coefficients in (7). Following the proof of Proposition C.1, we now use the implicit function theorem to show that there exists a unique solution to (9a), which can be expressed as an expansion (13) up to arbitrary order, by showing that the relevant Jacobian matrix is invertible: In the case $\sigma_{21}=0$ , $S(\varepsilon)=\varepsilon$ . Therefore, the Jacobian matrix ${\bf J}$ of (9a) at $\varepsilon=0$ reduces to an invertible $2\times 2$ matrix

\left.{\bf J}\right|_{(I_{1}^{*}(0),I_{2}^{*}(0))}=\left[\begin{array}[]{cc}-(\gamma_{1}+\mu)-\delta_{3}\frac{\gamma_{1}+\gamma_{2}+\mu}{\gamma_{2}+\mu}&\delta_{3}\frac{-(\gamma_{1}+\gamma_{2}+\mu)\delta_{3}+\gamma_{2}\delta_{2}}{\delta_{2}+\mu}\\ 0&-(\gamma_{2}+\mu)\end{array}\right].

∎

Proposition C.4 (Stability of $\phi^{*}$ for $\sigma_{21}>0$ ).

Let $\sigma_{21}=0$ , and let $\phi^{*}=(S^{*},I_{1}^{*},I_{2}^{*},J_{1}^{*},J_{2}^{*},R_{1}^{*},R_{2}^{*})$ be the steady-state solution of the system (1) with $I_{1}^{*}>0$ and $I_{2}^{*}>0$ . Then, for sufficiently large $\beta_{1}>\gamma_{1}$ , $\phi^{*}$ is locally stable.

Proof.

Substituting (7) in the Jacobian matrix of (1) and computing the corresponding characteristic polynomial $P(\lambda)$ yields

	$\varepsilon P(\lambda)=a_{1}(\gamma_{1}+\mu)P_{0}(\lambda)+\varepsilon P_{1}(\lambda)+\varepsilon^{2}P_{2}(\lambda)+\mathcal{O}(\varepsilon^{3}),\qquad\varepsilon=\frac{1}{\mathcal{R}_{1}},\quad 0<\varepsilon\ll 1,$	(14a)
where
	$\begin{split}P_{0}(\lambda)&=(\lambda+\gamma_{1}+\mu)(\lambda+\gamma_{2}+\mu)(\lambda+\delta_{2}+\mu)Q(\lambda),\quad Q(\lambda)=\lambda^{3}+\alpha_{2}\lambda^{2}+\alpha_{1}\lambda+\alpha_{0},\\ P_{1}(\lambda)&=\lambda^{7}+c_{1,6}\lambda^{6}+c_{1,5}\lambda^{5}+\cdots+c_{1,0},\\ P_{2}(\lambda)&=c_{2,6}\lambda^{6}+\cdots+c_{2,0},\end{split}$	(14b)
and
	$\begin{split}\alpha_{1}&=\gamma_{1}+\delta_{1}+\delta_{3}+2\mu+(\gamma_{2}+\mu)\sigma_{12}b_{2},\\ \alpha_{2}&=(\gamma_{2}+\mu)(\gamma_{1}+\gamma_{2}+\delta_{3}+2\mu)\sigma_{12}b_{2}+(\gamma_{1}+\delta_{1}+\mu)(\delta_{3}+\mu),\\ \alpha_{3}&=(\gamma_{2}+\mu)((\gamma_{1}+\delta_{3})\gamma_{2}+\gamma_{1}\delta_{3}+(\gamma_{1}+\gamma_{2}+\delta_{3})\mu+\mu^{2})\sigma_{12}b_{2}.\end{split}$	(14c)

The coefficients of $Q(x)$ are all strictly positive, and direct computation shows that $\alpha_{2}\alpha 1>\alpha_{0}$ . Thus, the Routh–Hurwitz criterion implies that all roots of $Q(x)$ have a negative real part.

Dominant balance and equation of orders shows that $P(\lambda)$ has an additional root

\lambda_{7}^{\rm approx}=-\frac{a_{1}(\gamma_{1}+\mu)}{\varepsilon}<0,\qquad\frac{|\lambda_{7}^{\rm approx}-\lambda_{7}|}{|\lambda_{7}|}=\mathcal{O(\varepsilon)}.

(15)

Overall, all roots have a negative real part for sufficiently large $\beta_{1}>\gamma_{1}$ . Therefore, $\phi^{*}$ is linearly stable ∎

Appendix D Numerical validation of Theorem 2.2

The asymptotic approximations presented in Theorem 2.2 are valid for sufficiently large values of $\mathcal{R}_{1}$ . In what follows we validate the asymptotic approximation (7) of Theorem 2.2, and show it is valid already for moderate values of $\mathcal{R}_{1}$ . To do so, we compute numerically $\phi^{*}$ for the values of the parameters used in the example presented in Figure 2, and then compute the error in the approximation (7). We note that we have considered additional sets of parameters and have reached similar results (data not shown).

Figure 5 presents the error in the approximation (7) for $\phi^{*}$ . As expected, we observe that $|S^{*}-\varepsilon-s_{2}\varepsilon^{2}|=O(\varepsilon^{3})$ and $|I_{i}^{*}-a_{i}-b_{i}\varepsilon|=O(\varepsilon^{2})$ . Formally the asymptotic results are valid for sufficiently small values of $\varepsilon$ . Yet, we observe that these results are valid also for $\varepsilon\approx 0.5$ , i.e., for $\mathcal{R}_{1}\approx 2$ .

References

[1] James D Bever, Thomas G Platt, and Elise R Morton. Microbial population and community dynamics on plant roots and their feedbacks on plant communities. Annual review of microbiology, 66:265–283, 2012.
[2] Hans J Bremermann and HR Thieme. A competitive exclusion principle for pathogen virulence. Journal of Mathematical Biology, 27(2):179–190, 1989.
[3] Li-Ming Cai, Maia Martcheva, and Xue-Zhi Li. Competitive exclusion in a vector–host epidemic model with distributed delay. Journal of biological dynamics, 7(sup1):47–67, 2013.
[4] MG Candela, E Serrano, C Martinez-Carrasco, P Martín-Atance, MJ Cubero, F Alonso, and L Leon. Coinfection is an important factor in epidemiological studies: the first serosurvey of the aoudad (ammotragus lervia). European journal of clinical microbiology & infectious diseases, 28:481–489, 2009.
[5] Carlos Castillo-Chavez, Herbert W Hethcote, Viggo Andreasen, Simon A Levin, and Wei M Liu. Epidemiological models with age structure, proportionate mixing, and cross-immunity. Journal of Mathematical Biology, 27:233–258, 1989.
[6] KW Chung and Roger Lui. Dynamics of two-strain influenza model with cross-immunity and no quarantine class. Journal of Mathematical Biology, 73:1467–1489, 2016.
[7] Yan-Xia Dang, Xue-Zhi Li, and Maia Martcheva. Competitive exclusion in a multi-strain immuno-epidemiological influenza model with environmental transmission. Journal of biological dynamics, 10(1):416–456, 2016.
[8] Odo Diekmann, Johan Andre Peter Heesterbeek, and Johan AJ Metz. On the definition and the computation of the basic reproduction ratio in models for infectious diseases in heterogeneous populations. Journal of Mathematical Biology, 28(4):365–382, 1990.
[9] Zhilan Feng and Jorge X Velasco-Hernández. Competitive exclusion in a vector-host model for the dengue fever. Journal of Mathematical Biology, 35:523–544, 1997.
[10] Georgy F Gause. Experimental analysis of Vito Volterra’s mathematical theory of the struggle for existence. Science, 79(2036):16–17, 1934.
[11] Nir Gavish and Guy Katriel. Optimal vaccination at high reproductive numbers: sharp transitions and counterintuitive allocations. Proceedings of the Royal Society B, 289(1983):20221525, 2022.
[12] Nir Gavish, Rami Yaari, Amit Huppert, and Guy Katriel. Population-level implications of the israeli booster campaign to curtail covid-19 resurgence. Science translational medicine, 14(647):eabn9836, 2022.
[13] Yair Goldberg, Micha Mandel, Yinon M. Bar-On, Omri Bodenheimer, Laurence S. Freedman, Nachman Ash, Sharon Alroy-Preis, Amit Huppert, and Ron Milo. Protection and waning of natural and hybrid immunity to sars-cov-2. New England Journal of Medicine, 386(23):2201–2212, 2022.
[14] Garrett Hardin. The competitive exclusion principle. Science, 131(3409):1292–1297, 1960.
[15] Matt J. Keeling and Pejman Rohani. Modeling Infectious Diseases in Humans and Animals. Princeton University Press, 2008.
[16] Md Abdul Kuddus, Emma S McBryde, Adeshina I Adekunle, and Michael T Meehan. Analysis and simulation of a two-strain disease model with nonlinear incidence. Chaos, Solitons & Fractals, 155:111637, 2022.
[17] Simon A Levin. Community equilibria and stability, and an extension of the competitive exclusion principle. The American Naturalist, 104(939):413–423, 1970.
[18] A Margheri and C Rebelo. Some examples of persistence in epidemiological models. Journal of mathematical biology, 46(6):564, 2003.
[19] Maia Martcheva. An introduction to mathematical epidemiology, volume 61. Springer, New York, NY, 2015.
[20] Maia Martcheva, Benjamin M Bolker, and Robert D Holt. Vaccine-induced pathogen strain replacement: what are the mechanisms? Journal of the Royal Society Interface, 5(18):3–13, 2008.
[21] Maia Martcheva and Xue-Zhi Li. Competitive exclusion in an infection-age structured model with environmental transmission. Journal of Mathematical Analysis and Applications, 408(1):225–246, 2013.
[22] Maia Martcheva and Xue-Zhi Li. Linking immunological and epidemiological dynamics of hiv: the case of super-infection. Journal of biological dynamics, 7(1):161–182, 2013.
[23] Martin A Nowak and Robert Mccredie May. Superinfection and the evolution of parasite virulence. Proceedings of the Royal Society of London. Series B: Biological Sciences, 255(1342):81–89, 1994.
[24] M Nuno, Carlos Castillo-Chavez, Z Feng, and M Martcheva. Mathematical models of influenza: the role of cross-immunity, quarantine and age-structure. Mathematical Epidemiology, pages 349–364, 2008.
[25] Miriam Nuño, Zhilan Feng, Maia Martcheva, and Carlos Castillo-Chavez. Dynamics of two-strain influenza with isolation and partial cross-immunity. SIAM Journal on Applied Mathematics, 65(3):964–982, 2005.
[26] Horst R Thieme. Pathogen competition and coexistence and the evolution of virulence. Mathematics for Life Science and Medicine, pages 123–153, 2007.
[27] Pauline Van den Driessche and James Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences, 180(1-2):29–48, 2002.
[28] Vito Volterra. Variazioni e fluttuazioni del numero d’individui in specie animali conviventi, volume 2. Societá anonima tipografica ”Leonardo da Vinci”, Cittá di Castello, 1927.
[29] Xiaoyan Wang, Junyuan Yang, and Xiaofeng Luo. Competitive exclusion and coexistence phenomena of a two-strain sis model on complex networks from global perspectives. Journal of Applied Mathematics and Computing, 68(6):4415–4433, 2022.
[30] LJ White, MJ Cox, and GF Medley. Cross immunity and vaccination against multiple microparasite strains. Mathematical Medicine and Biology: A Journal of the IMA, 15(3):211–233, 1998.
[31] Xiahong Zhao, Yilin Ning, Mark I-Cheng Chen, and Alex R Cook. Individual and population trajectories of influenza antibody titers over multiple seasons in a tropical country. American Journal of Epidemiology, 187(1):135–143, 2018.
[32] Tingting Zheng, Lin-Fei Nie, Zhidong Teng, and Yantao Luo. Competitive exclusion in a multi-strain malaria transmission model with incubation period. Chaos, Solitons & Fractals, 131:109545, 2020.

Revisiting the exclusion principle in epidemiology at its ultimate limit

Abstract

1 Introduction

2 Methods

2.1 Transmission model

2.1.1 Assumptions on model parameters

2.1.2 Basic reproduction number

2.2 Analysis in the case of a highly transmissible strain

Theorem 2.1.

Proof.

Theorem 2.2.

Proof.

3 Results

4 Conclusions

Acknowledgements

Appendix A Computation of ℛ0\mathcal{R}_{0}

Appendix B Proof of Theorem 2.1

Appendix C Proof of Theorem 2.2

Lemma C.1.

Proof.

Proposition C.1 (Existence, uniqueness and approximation of ϕ∗\phi^{*} for σ21>0\sigma_{21}>0).

Proof.

Proposition C.2 (Stability of ϕ∗\phi^{*} for σ21>0\sigma_{21}>0).

Proof.

Proposition C.3 (Existence, uniqueness and approximation of ϕ∗\phi^{*} for σ21=0\sigma_{21}=0).

Proof.

Proposition C.4 (Stability of ϕ∗\phi^{*} for σ21>0\sigma_{21}>0).

Proof.

Appendix D Numerical validation of Theorem 2.2

References

Appendix A Computation of $\mathcal{R}_{0}$

Proposition C.1 (Existence, uniqueness and approximation of $\phi^{*}$ for $\sigma_{21}>0$ ).

Proposition C.2 (Stability of $\phi^{*}$ for $\sigma_{21}>0$ ).

Proposition C.3 (Existence, uniqueness and approximation of $\phi^{*}$ for $\sigma_{21}=0$ ).

Proposition C.4 (Stability of $\phi^{*}$ for $\sigma_{21}>0$ ).