1. Introduction {#sec1}
===============

Infectious diseases cause huge morbidity and mortality in human populations and remain as a dramatic threat to public health \[[@B1], [@B2]\]. The incidence of many infectious diseases varies periodically. Cyclical infections of humans range from childhood infections such as measles, pertussis, mumps, diphtheria, varicella, and scarlet fever \[[@B3]\] to faecal-oral infections such as cholera and rotavirus and vector-borne diseases such as malaria and dengue fever \[[@B4], [@B5]\]. For example, seasonal flu develops as an epidemic during winter in temperate regions but it remains at very low levels during summer \[[@B6]\] while pertussis has 3-4-year cycles \[[@B7]\]. Despite the pervasive nature of oscillation, its underlying mechanisms are not well understood. From detailed analyses of empirical data, several distinct mechanisms have been proposed \[[@B4], [@B6], [@B8]\]: survival of disease pathogen outside host; host behaviour; seasonal changes in host susceptibility and immune defence and in vector abundance. Overall, those mechanisms have accounted only for the variation caused by external factors. Mechanically, they can be explained by a standard Susceptible-Infected-Recovered (SIR) or Susceptible-Infected-Recovered-Susceptible (SIRS) transmission model that was embedded with a seasonal transmission rate (i.e., the so-called "seasonal forcing") \[[@B4], [@B5], [@B8]--[@B12]\]. However, these models cannot provide a reasonable explanation for, say, a quickly globally spreading pattern of influenza \[[@B13]\] and require unrealistically high reproduction rates for the observational seasonal patterns in influenza \[[@B14]\]. Moreover, Grassly et al. \[[@B15]\] show that it is the intrinsic factors (e.g., immunity following the recovery from syphilis infection), rather than the external factors (e.g., changes in human sexual behaviour), that cause the 8--11-year cycles in syphilis incidence from the 1960s to the 1980s in the United States cities. Other mechanisms were also sought \[[@B16]\]. A single-strain standard SIR model that was not embedded with explicit seasonal forcing cannot generate sustained oscillations in incidence unless other complications are included. For example, variants of the standard SIR model such as cyclic models of SIRS and SEIRS type with a large time delay in the removed class, models with nonlinear force of infection, models with variable population size and disease-related deaths \[[@B17]\], or models including an isolated class generate sustained oscillations albeit under unrealistic region of parameter space \[[@B18]\]. However, it has been recognised that, without consideration of the intrinsic aspects of infectious diseases such as interactions between strains of an infectious agent, the understanding of epidemic cycling cannot be complete \[[@B10], [@B19]\] or is impossible \[[@B20]\].

Due to the high mutation rate and other genetic processes (such as recombination and reassortment), pathogens can quickly generate a lot of different genotypes. The striking fact of infectious agents is their antigenically diverse strains that are cocirculating within host populations \[[@B21]--[@B26]\] and concurrent infections with multiple strains \[[@B27]--[@B31]\]. Polymorphic pathogens such as dengue virus \[[@B22]\], malaria parasite*Plasmodium falciparum* \[[@B21]\], influenza virus \[[@B23], [@B24]\], and*Streptococcus pneumoniae* \[[@B26]\] are responsible for enormous health burdens of human populations \[[@B26], [@B32], [@B33]\]. Understanding the underlying mechanisms of antigenic diversity and cyclical epidemics will provide theoretical bases for their evolution under selective pressure such as the replacement of strains and the emergence of resistance strains under pharmacological control measures \[[@B25]\]. In practice, this will also be supplied with the basis for the design of effective control strategies such as vaccination \[[@B34]\].

In general, once two different pathogen strains infect the same host individuals, they will interact with each other \[[@B32]\]. The interaction can be classified into direct interactions of pathogen genes or gene products (such as recombination in bacterium*Streptococcus* and reassortment in influenza A viruses) and indirect interactions that result from alteration in the host environment such as immunological responses \[[@B35]\]. The indirect interactions may be synergistic, neutral, or antagonistic in view of changes in transmissibility and virulence \[[@B2], [@B33], [@B35]--[@B37]\]. For example, epidemiological evidence suggests that primary infection with one strain of dengue virus can enhance transmission of subsequent infection with another strain \[[@B38]--[@B40]\]. Further, the strains present in coinfections can transmit simultaneously and separately \[[@B29], [@B41], [@B42]\]. The direct observation of transmission process from one to another person is hard even for single transmission, let alone simultaneous transmission of multiple strains. To our knowledge, the one direct evidence of simultaneous transmission (or cotransmission) came from the field study during 2009 pandemic influenza by Liu et al. \[[@B29]\]. Cross-immunity between different subtypes of influenza has been observed \[[@B43]--[@B46]\]. Strong cross-immunity has been observed among antigenic drift variants of the same influenza subtype; for example, the first wave of the 1918 (H1N1) pandemic provided up to 94% protection against clinical illness during the second wave \[[@B47]\]. Cross-immunity was also found between different strains of other infectious agents \[[@B48], [@B49]\]. These characteristics surely affect pathogen diversity and dynamical behaviour of polymorphic pathogens. Though a large body of epidemiological researches are devoted to understanding the mechanisms of host single-strain pathogen dynamics, the increased evidence has gathered that the characteristics of polymorphic pathogen causing infections prove deviant from the predictions from one-strain models. For example, estimation of characteristics of transmission dynamics of polymorphic malaria*Plasmodium falciparum* that is based on one-strain epidemic model is misleading \[[@B50]\]; further, coinfections may imply that the survivorship of strains is not necessarily maximized by maximizing the basic reproductive number \[[@B51]\].

Transmission consequences of coinfection and superinfection have received a great deal of interest over the last 30 years \[[@B37], [@B52]--[@B58]\]. Theoretical modelling and empirical analyses show that pathogen interactions could have profound implications for the pathogen communities \[[@B51], [@B59]--[@B61]\]. For example, it has been shown that cocirculation of different diseases and cross-immunity \[[@B59]\] and antibody-dependent enhancement in transmission \[[@B38], [@B39], [@B62]\] can induce cyclical and chaotic epidemics. Cross-immunity alone within a two-strain SIR or SIRS model cannot generate sustained oscillation \[[@B12], [@B63]--[@B65]\] unless other factors such as age-dependent survival \[[@B53], [@B57]\] or isolation \[[@B66]\] are included. Following the traditional assumption that individuals, while being infected with one strain, are temporarily not susceptible to the other \[[@B52], [@B53], [@B62], [@B67]\], most of all these models consider only subsequent infection (i.e., individuals that recovered from primary infection with one strain were further infected by others). Even in those that included coinfections, different strains are not assumed to simultaneously transmit \[[@B12], [@B58], [@B61], [@B67], [@B68]\].

Recent empirical evidence has challenged the traditional assumption and indicated the concurrent infections during the infectious period of primary infection \[[@B29], [@B32], [@B41], [@B69], [@B70]\]. Due to technical difficulty and complexity involved in distinguishing multiple strains, multiple-strain infections are surely underestimated as argued by Balmer and Tanner \[[@B32]\]. The available empirical studies suggest that the occurrence of multiple strains was significantly more common than expected by chance. For example, Vaccarella et al. \[[@B70]\] found that the observed-to-expected ratio for concurrent infections with two human papillomavirus (HPV) types was 1.16 and for 3--7 types was 1.04. This implies that prior infection by one strain may facilitate infection by the other and strain interactions within coinfections can be synergism \[[@B42], [@B61]\]. For influenza that spreads from person-to-person through aerosols of droplets nuclei, for example, infection with one strain may increase the aerosolisation and hence the transmissibility of coinfecting other strains \[[@B71]\]. It is hence necessary to explore the possible role that strain interactions within concurrent infections might play in maintaining pathogen diversity and generating oscillatory diseases.

In this study we focus on the intrinsic factors of cyclical epidemics within a multistrain SIRS model by assuming that coinfections can directly transmit from people to people but at a low rate in relation to single transmission. We consider the indirect strain interactions that do not produce novel strains (cf. \[[@B59], [@B72]\]). In view of evidence of strong cross-immunity among strains (e.g., \[[@B47]\]), we assume a special situation where the patients recovered from infection will become completely cross-immune against any subsequent infections during their protective immunity periods. We thus exclude subsequent infection (cf. \[[@B38], [@B62], [@B67]\]). For concurrent infection of strains of a polymorphic pathogen, we have not found any informative data about the impact of order of infection and therefore ignore the impact of infection passage and superinfection \[[@B37], [@B42], [@B54]\]. We investigate how the interplay of complete cross-immunity during subsequent infection and strain interactions within coinfection affects antigenic diversity structure and generates cyclical epidemics.

2. Models and Methods {#sec2}
=====================

We consider a human population in which there are *n* cocirculating strains and every individual is capable of being simultaneously infected by at most *m*  (≤*n*) strains. For convenience, we denote the coinfection of multiplicity *l* as *i* ~1~ *i* ~2~ ⋯ *i* ~*l*~, here *i* ~1~, *i* ~2~,..., *i* ~*l*~ ∈ *H* ≡ {1,2, 3,..., *n*} and *i* ~1~ \< *i* ~2~ \< ⋯\<*i* ~*l*~ representing the strains present in the coinfection, which is independent of order of such infections. Let *I* ~*l*,*i*~1~⋯*i*~*l*~~ denote the fraction of individuals that carry the coinfection *i* ~1~ *i* ~2~ ⋯ *i* ~*l*~. In this study we define the multiplicity of strains in a coinfection as the within-in host antigenic diversity and the number (*n*) of cocirculating strains in the host population as the antigenic diversity of the population. We consider a SIRS epidemic model by assuming the coinfection *i* ~1~ *i* ~2~ ⋯ *i* ~*l*~ can simultaneously transmit among individuals at rate *β* ~*l*,*i*~1~⋯*i*~*l*~~; the coinfection simultaneously progresses to the recovered at rate *γ* ~*l*,*i*~1~⋯*i*~*l*~~ (cf. \[[@B34], [@B73]\]); strain interactions within coinfection occur via the change in infectivity (cf. \[[@B34]\]). Infections are acquired through a mass action term with the force of infection for the coinfection *i* ~1~ *i* ~2~ ⋯ *i* ~*l*~ given by $$\begin{matrix}
{\Lambda_{l,j_{1}\cdots j_{l}} = \beta_{l,j_{1}\cdots j_{l}}\left\lbrack {I_{l,j_{1}\cdots j_{l}} + {\sum\limits_{v = l + 1}^{m}{\phi_{v,l}{\sum_{\{ j_{1},\ldots,j_{l}\} \subseteq \{ k_{1},\ldots,k_{v}\}}I_{v,k_{1}\cdots k_{v}}}}}} \right\rbrack.} \\
\end{matrix}$$ Here *ϕ* ~*v*,*l*~ represents the factor that *l* strains transmit simultaneously from a patient coinfected with *v*  (\>*l*) strains, which reflects the strain interaction within coinfection. If strains interfere with each other, *ϕ* ~*v*,*l*~ \< 1; otherwise, *ϕ* ~*v*,*l*~ \> 1 implies enhancement in infectivity of strains. The factor *ϕ* ~*v*,*l*~ here is distinct from that proposed in dengue virus models (e.g., \[[@B38], [@B62]\]) where it reflects the change in transmissibility during secondary infection.

Further we assume that strain interaction also occurs via the change in susceptibility: once recovered from primary infection, patients become completely immune against any strain during the immunity period (1/*σ*), which is assumed to be the same for people recovered from different infections (cf. \[[@B51], [@B61]\]). This assumption neglects the complicated history of infection and thus largely simplifies the system (cf. \[[@B51], [@B63], [@B67], [@B74]\]). Further we assume that the birth to the susceptible and the loss of population by death are introduced through the same constant rate *μ*, therefore a constant population size. These assumptions lead to a system of coupled differential equations $$\begin{matrix}
{\frac{dS}{dt} = - {\sum\limits_{v = 1}^{m}{\sum_{\{ k_{1},\ldots,k_{v}\}}\Lambda_{v,k_{1}\cdots k_{v}}}}S + \mu\left( 1 - S \right)} \\
{\quad\quad   + \sigma\left( {1 - S - {\sum\limits_{l = 1}^{m}{\sum_{i_{1}\ldots i_{l} \in H}I_{l,i_{1}\cdots i_{l}}}}} \right),} \\
\frac{dI_{l,i_{1}\cdots i_{l}}}{dt} \\
{\quad = {\sum\limits_{v = 1}^{l}{\sum_{\{ j_{1},\ldots,j_{v}\} \cup \{ j_{v + 1},\ldots,j_{l}\} = \{ i_{1},\ldots,i_{l}\}}{\Lambda_{v,j_{1}\cdots j_{v}}I_{l - v,j_{v + 1}\cdots j_{l}}}}}} \\
{\quad\quad - \left( {\sum\limits_{v = 1}^{m - l}{\sum_{\{{j_{1},\ldots,j_{v}}\} \cap \{{i_{1}\cdots i_{l}}\} = Ø}\Lambda_{v,j_{1}\cdots j_{v}}}} \right)I_{l,i_{1}\cdots i_{l}}} \\
{\quad\quad - \left( {\gamma_{l,i_{1}\cdots i_{l}} + \mu} \right)I_{l,i_{1}\cdots i_{l}}.} \\
\end{matrix}$$ Here *Ø* is used to stand for empty set. In the above equations, *l* = 1,..., *m*, and *S*  ( = *I* ~0~) represents the fraction of the susceptible. For coinfection of multiplicity *l*, there are *n* ~*l*~ = *n*!/*l*!(*n* − *l*)! different combinations. A very large population is assumed to ignore the stochasticity.

The system at an arbitrary number of strains is intractable. We simplify the problem to investigate only the symmetric multistrain system where all strains share the same epidemiological properties: *β* ~1,*i*~ = *β* ~1~ and *γ* ~1,*i*~ = *γ* ~1~,  *i* = 1,..., *n* (cf. \[[@B51], [@B67]\]). Further, strain interactions are identical and depend only on the multiplicity of coinfections: *β* ~2,*ij*~ = *β* ~2~,..., *β* ~*m*,*i*~1~*i*~2~⋯*i*~*m*~~ = *β* ~*m*~,  *ϕ* ~*j*,*l*~ = *ϕ* ~*l*~ for *j* \> *l* and *γ* ~2,*ij*~ = *γ* ~2~,..., *γ* ~*m*,*i*~1~*i*~2~⋯*i*~*m*~~ = *γ* ~*m*~. At the equilibrium of such symmetric system the fraction of individuals that were coinfected with *m* strains is the same among different combinations; for example, *I* ~1,*i*~ ~1~ = *I* ~1~,..., *I* ~*m*,*i*~1~*i*~2~⋯*i*~*m*~~ = *I* ~*m*~, *i* ~1~,..., *i* ~*m*~ ∈ *H*. The model system reduces to that of 1 + *m* equations $$\begin{matrix}
{\frac{dS}{dt} = - {\sum\limits_{i = 1}^{m}{\begin{pmatrix}
n \\
i \\
\end{pmatrix}\Lambda_{i}}}S + \mu\left( {1 - S} \right) + \sigma\left( {1 - S - {\sum\limits_{l = 1}^{m}{n_{l}I_{l}}}} \right),} \\
{\frac{dI_{l}}{dt} = {\sum\limits_{i = 1}^{l}{\begin{pmatrix}
l \\
i \\
\end{pmatrix}\Lambda_{i}I_{l - i}}} - {\sum\limits_{j = 1}^{m - l}{\begin{pmatrix}
{n - l} \\
j \\
\end{pmatrix}\Lambda_{j}}}I_{l} - \left( {\gamma_{l} + \mu} \right)I_{l},} \\
{\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad l = 1,\ldots,m,} \\
\end{matrix}$$ with the force of infection, $$\begin{matrix}
{\Lambda_{l} = \beta_{l}\left\lbrack {I_{l} + \phi_{l}{\sum\limits_{j = l + 1}^{m}{\begin{pmatrix}
{n - l} \\
{j - l} \\
\end{pmatrix}I_{j}}}} \right\rbrack,} \\
{\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad l = 1,\ldots,m.} \\
\end{matrix}$$ The prevalence of coinfections of multiplicity *l* is *p* ~*l*~ = *n* ~*l*~ *I* ~*l*~ and the total prevalence is *p* = ∑~*l*=1~ ^*m*^ *p* ~*l*~. For further simplicity we assume that all the efficiencies of partial transmission from coinfections are the same (i.e., *ϕ* ~1~ = *ϕ* ~2~ = ⋯ = *ϕ* ~*m*~ = *ϕ*) and the infectious periods of all infections are the same (i.e., 1/*γ* ~1~ = 1/*γ* ~2~ = ⋯ = 1/*γ* ~*m*~ = 1/*γ*).

3. Results {#sec3}
==========

Though analytical results are hard to obtain for the general model system, they can be obtained for a two-strain model which was given in Appendix. Otherwise, our results are mainly based on the numerical integrations of ([3](#EEq3){ref-type="disp-formula"}) and ([4](#EEq4){ref-type="disp-formula"}). We first illustrate the invasion conditions for a pathogen strain to the host population in which other strains were established and then demonstrate how coinfections affect the antigenic diversity structure and facilitate the emergence of cyclical epidemics.

3.1. Dynamics of Invasion and Coexistence of Multiple Strains {#sec3.1}
-------------------------------------------------------------

Consider a situation where a host population was already at equilibrium with *n* ~*e*~ endemic strains, each of which has a basic reproductive number *R* ~0~ = *β* ~1~/(*μ* + *γ*). What is the condition for a new indistinguishable strain to invade into the population? In the case of *n* ~*e*~ = 1, as shown in Appendix, it is that the efficiency of partial transmission from coinfections (*ϕ*) and the simultaneous transmission rate of double infections (*β* ~2~) satisfy the inequality $$\begin{matrix}
{\phi \geq \phi_{s} \equiv 0.5 - \frac{\beta_{2}}{2\beta_{1}}.} \\
\end{matrix}$$ This approximated relationship holds for pathogens that have a much shorter infectious period than the duration of immunity such as influenza virus.

Though condition ([5](#EEq5){ref-type="disp-formula"}) was analytically obtained for the system of two strains, numerical simulations confirm it as a good approximation for multistrain systems (see [Table 1](#tab1){ref-type="table"}). Results show that when the efficiency of partial transmission from coinfections is less than the critical value (*ϕ* ~*s*~), only a single strain persists stably. Otherwise, a new indistinguishable strain can successfully invade and coexist with the strains in the previous endemic. That is, once condition ([5](#EEq5){ref-type="disp-formula"}) is satisfied, arbitrary number of strains can be maintained. [Table 1](#tab1){ref-type="table"} also indicates that simultaneous transmissions of triple and quadruple infections have a very weak impact on the invasion and pathogen diversity. This implies that, with simultaneous transmission of double infections, pathogen diversity can increase even under fairly strong interference between strains (i.e., *ϕ* ~*s*~ is much smaller than the unit) (cf. \[[@B73]\]). The critical value (*ϕ* ~*s*~) is 50% in the absence of simultaneous transmission of double infections; with simultaneous transmission (*β* ~2~), the critical value (*ϕ* ~*s*~) further reduces as predicted in ([5](#EEq5){ref-type="disp-formula"}) ([Table 1](#tab1){ref-type="table"}). Therefore, given that the basic reproductive number of each strain is larger than unit, equilibrium dynamics fall into two classes depending on the efficiency of partial transmission from coinfections and simultaneous transmission: stable single-strain persistence and multistrain persistence with stable diversity (cf. \[[@B67]\]). Hence it is the strain interaction within coinfection (i.e., *β* ~2~ and *ϕ*) that changes the traditional competitive exclusion \[[@B67], [@B73]\] and facilitates the increase in antigenic diversity.

3.2. Structure of Antigenic Diversity {#sec3.2}
-------------------------------------

Under the traditional circumstance where only one strain can transmit (i.e., *β* ~1~ \> 0 but *β* ~*l*~ = 0 for *l* \> 1), the antigenic diversity structure is shown in [Figure 1](#fig1){ref-type="fig"}. It shows that as the antigenic diversity (*l*) of coinfections increases, the proportion of coinfections of the same *l* strains (i.e., *I* ~*l*~) decrease as 1/*n* ~*l*~; however, their prevalence (i.e., *p* ~*l*~) remains roughly the same. The reduction in *I* ~*l*~ with multiplicity *l* is so remarkable that single infections are the most common infections (Figures [1(a)](#fig1){ref-type="fig"} and [1(b)](#fig1){ref-type="fig"}). When transmissibility is relatively low, triple and quadruple infections are prohibited (Figures [1(a)](#fig1){ref-type="fig"} and [1(c)](#fig1){ref-type="fig"}). However, with a high transmissibility, high diverse infections can safely survive (Figures [1(b)](#fig1){ref-type="fig"} and [1(d)](#fig1){ref-type="fig"}). Nevertheless, the total prevalence of all different forms of infections is nearly independent of the number of cocirculating strains and approximated as in the one-strain model system $$\begin{matrix}
{p = {\sum\limits_{l = 1}^{m}p_{l}} \approx \frac{\sigma + \mu}{\sigma + \mu + \gamma}\left( 1 - \frac{1}{R_{0}} \right)} \\
\end{matrix}$$ (see [Figure 1(e)](#fig1){ref-type="fig"}). This approximation can be simply confirmed from ([3](#EEq3){ref-type="disp-formula"}) and ([4](#EEq4){ref-type="disp-formula"}) under the situation where coinfections are rare. Ignoring coinfections, ([3](#EEq3){ref-type="disp-formula"}) and ([4](#EEq4){ref-type="disp-formula"}) reduce to that of one-strain SIR model with the prevalence of one strain replaced by total prevalence of *n* strains (cf. \[[@B67]\]). The derivation of approximation ([6](#EEq6){ref-type="disp-formula"}) for a special situation of two strains is given in Appendix.

The influences of simultaneous transmission of coinfections are shown in Figures [2](#fig2){ref-type="fig"} and [3](#fig3){ref-type="fig"} under two different immunity periods. Increasing *β* ~2~ increases the prevalence of coinfections but reduces single infections when strains interfere with each other (*ϕ* \< 1) (Figures [2(a)](#fig2){ref-type="fig"} and [3(a)](#fig3){ref-type="fig"}). However, if strains interact to enhance the infectivity within coinfections (*ϕ* \> 1), cyclical epidemics may emerge when immunity lasts longer ([Figure 2(b)](#fig2){ref-type="fig"}) or the prevalence of highly diverse coinfections increases substantially while single infections decrease greatly ([Figure 3(b)](#fig3){ref-type="fig"}). However, the impact of simultaneous transmission of more complicated coinfections is much weaker such that *β* ~4~ hardly have any effect (Figures [2(e)](#fig2){ref-type="fig"}, [2(f)](#fig2){ref-type="fig"}, [3(e)](#fig3){ref-type="fig"}, and [3(f)](#fig3){ref-type="fig"}). Under the circumstance of *β* ~2~ \> *β* ~3~ \> *β* ~4~, the total prevalence is dominated by single infections while the triple infections or higher complex coinfections increase but still remain at very low levels if the immunity period is long (cf. [Figure 1(d)](#fig1){ref-type="fig"}) and *ϕ* \< 1 (Figures [2(a)](#fig2){ref-type="fig"}, [2(c)](#fig2){ref-type="fig"}, [2(e)](#fig2){ref-type="fig"}, and [3(a)](#fig3){ref-type="fig"}). However, when strains facilitate each other to enhance their infectivity (*ϕ* \> 1), single infections cannot dominate (Figures [2(b)](#fig2){ref-type="fig"}, [2(d)](#fig2){ref-type="fig"}, and [2(f)](#fig2){ref-type="fig"}). Furthermore, in the situation of short immunity period, complicated infections such as double, triple, and quadruple infections become more popular than single infections (Figures [3(b)](#fig3){ref-type="fig"}, [3(d)](#fig3){ref-type="fig"}, and [3(f)](#fig3){ref-type="fig"}). This indicates that decreasing the duration of immunity can increase coinfection and stabilize the coexistence (cf. Figures [2(b)](#fig2){ref-type="fig"} and [3(b)](#fig3){ref-type="fig"}; see also [Figure 5(a)](#fig5){ref-type="fig"}). When the coexistence is steady, the total prevalence (*p*) slightly increases with these characteristics. Given all others being the same, increasing *β* ~2~ from 0 to (*β* ~1~ = )  1 only increases *p* less than 3% while increasing *β* ~3~ from 0 to 0.75 and increasing *β* ~4~ from 0 to 0.6 only increase *p* by 0.9% and 0.3%, respectively.

The interaction between strains is also reflected in the efficiency of partial transmission from coinfection (*ϕ*). When *ϕ* is not strong (≤1), single infections dominate ([Figure 1](#fig1){ref-type="fig"}). When *ϕ* is strong, given all others being the same, coinfections such as double and triple coinfections become common and predominate ([Figure 3](#fig3){ref-type="fig"}) and may also induce cyclical epidemics ([Figure 2(b)](#fig2){ref-type="fig"}). Nevertheless, as *ϕ* increases from 0 to 2, coinfections predominate but the total prevalence *p* only increases about 5%. Overall, the results cannot exceed 106% of ([6](#EEq6){ref-type="disp-formula"}) for all the steady scenarios examined in Figures [2](#fig2){ref-type="fig"} and [3](#fig3){ref-type="fig"}. These observations are also supported by the results listed in [Table 2](#tab2){ref-type="table"}. This indicates that ([6](#EEq6){ref-type="disp-formula"}) is a good approximation for the total prevalence of steady coexistence. Thus strain interactions (i.e., *β* ~2~ and *ϕ*) determine the structure of antigenic diversity but cannot much alter the total prevalence.

3.3. Cyclic Epidemics {#sec3.3}
---------------------

As shown above, coinfections facilitate coexistence; further they also aid the generation of sustained epidemic oscillations. Stability analysis and simulations show that with sufficiently enhanced transmissibility from coinfections the dynamic system will experience a Hopf bifurcation at which steady coexistences develop into limit cycles (Appendix). [Figure 4](#fig4){ref-type="fig"} shows the impact of simultaneous transmission of coinfections and antigenic diversity on the critical value of efficiency of partial transmission from coinfections (*ϕ* ~*c*~) for emergence of cyclical epidemics. In the absence of simultaneous transmission of coinfections, *ϕ* ~*c*~ \> 1. The observation under this special situation is similar to the current view that multistrain oscillations can occur only in the case of enhancement if cross-immunity during subsequent infection is very strong \[[@B38], [@B62], [@B75]\]. Two previous multistrain models \[[@B59], [@B63]\] which ignored current infection concluded that as cross-immunity becomes very strong, cyclical and chaotic epidemics vanish. In the presence of simultaneous transmission of double infections (*β* ~2~) and high antigenic diversity, *ϕ* ~*c*~ reduces to below unit. For example, when *β* ~2~ = *β* ~1~/2 = 0.5,  *ϕ* ~*c*~ = 0.81 for the situation of ten cocirculating strains. Comparatively, the simultaneous transmissions of triple and quadruple infections only have a much weak effect (Figures [4(b)](#fig4){ref-type="fig"} and [4(c)](#fig4){ref-type="fig"}). These indicate that *β* ~2~ has a predominant effect on the critical value *ϕ* ~*c*~. This implies that even with interference among strains, cyclic epidemics can still be generated provided there is enough simultaneous transmission of double infections and sufficient antigenic diversity.

[Figure 5](#fig5){ref-type="fig"} shows that cyclic epidemics are also controlled by other model parameters. For a given immunity period, [Figure 5(a)](#fig5){ref-type="fig"} shows that cyclic epidemics become impossible for a too small or too large reproductive number *R* ~0~. While for a given basic reproductive number, it also shows that the oscillation is prohibited under short immunity. Once emerged, the period of cyclic epidemics increases with immunity period and efficiency of partial transmission from coinfections but decreases with the basic reproductive number (Figures [5(b)](#fig5){ref-type="fig"} and [5(c)](#fig5){ref-type="fig"} and [Figure 6](#fig6){ref-type="fig"}). However, infectious period can only slightly increase the interepidemic period (data not shown). [Figure 6](#fig6){ref-type="fig"} shows some example time series of cyclical epidemics under different strain interactions. When strains interfere with each other (*ϕ* \< 1), single infections dominate the total prevalence; however, when they enhance each other to increase infectivity (*ϕ* \> 1), the fluctuations in prevalence become large (cf. [Figure 10](#fig10){ref-type="fig"}) and the fraction of complicated coinfections can become popular as seen in steady coexistence shown in [Figure 1(d)](#fig1){ref-type="fig"}.

4. Discussion {#sec4}
=============

In this study we investigate a multistrain SIRS epidemic model that includes coinfections and interactions between strains within coinfections. Though remaining at very low prevalence relative to single infections, coinfections and their simultaneous transmission substantially change the behaviour of infectious dynamics. They enable multiple strains to readily coexist and therefore maintain antigenic diversity. Further, interplay of strain interactions within coinfections and complete cross-immunity during subsequent infection induces sustained cyclical epidemics over a wide range of biologically realistic parameter space. In contrast to the heavily sought external seasonal forcing, this simple model provides a possible intrinsic mechanism for epidemic cycling of infectious diseases.

4.1. Invasion of New Strains and Antigenic Diversity {#sec4.1}
----------------------------------------------------

In contrast to other coinfection models \[[@B51]\], we include simultaneous transmission (or cotransmission) of multiple strains within concurrent infection. Even under total cross-immunity, coinfections increase the chance of invasion by a new strain into a population at endemic with other strains (Figures [7](#fig7){ref-type="fig"} and [8](#fig8){ref-type="fig"}; \[[@B40]\]) and simultaneous transmission enables a less transmissible strain to successfully invade ([Table 4](#tab4){ref-type="table"}). This allows multiple strains to coexist and increases the pathogen diversity though, in its absence, only one strain would persist, as suggested in traditional competitive exclusion ([Table 1](#tab1){ref-type="table"}; \[[@B67], [@B73]\]). This is comparable with multistrain models for dengue virus which include cross-reactive interaction during subsequent infection by other strains \[[@B76], [@B77]\]. Kooi et al. \[[@B76]\] show that though heterogeneity in transmission among strains reduces coexistence, the cross-reaction in second infection increases it. Interestingly, Mier-y-T-Teran-Romero et al. \[[@B77]\] show that combination of antibody-dependent enhancement and slight heterogeneity in transmission among strains increases persistence of four strains, suggesting evolutionary advantage of antibody-dependent enhancement. In contrast to Abu-Raddad et al. \[[@B51]\], our results show that, during the establishment process of a new strain, the total prevalence remains roughly unchanged. This observation suggests invisible invasion of genotypically different but antigenic indistinguishable new strains. When strains within coinfections can be transmitted simultaneously, coinfections can become prevalent if pathogen strains are highly transmissible (Figures [1](#fig1){ref-type="fig"}, [2](#fig2){ref-type="fig"}, and [3](#fig3){ref-type="fig"}). However, the long immunity period and interference between strains within coinfection will reduce the opportunities of high complicated coinfection and the antigenic diversity. This is comparable with predicted patterns of antigenic diversity by other models \[[@B25]\].

The observed levels of antigenic diversity may be large but surely are not infinite, say, influenza A \[[@B78]\]. So what limits the antigenic diversity? As the total prevalence is roughly constrained by ([6](#EEq6){ref-type="disp-formula"}) when multiple strains are in steady coexistence, the fraction of each single infection decreases proportionally with the antigenic diversity. In view of the actual limited size of host populations, stochastic extinction may be one of the forces that can limit the antigenic diversity \[[@B25], [@B67], [@B74]\]. Further, when infections experience oscillations, trough extinction ([Figure 6](#fig6){ref-type="fig"} and [Figure 12](#fig12){ref-type="fig"}; \[[@B65]\]) may be another mechanism. One important issue is that many newcomers, such as mutants, recombinants, reassortants, and immigrants, may have different epidemiological properties. Simulations show that any less transmissible newcomer will vanish no matter what the values of *β* ~2~ and *ϕ* are. The newcomers will replace endemic strains once they are more transmissible, which may also reduce pathogen diversity. Because of quick introduction of newcomers, the infections that were observed may never reach a static status. These factors may work together to restrain antigenic diversity.

4.2. Epidemic Cycling {#sec4.2}
---------------------

Though coinfection enhances coexistence of pathogen strains it alone cannot induce cyclical epidemics ([Figure 5(a)](#fig5){ref-type="fig"}). For example, Zhang et al. \[[@B36]\] considered a two-plant virus model in which the coinfection was included. No cyclical epidemic emerges because the infected plants die and there is no recovered compartment. The two-disease SIR models investigated by Blyuss and Kyrychko \[[@B79]\] and Martcheva and Pilyugin \[[@B57]\] included coinfections but neglected simultaneous transmission and cross-immunity. The sustained cycles were not generated in the former though induced in the latter, which is due to enhanced transmissibility of the strain that is weaker on its own.

Immunity is a common course of infection \[[@B2], [@B3]\] and is an important factor that leads to sustained oscillation \[[@B11], [@B15], [@B20]\]. Empirical analyses show that, for example, vaccination increases the degree of seasonal variation in the incidence of rotavirus \[[@B5]\]. However, theoretical investigations show that cross-immunity alone within a two-strain SIR model cannot support cyclical epidemics \[[@B63], [@B66], [@B80]\]. To simplify the study, we assume complete immunity and complete cross-immunity but the immunity decays. Therefore coinfection in this study is only referred to as concurrent infection (cf. \[[@B2], [@B38]\]). This is comparable with the assumption of transient strain-transcending immunity \[[@B74], [@B78]\] which states that patients who recovered from infection can be temporally protected against other strains and permanently against the previous infection strain. The decay of immunity in our model is assumed to reflect genotype changes in strains due to antigenic drift.

In some sense the immunity period can be regarded as the "strength" of the immunity: a short immunity period implying a weak immunity. It was demonstrated here that cyclical epidemics are prohibited under weak immunity ([Figure 5(a)](#fig5){ref-type="fig"}; \[[@B59]\]). [Figure 11](#fig11){ref-type="fig"} shows that pathogens that induce very short immunity can readily coexist and they remain steady if there is no immunity. Whilst those that induce lifelong immunity in host more likely lead to competitive exclusion, which suggests that extremely strong cross-immunity avoids cyclical epidemics (cf. \[[@B59], [@B63]\]). In the presence of simultaneous transmission, strain interactions within coinfections interplay with cross-immunity to generate cyclical epidemics (Figures [4](#fig4){ref-type="fig"} and [5](#fig5){ref-type="fig"}). The prevalent view of intrinsic mechanisms for cyclical epidemics in the situation of strong cross-immunity is the enhancement in sequential infections between strains \[[@B38], [@B59], [@B62], [@B63], [@B75]\]. In contrast with this, our results reveal that, with simultaneous transmissibility of strains and antigenic diversity, cyclical epidemics could be induced even if strains within coinfections interfere with each other. Furthermore, our results show that it is the simultaneous transmission of double infections that enables this, while the simultaneous transmissibility of more highly diverse coinfections only has minor effects (Figures [4](#fig4){ref-type="fig"} and [5](#fig5){ref-type="fig"}).

In our model, the interepidemic period is collectively determined by model parameters. Our analyses reveal that the combination of a short immunity period and a high *R* ~0~ gives rise to cyclical epidemics of short periods ([Figure 5](#fig5){ref-type="fig"}). This prediction is consistent with the results of Ferguson et al. \[[@B38]\] and is qualitatively in agreement with the empirically periodic patterns of childhood infection dynamics: seasonal flu that has a short immunity duration of 3--8 years \[[@B12], [@B81]\] and *R* ~0~ about 2 \[[@B82]\] was bound to annual cycles \[[@B24]\] while pertussis of a longer immunity period of about 30 years \[[@B83]\] leads to 3-4-year cycles \[[@B7]\]. Moreover, the interepidemic period initially increases with efficiencies of partial transmission from coinfections but quickly saturates ([Figure 5(c)](#fig5){ref-type="fig"} and [Figure 6](#fig6){ref-type="fig"}). It first increases with simultaneous transmission rate of double infections (*β* ~2~) and then decreases ([Figure 12](#fig12){ref-type="fig"}) while it is insensitive to simultaneous transmission of higher diverse coinfections when epidemic cycle first emerges ([Figure 4(d)](#fig4){ref-type="fig"}). Unfortunately, we currently do not have empirical data to test these relationships.

4.3. Comparison to Other Models of Epidemic Cycling {#sec4.3}
---------------------------------------------------

Oscillations in incidence were predicted in many previous epidemic models with different underlying mechanisms, for example, in one-strain SIR models \[[@B4], [@B8]--[@B11]\], two-strain SIR models \[[@B68]\], and two-strain SIRS models \[[@B12]\] where the seasonal forcing was thought to be the mechanism. Castillo-Chavez et al. \[[@B53]\] established that age-structure and cross-immunity generate sustained oscillations in a two-strain SIR model while Nuño et al. \[[@B66]\] showed that host isolation and cross-immunity lead to periodic epidemic outbreaks. Gupta et al. \[[@B59]\] investigated a multistrain SIR model in which coinfections included in the force of infection to reflect the generation of new strains through recombination are not set up as separate compartments and do not contribute to transmission. This is equivalent to the situation of *ϕ* = 0 and *β* ~*l*~ = 0 with *l* \> 1 within our model. Gupta et al. \[[@B59]\] showed that cyclical and chaotic epidemic is generated at intermediate cross-immunity but disappeared at strong cross-immunity. This is in agreement with our results shown in [Figure 4(a)](#fig4){ref-type="fig"}: cyclical and chaotic epidemics are prohibited when simultaneous transmission of and partial transmission from coinfections are absent. Ferguson et al. \[[@B38]\] and others \[[@B39], [@B62], [@B75]\] considered a similar SIR model and showed the enhancement of transmission during subsequent infection (i.e., so-called antibody-dependent enhancement) induces the cyclical and chaotic epidemics. Wearing and Rohani \[[@B84]\] and Bianco et al. \[[@B85]\] both consider multistrain models for dengue disease and combine antibody-dependent enhancement and temporal cross-immunity after recovering from primary infection. In contrast to Ferguson et al. \[[@B38]\] and other similar models, they can generate sustained oscillation in incidence even without antibody-dependent enhancement. However, the underlying mechanism for the cyclical epidemics is similar to what was shown by Nuño et al. \[[@B66]\]: the temporal cross-immunity class acts as the isolation class. All these dengue diseases models considered strain interaction during sequential infection, which is different from the concurrent infection we study here. Vasco et al. \[[@B61]\] considered a complicated variant of two-pathogen SIR model. It can produce sustained cyclical epidemics through a complicated combination of interactions between pathogens: competition due to cross-immunity, quarantine, or disease-induced mortality and cooperation due to immunosuppression and cross-enhancement. However, their model is different from ours in that the interaction between pathogens is associated with the changed susceptibility only and two pathogens cannot be simultaneously transmitted. Among these different mechanisms of epidemic cycling, which should be the truth? For example, inclusion of either age-structure \[[@B53], [@B57]\] or strain interaction within coinfection enables two-strain models to produce sustained oscillations in incidence. To discern which is correct, we need epidemiological data that include the additional information of age and antigenic diversity of infected people. These valuable data might be possible in near future, for example, for influenza viruses.

4.4. Simplifications and Possible Developments {#sec4.4}
----------------------------------------------

Our model is simplified in several aspects. For instance, the transmission rate is assumed to be constant, which is not biologically true. The transmissibility depends on infectiousness of disease agents, contact patterns, and susceptibility of host \[[@B1], [@B73], [@B86], [@B87]\] and these factors alter with environmental conditions. It has been shown that the inclusion of a stochastic environmental transmission in a one-strain SIR model can induce periodic outbreaks of infections \[[@B15], [@B88]\]. Actually, much effort has been made to search for the seasonal forcing that was induced from the variation in contact patterns and susceptibility of host populations \[[@B4], [@B8], [@B19]\]. These external mechanisms are certainly important, especially for explaining the consistent peak timing of cyclical epidemics, which appears as a weakness of the model presented here. However, they have difficulty in explaining the globally quick spread of influenza \[[@B13]\] and especially for the cyclical dynamics of syphilis \[[@B15]\]. It has been suggested that the mechanisms for the common oscillation of infectious diseases must be multiple \[[@B6], [@B10], [@B19]\]. The true mechanism might lie in the interplay of the intrinsic (as suggested by this investigation and Grassly et al. \[[@B15]\]) and external periodic forcing \[[@B4], [@B8]\].

For simplicity, we assume complete cross-immunity during subsequent infection. In reality, cross-immunity is more likely to be partial \[[@B43]--[@B45], [@B89]\] and therefore infection by different strains is likely to occur sequentially. Relaxing these assumptions will improve the predictions but will not change the essential findings: coinfections facilitate pathogen diversity, and simultaneous transmission from coinfections collaborates with cross-immunity to induce cyclical epidemics.

Xu-Sheng Zhang was funded by the Public Health England and Ke-Fei Cao was funded by the National Natural Science Foundation of China (NSFC) (Grant nos. 11365023 and 10565004) and the Specialized Research Fund for the Doctoral Program of Higher Education of China (SRFDP) (Grant no. 20050673001). The authors would like to thank the anonymous reviewers for their helpful comments and constructive suggestions on the paper.

A. Two-Strain SIRS Model {#secA}
========================

In this appendix we consider a special scenario with *n* = *m* = 2 but allow two strains to differ in epidemiological properties (i.e., so-called asymmetric two-strain model). From the model systems ([1](#EEq1){ref-type="disp-formula"}) and ([2](#EEq2){ref-type="disp-formula"}), the two-strain SIRS model is given by $$\begin{matrix}
{\frac{dS}{dt} = \mu - \left\lbrack {\beta_{a}\left( I_{a} + \phi_{a}I_{d} \right) + \beta_{b}\left( I_{b} + \phi_{b}I_{d} \right)} \right.} \\
{\quad\quad  \left. {+ \beta_{d}I_{d} + \mu} \right\rbrack S + \sigma\left( {1 - S - I_{a} - I_{b} - I_{d}} \right),} \\
{\frac{dI_{a}}{dt} = \beta_{a}\left( {I_{a} + \phi_{a}I_{d}} \right)S - \beta_{b}\left( {I_{b} + \phi_{b}I_{d}} \right)I_{a} - \left( {\gamma_{a} + \mu} \right)I_{a},} \\
{\frac{dI_{b}}{dt} = \beta_{b}\left( {I_{b} + \phi_{b}I_{d}} \right)S - \beta_{a}\left( {I_{a} + \phi_{a}I_{d}} \right)I_{b} - \left( {\gamma_{b} + \mu} \right)I_{b},} \\
{\frac{dI_{d}}{dt} = \beta_{b}\left( {I_{b} + \phi_{b}I_{d}} \right)I_{a} + \beta_{a}\left( {I_{a} + \phi_{a}I_{d}} \right)I_{b} + \beta_{d}I_{d}S} \\
{- \left( {\gamma_{d} + \mu} \right)I_{d},} \\
\end{matrix}$$ where *I* ~*a*~,  *I* ~*b*~, and *I* ~*d*~ are the proportions in infections with strain *a* alone, strain *b* alone, and both. *β* ~*a*~,  *β* ~*b*~, and *β* ~*d*~ are transmission rate of strain *a*, strain *b*, and both, respectively, and *γ* ~*a*~,  *γ* ~*b*~, and *γ* ~*d*~ are their corresponding recovery rates.

There are four equilibria (see [Table 3](#tab3){ref-type="table"}); in the following we use stability analysis to establish conditions for successful invasion and coexistence and the condition for Hopf bifurcation from steady to cyclical epidemics.

A.1. Invasion Conditions and Stability Analysis of Equilibria {#secA.1}
-------------------------------------------------------------

The threshold condition for pathogen strain *a* to successfully invade a naive host population was that its basic reproductive number (*R* ~0~ ^*a*^) must be larger than one where *R* ~0~ ^*a*^ is defined as the average number of secondary infections that result from the introduction of a single infectious individual into an entirely susceptible population. Similarly, the critical condition for strain *b* to invade and persist in a host population in which strain *a* was established at equilibrium (*Q* ~*a*~\*) is its invasion (or effective) reproductive number being larger than unit; that is, *R* ~0~ ^*b*^(*Q* ~*a*~\*) \> 1 \[[@B57]\]. Applying the method of van den Driessche and Watmough \[[@B91]\], we obtain $$\begin{matrix}
{R_{0}^{b}\left( Q_{a}^{\ast} \right) = \frac{1}{2}\left\lbrack {\Phi + \sqrt{\Phi^{2} + 4\Delta}} \right\rbrack,} \\
\end{matrix}$$ where $$\begin{matrix}
{\Phi \equiv \frac{\beta_{b}}{R_{0}^{a}}\frac{1}{\gamma_{b} + \mu + \beta_{a}{\overset{-}{I}}_{a}} + \frac{\beta_{b}\phi_{b}{\overset{-}{I}}_{a} + \beta_{d}/R_{0}^{a}}{\gamma_{d} + \mu},} \\
\end{matrix}$$ $$\begin{matrix}
{\Delta \equiv \frac{\beta_{b}}{R_{0}^{a}}\frac{\phi_{b}\beta_{a}{\overset{-}{I}}_{a} - \beta_{d}/R_{0}^{a}}{\left( \gamma_{d} + \mu \right)\left( \gamma_{b} + \mu + \beta_{a}{\overset{-}{I}}_{a} \right)}.} \\
\end{matrix}$$ In expressions ([A.3](#EEqA.3){ref-type="disp-formula"})-([A.4](#EEqA.4){ref-type="disp-formula"}), ${\overset{-}{I}}_{a}$ is the fraction of infective individuals with strain *a* at equilibrium *Q* ~*a*~\*;  *R* ~0~ ^*d*^ represents the basic reproductive number of double infection. It is worth mentioning that *R* ~0~ ^*b*^(*Q* ~*a*~\*) = 1 when Φ + Δ = 1; thus, when Φ + Δ \> 1,  *R* ~0~ ^*b*^(*Q* ~*a*~\*) \> 1. To demonstrate how the interaction between strains determines the successful invasion and coexistence, we rearrange the inequality Φ + Δ \> 1 into $$\begin{matrix}
{\beta_{d} > \left\lbrack {\left( \gamma_{d} + \mu \right) - \phi_{b}\beta_{a}{\overset{-}{I}}_{a}\frac{1 + \Gamma}{1 - \Gamma}} \right\rbrack R_{0}^{a}} \\
\end{matrix}$$ with $\Gamma \equiv {({R_{0}^{b}/R_{0}^{a}})}{(1 + ({{{\overset{-}{I}}_{a}\beta_{a}}/{(\gamma_{b} + \mu)}}))}^{- 1}$. Under the situation where *R* ~0~ ^*a*^ = *R* ~0~ ^*b*^ and $\gamma_{b} + \mu \gg \beta_{a}{\overset{-}{I}}_{a}$ which equivalently requires that the infectious period is much shorter than the duration of immunity and life span, we have $\Gamma \approx 1 - {\overset{-}{I}}_{a}\beta_{a}/(\gamma_{b} + \mu)$, and the above inequality is approximated as follows: $$\begin{matrix}
{\beta_{d} > \left\lbrack {\left( \gamma_{d} + \mu \right) - 2\phi_{b}\left( \gamma_{b} + \mu \right)} \right\rbrack R_{0}^{a}.} \\
\end{matrix}$$

By the symmetry between strain *a* and strain *b*, *R* ~0~ ^*a*^(*Q* ~*b*~\*) can be obtained by swapping the indexes *a* and *b* in the above equations. When *R* ~0~ ^*a*^(*Q* ~*b*~\*) \> 1 and *R* ~0~ ^*b*^(*Q* ~*a*~\*) \> 1 hold simultaneously, two strains coexist. To illustrate these relationships, two examples are shown: one for how strain interaction affects the invasion and coexistence ([Figure 7](#fig7){ref-type="fig"}) and the other for distribution of four equilibria and the effect of simultaneous transmission ([Figure 8](#fig8){ref-type="fig"}).

The stability conditions for equilibria can be obtained by analysing the eigenvalues of Jacobian matrix of the system ([A.1](#EEqA.1){ref-type="disp-formula"}) at equilibria. Consider the stability condition of equilibrium *Q* ~*a*~\* of single infection with strain *a*; its characteristic equation can be reduced to $$\begin{matrix}
\left| {\mathbf{M} - \lambda\mathbf{I}} \right| \\
{\quad = \left| \begin{matrix}
{- \beta_{a}{\overset{-}{I}}_{a} - \mu - \sigma - \lambda} & {- \beta_{a}\overset{-}{S} - \sigma} \\
{\beta_{a}{\overset{-}{I}}_{a}} & {- \lambda} \\
\end{matrix} \right|} \\
{\quad\quad \times \left| \begin{matrix}
{\beta_{b}\overset{-}{S} - \beta_{a}{\overset{-}{I}}_{a} - \gamma_{b} - \mu - \lambda} & {\phi_{b}\beta_{b}\overset{-}{S}} \\
{\left( \beta_{a} + \beta_{b} \right){\overset{-}{I}}_{a}} & {\phi_{b}\beta_{b}{\overset{-}{I}}_{a} + \beta_{d}\overset{-}{S} - \gamma_{d} - \mu - \lambda} \\
\end{matrix} \right|} \\
{\quad \equiv \left| \mathbf{C} \right| \times \left| \mathbf{D} \right| = 0.} \\
\end{matrix}$$ Here $\overset{-}{S}$ and ${\overset{-}{I}}_{a}$ are given in [Table 3](#tab3){ref-type="table"}. Matrix**C** always has two negative eigenvalues while matrix**D** has negative eigenvalues only if the following inequalities hold: $$\begin{matrix}
{\beta_{b}\left( {\frac{1}{R_{0}^{b}} - \frac{1}{R_{0}^{a}}} \right) + \beta_{d}\left( {\frac{1}{R_{0}^{d}} - \frac{1}{R_{0}^{a}}} \right) > - \left( {\beta_{a} - \beta_{b}\phi_{b}} \right){\overset{-}{I}}_{a},} \\
\end{matrix}$$ $$\begin{matrix}
{1 > \frac{\beta_{b}\overset{-}{S}}{\gamma_{b} + \mu + \beta_{a}{\overset{-}{I}}_{a}} + \frac{\beta_{b}\phi_{b}{\overset{-}{I}}_{a} + \beta_{d}\overset{-}{S}}{\gamma_{d} + \mu}} \\
{+ \frac{\left( {\beta_{a}\phi_{b}{\overset{-}{I}}_{a} - \beta_{d}\overset{-}{S}} \right)}{\left( {\gamma_{b} + \mu + \beta_{a}{\overset{-}{I}}_{a}} \right)\left( {\gamma_{d} + \mu} \right)}\frac{\beta_{b}}{R_{0}^{a}} \equiv \Phi + \Delta.} \\
\end{matrix}$$ Inequality ([A.9](#EEqA.7){ref-type="disp-formula"}) is equivalent to *R* ~0~ ^*b*^(*Q* ~*a*~\*) \< 1 and it can also be expressed as follows: $$\begin{matrix}
{R_{0}^{a} > R_{0}^{b}\left( {1 - \frac{R_{0}^{d}}{R_{0}^{a}}} \right)\frac{\gamma_{b} + \mu}{\gamma_{b} + \mu + \beta_{a}{\overset{-}{I}}_{a}} + R_{0}^{d}} \\
{+ R_{0}^{1}\frac{\beta_{b}\phi_{b}{\overset{-}{I}}_{a}}{\gamma_{d} + \mu}\left( {1 + \frac{\gamma_{a} + \mu}{\gamma_{b} + \mu + \beta_{a}{\overset{-}{I}}_{a}}} \right).} \\
\end{matrix}$$ Combining this with inequality ([A.8](#EEqA.6){ref-type="disp-formula"}), we have *R* ~0~ ^*a*^ \> 1,  *R* ~0~ ^*a*^ \> *R* ~0~ ^*b*^, and *R* ~0~ ^*a*^ \> *R* ~0~ ^*d*^. These are stability conditions listed in [Table 3](#tab3){ref-type="table"}. Similarly we can obtain the stability conditions of equilibria *Q* ~*b*~\* and *Q* ~0~\*. Due to the complexity, the stability condition for coexistence equilibrium *Q* ~*ab*~\* is hard to obtain analytically though in principle it can be determined through the negativeness of the eigenvalues of the characteristic equation (see below). Interestingly, numerical simulations (data not shown) demonstrate that even when *R* ~0~ ^*a*^ \< 1 and *R* ~0~ ^*b*^ \< 1, conditions that *R* ~0~ ^*d*^ \> 1 plus *ϕ* ~*a*~ \> 0 and *ϕ* ~*b*~ \> 0 can guarantee the coexistence of both strains (cf. \[[@B91]\]).

Some special situations are listed in [Table 4](#tab4){ref-type="table"} to illustrate the effect of coinfections on the successful invasion. Without simultaneous transmission of coinfections, strain *b* can invade only when it has a higher basic reproductive number *R* ~0~ ^*b*^ (see also [Figure 8](#fig8){ref-type="fig"}). However, strain interactions within coinfections change this outcome. If *R* ~0~ ^*b*^ = *R* ~0~ ^*a*^, strain *b* can invade if the efficiency of transmitting strain *b* from coinfection (*ϕ* ~*b*~) is larger than (1/2)(1/*γ* ~*b*~)/(1/*γ* ~*d*~), suggesting that if coinfections have doubled infectious periods, strain *b* can invade if *ϕ* ~*b*~ \> 25%. On the other extreme scenario where *ϕ* ~*b*~ is nil, strain *b* can invade only when *R* ~0~ ^*d*^ \> *R* ~0~ ^*a*^ and the previous single infection will be replaced by coinfections. In particular, it is possible for a less transmissible strain *b* (i.e., *R* ~0~ ^*b*^ \< *R* ~0~ ^*a*^) to invade if both *β* ~*d*~ and *ϕ* ~*b*~ are strong enough.

Next we consider a symmetric situation where two strains share the same properties: *β* ≡ *β* ~*a*~ = *β* ~*b*~,  *γ* ≡ *γ* ~*a*~ = *γ* ~*b*~ = *γ* ~*d*~ and *ϕ* ≡ *ϕ* ~*a*~ = *ϕ* ~*b*~ so that from ([A.1](#EEqA.1){ref-type="disp-formula"}) we have at equilibrium ${\overset{-}{I}}_{1} = {\overset{-}{I}}_{2}$ and $$\begin{matrix}
{{\overset{-}{I}}_{1} + {\overset{-}{I}}_{2} + {\overset{-}{I}}_{d} = \frac{\sigma + \mu}{\sigma + \mu + \gamma}\left( {1 - \overset{-}{S}} \right).} \\
\end{matrix}$$ When *ϕ* is not very strong and *β* \> *β* ~*d*~ so that only the steady coexistence is possible, numerical results show that $\overset{-}{S} \gg {\overset{-}{I}}_{1} = {\overset{-}{I}}_{2} \gg I_{d}$ (Tables [1](#tab1){ref-type="table"} and [2](#tab2){ref-type="table"}); from the second and third equations of ([A.1](#EEqA.1){ref-type="disp-formula"}), we obtain the approximation $$\begin{matrix}
{\overset{-}{S} \approx \frac{\left( {\mu + \gamma} \right)}{\beta} = \frac{1}{R_{0}^{a}} = \frac{1}{R_{0}^{b}}.} \\
\end{matrix}$$ Putting ([A.12](#EEqA.9){ref-type="disp-formula"}) into ([A.11](#EEqA.8){ref-type="disp-formula"}) gives rise to approximation ([6](#EEq6){ref-type="disp-formula"}) in the text for the situation of two symmetric strains.

A.2. Analysis of Hopf Bifurcation {#secA.2}
---------------------------------

Because of the complexity in model system ([A.1](#EEqA.1){ref-type="disp-formula"}), the analytical expression for the condition of Hopf bifurcation is not obtained. Taking the advantage of R computing language (<http://www.r-project.org/>), however, we can numerically test and obtain the condition under which the Hopf bifurcation occurs. The characteristic equation at the coexistence equilibrium (*Q* ~*ab*~\*) is as follows: $$\begin{matrix}
{\left| {\mathbf{M} - \lambda\mathbf{I}} \right| = \left| \begin{matrix}
{m_{11} - \lambda} & m_{21} & m_{31} & m_{41} \\
m_{12} & {m_{22} - \lambda} & m_{32} & m_{42} \\
m_{13} & m_{23} & {m_{33} - \lambda} & m_{43} \\
m_{14} & m_{24} & m_{34} & {m_{44} - \lambda} \\
\end{matrix} \right| = 0,} \\
\end{matrix}$$ where elements *m* ~*ij*~ are functions of the interior equilibrium $(\overset{-}{S},{\overset{-}{I}}_{a},{\overset{-}{I}}_{b},{\overset{-}{I}}_{d},\overset{-}{R})$ and model parameters. The characteristic equation can be written as follows: $$\begin{matrix}
{\lambda^{4} + a_{3}\lambda^{3} + a_{2}\lambda^{2} + a_{1}\lambda + a_{0} = 0.} \\
\end{matrix}$$ The coefficients *a* ~*k*~,  *k* = 0,1,..., 3, are the functions of matrix elements *m* ~*ij*~,  *i*, *j* = 1,..., 4.

The quartic equation ([A.14](#EEqA.10){ref-type="disp-formula"}) can be solved by Ferrari\'s method. To investigate the Hopf bifurcation from a steady state to a periodic cycle in incidence, we examine the eigenvalues against the efficiency (*ϕ* ~*a*~, *ϕ* ~*b*~) of single transmission from the coinfection. In general, the four eigenvalues can be expressed as follows: $$\begin{matrix}
{\lambda_{j}\left( {\phi_{a},\phi_{b}} \right) = \eta_{j}\left( {\phi_{a},\phi_{b}} \right) + i\omega_{j}\left( {\phi_{a},\phi_{b}} \right),} \\
{\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad j = 1,\ldots,4.} \\
\end{matrix}$$ When the real parts *η* ~*j*~ are negative, the steady coexistence equilibrium is stable. However, when values of (*ϕ* ~*a*~, *ϕ* ~*b*~) increase to some critical values (*ϕ* ~*a*0~, *ϕ* ~*b*0~), real part *η* ~*j*~ of some eigenvalues will increase to zero and the equilibrium loses its stability and a Hopf bifurcation occurs \[[@B92]\], from which oscillation in incidence will be induced. The imaginary part of the dominant eigenvalue, $\hat{\omega}(\phi_{a},\phi_{b})$, is the frequency of the oscillation when the real part $\hat{\eta}(\phi_{a},\phi_{b})$ passes through zero to become positive and the interepidemic period can be calculated as follows: $$\begin{matrix}
{T = \frac{2\pi}{\hat{\omega}\left( \phi_{a},\phi_{b} \right)}.} \\
\end{matrix}$$ One numerical example was shown in [Figure 9](#fig9){ref-type="fig"}. As the cooperation between strains within coinfections becomes strong enough, complicated dynamical behaviour such as period-doubling bifurcation and chaotic epidemics can be generated ([Figure 10](#fig10){ref-type="fig"}). Similarly, this can also be tested against other model parameters.

Interestingly, the fraction of infections at the peak of a cyclical epidemic is much larger than that at equilibrium of a single strain established alone. As illustrated in [Figure 10(b)](#fig10){ref-type="fig"}, the peak fraction of single infections with dominant strain *a* is 0.24% while it is 0.04% at equilibrium if it alone infects the host population (from [Table 1](#tab1){ref-type="table"} and see also [Figure 10(a)](#fig10){ref-type="fig"}). Therefore we should be cautious when explaining the meanings and estimating values of the basic reproductive number (*R* ~0~) for cyclical infections. Moreover, the fraction of single infection with a slightly less transmissible strain is much lower while the fraction of coinfections remains extremely low even at the peak. In [Figure 10(b)](#fig10){ref-type="fig"}, they are 0.024% and 0.001%, respectively, which indicates the chance for coinfections to be detected being 1/264 ≈ 0.4% in relation to single infections. Therefore it is understandable that coinfections are practically hard to attract sensible attention, especially if advantageous technology such as genotype typing is not readily available.

The cyclical epidemics become hard to generate when immunity period becomes either too short or too long ([Figure 11](#fig11){ref-type="fig"}). Given low reproductive numbers for two different strains but relatively high transmissibility for double infection as in [Figure 11](#fig11){ref-type="fig"}, the critical efficiency of single transmission from double infection for cyclical epidemics is *ϕ* ~*a*0~ = *ϕ* ~*b*0~ = 4.0 for immunity that lasts for one month or *ϕ* ~*a*0~ = *ϕ* ~*b*0~ = 9.2 for lifelong immunity. This may be too high for most pathogens. However, when the duration of immunity takes some intermediate values, a low efficiency of single transmission from double infection should be sufficient to maintain the cyclical epidemics. When oscillation is induced, the interepidemic period increases with the efficiency (*ϕ*) of single transmission from coinfections (cf. Figures [10(b)](#fig10){ref-type="fig"} and [10(c)](#fig10){ref-type="fig"}). While it first increases with the simultaneous transmissibility (*R* ~0~ ^*d*^), after reaching peak at a critical value, it decreases ([Figure 12](#fig12){ref-type="fig"}).

Conflict of Interests
=====================

The authors declare that there is no conflict of interests regarding the publication of this paper.

![Influence of the number of strains, the immunity period, and the basic reproductive number on the antigenic diversity structure. Model parameters: transmission rate: *β* ~1~ = 1.0 day^−1^ and *β* ~*l*~ = 0,  *l* ≥ 2,  *ϕ* = 1, and life span (1/*μ*) = 70 years. In panels (a) and (c) infectious period = 2 days and *R* ~0~ = 2.0; in panels (b) and (d) infectious period = 20 days and *R* ~0~ = 20.0. In panels (a) and (b) the total prevalence overlaps with the prevalence of single infections, implying most infections are single. In panel (e) the total prevalence *p* is well approximated by ([6](#EEq6){ref-type="disp-formula"}). Note that in panel (d) the prevalence *p* ~4~ exceeds *p* ~3~ when immunity period is around one year. This is the consequence of the assumption made in the example that the most complex coinfections are of multiplicity 4 so quadruple infections will not be further infected to become other more complex coinfections and hence its prevalence accumulated. The same phenomena also occur in [Figure 3](#fig3){ref-type="fig"}.](BMRI2014-375862.001){#fig1}

![Influence of simultaneous transmission of coinfections and efficiency of partial transmission from coinfections on the prevalence of single infections and coinfections. The effects of each simultaneous transmission rate (*β* ~*l*~) of coinfections of multiplicity are shown: panels (a) and (b) *l* = 2; (c) and (d) *l* = 3; and (e) and (f) *l* = 4 while others remain fixed. Two efficiencies of partial transmission from coinfections are shown: *ϕ* = 0.5 (panels (a), (c), and (e)) and =2 (panels (b), (d), and (f)). Model parameter values: immunity period (1/*σ*) = 10 years, transmission rate of single infection (*β* ~1~) = 1.0 day^−1^, and infectious period (1/*γ*) = 20 days; hence *R* ~0~ = 20; life span (1/*μ*) = 70 years, *n* = 10, and *m* = 4. Note that within panel (b) cyclical epidemics occur when *β* ~2~ = 0, 0.25, 0.5 day^−1^ with interepidemic periods around 420 days. The maximum (with subscript *M*) and minimum (with subscript *m*) prevalence are shown.](BMRI2014-375862.002){#fig2}

![As in [Figure 2](#fig2){ref-type="fig"} except immunity period = 3 years. Under this short immunity period, no cyclical epidemics appear within the ranges of parameters investigated.](BMRI2014-375862.003){#fig3}

![Impact of simultaneous transmission of coinfections and antigenic diversity on the critical value of the efficiency of partial transmission from coinfections (*ϕ* ~*c*~) for cyclic epidemics. Other model parameter values are 1/*σ* = 3 years, *β* ~1~ = 1.0 day^−1^, and 1/*μ* = 70 years. In panel (a) two situations are shown: two and ten strains circulate within the human population. For *n* = 10, the solid blue line represents the results for *β* ~3~ = *β* ~4~ = 0 and the dashed red line for *β* ~3~ = 0.40 day^−1^ and *β* ~4~ = 0.32 day^−1^. The black line represents the threshold condition for the second strain to successfully invade and coexist with strain 1. Since empirical patterns of incidence time series are very irregular, distinguishing between chaotic and cyclical epidemics is not important so we lump them together as chaotic/cyclical epidemics and just mention them as cyclical epidemics in text. In panels (a) and (b) infectious period = 2 days and *R* ~0~ = 2.0 for single infection. In panels (b) and (c) different combinations of simultaneous transmission rates of coinfections are assumed. In panel (c) the basic reproductive number of single infection (*R* ~0~) varies by altering the infectious period. Panel (d) shows the interepidemic periods in years when they first emerge under the three different simultaneous transmissions of various coinfections in panel (c).](BMRI2014-375862.004){#fig4}

![Impact of the basic reproductive number and immunity period on the critical value of the efficiency of partial transmission from coinfections (*ϕ* ~*c*~) for the cyclic epidemics and their periods. Five immunity periods are illustrated: 1, 3, 10, 20, and 50 years in panels (a) and (b). Panel (a) shows the dynamical maps with CCE representing chaotic/cyclical epidemics and panel (b) shows the interepidemic period when *ϕ* = *ϕ* ~*c*~. In panel (c) the interepidemic periods in years are shown against *ϕ* and *R* ~0~. Other parameter values are 1/*σ* = 10 years, 1/*μ* = 70 years, *β* ~1~ = 1.0 day^−1^, *β* ~2~ = 0.5 day^−1^, *β* ~3~ = 0.40 day^−1^, *β* ~4~ = 0.32 day^−1^, *n* = 10, and *m* = 4. The basic reproductive number of single infection (*R* ~0~) varies by altering the infectious period (1/*γ*).](BMRI2014-375862.005){#fig5}

![Example time series for the proportion of infected individuals. Other model parameter values are *β* ~1~ = 1.0 day^−1^, *R* ~0~ = 15,  1/*μ* = 70 years, *n* = 10, and *m* = 4; simultaneous transmissions of various coinfections as *β* ~2~ = 0.5 day^−1^, *β* ~3~ = 0.40 day^−1^, and *β* ~4~ = 0.32 day^−1^. Note that the steady total prevalence from ([6](#EEq6){ref-type="disp-formula"}) would be 4.37 × 10^−3^ if the system remained in steady state. This is much smaller than the maximal prevalence within cyclical epidemics shown here.](BMRI2014-375862.006){#fig6}

![Impact of strain interactions on coexistence and sustained oscillation in incidence of infections. The host population was initially endemic with strain *a* and then challenged with strain *b*. Two strains are assumed to be equally transmissible with *β* ~*a*~ = *β* ~*b*~ = 0.75 day^−1^ and infectious period (1/*γ* ~*a*~ = 1/*γ* ~*b*~) = 2 days (i.e., *R* ~0~ ^*a*^ = *R* ~0~ ^*b*^ = 1.5). Coinfection can simultaneously transmit at a rate (*β* ~*d*~) = 0.375 day^−1^ and with an infectious period (1/*γ* ~*d*~) = 2 days. Other parameters are the immunity period (1/*σ*) = 3 years and life span (1/*μ*) = 55 years.](BMRI2014-375862.007){#fig7}

![Impact of efficiency of single transmission (*ϕ*) and simultaneous transmissibility (*R* ~0~ ^*d*^) from dual infection on the distribution of equilibria. The four equilibria on the (*R* ~0~ ^*a*^, *R* ~0~ ^*b*^) plane are separated by four curves: *R* ~0~ ^*a*^ = 1, *R* ~0~ ^*b*^ = 1, *R* ~0~ ^*b*^(*Q* ~*a*~\*) = 1, and *R* ~0~ ^*a*^(*Q* ~*b*~\*) = 1. In panel (a) where *R* ~0~ ^*d*^ = 0 and *ϕ* ≤ 0.5 or *R* ~0~ ^*d*^ = *R* ~0~ ^*a*^/2 and *ϕ* ≤ 0.25,  *R* ~0~ ^*b*^ \> *R* ~0~ ^*a*^ is required to maintain *R* ~0~ ^*b*^(*Q* ~*a*~\*) ≥ 1 (and similarly requiring *R* ~0~ ^*a*^ \> *R* ~0~ ^*b*^ to maintain *R* ~0~ ^*a*^(*Q* ~*b*~\*) ≥ 1); hence no coexistence is possible. In panels (b) and (c) we consider the symmetric situation in which *ϕ* ~*a*~ = *ϕ* ~*b*~ = *ϕ* = 5 under different simultaneous transmissibility. The other parameters are infectious period (1/*γ* ~*a*~ = 1/*γ* ~*b*~ = 1/*γ* ~*d*~) = 2 days, immunity period = 3 years, and life span = 55 years.](BMRI2014-375862.008){#fig8}

![Numerical integration of ([A.1](#EEqA.1){ref-type="disp-formula"}) with the efficiency of single transmission from dual infection *ϕ*  ( = *ϕ* ~*a*~ = *ϕ* ~*b*~) = 1.75, 1.76, and 1.77. The real parts of two identical eigenvalues of characteristic equation ([A.6](#EEqA.5){ref-type="disp-formula"}) increase with *ϕ* and passed through 0 at 1.76 while the real parts of other two eigenvalues decrease with *ϕ* and remain negative. The imaginary part of the dominant eigenvalue is 0.0154, *T* = 410 days from ([A.16](#EEqA.11){ref-type="disp-formula"}), in agreement with the result shown in the graph. As we assume both strains are equal, the fractions of both single infections overlap. Other parameters are transmission rates *β* ~*a*~ = *β* ~*b*~ = 0.75 day^−1^, *β* ~*d*~ = 0.375 day^−1^, infectious period (1/*γ* ~*a*~ = 1/*γ* ~*b*~ = 1/*γ* ~*d*~) = 2 days, immunity period = 3 years, and life span = 55 years.](BMRI2014-375862.009){#fig9}

![(a) Bifurcation diagram plotting the local maxima of *I* ~*a*~ (the fraction of the infective individuals with strain *a*) against efficiency (*ϕ* = *ϕ* ~*a*~ = *ϕ* ~*b*~) of single transmission from coinfection. Numerical integration of the equations is illustrated at different values of *ϕ* ~*a*~ = *ϕ* ~*b*~: (b) 1.5 (period-2), (c) 3.9 (period-4), and (d) 5.0 (chaotic) for *I* ~*a*~. In panel (b) the fraction of infective individuals with strain *b* and that of coinfections are also included. The other parameters are immunity period = 5 years, *β* ~*a*~ = 0.75 day^−1^, *β* ~*b*~ = 0.74 day^−1^, *β* ~*d*~ = 0.7 day^−1^, and infectious period (1/*γ* ~*a*~ = 1/*γ* ~*b*~ = 1/*γ* ~*d*~) = 2 days. The first 100 years were discarded. Note that as *ϕ* increases, the trough that the infection experiences at valley becomes deeper and deeper while the epidemic interval shortens from 200 to 100 days. Here the epidemic interval is defined relatively as the interval during which the fraction of infections is maintained at \>10% of its peak size. This suggests that pathogen strains will suffer the trough extinction at large values of efficiency (*ϕ*).](BMRI2014-375862.010){#fig10}

![Impact of immunity period on coexistence and sustained oscillation in incidence. The narrow gap between two lines represents the area of steady coexistence. Note that when the cross-immunity remains permanent (i.e., *σ* = 0), two strains can coexist under strong enhancement in transmission (i.e., 9.1 \< *ϕ* ~*a*~ = *ϕ* ~*b*~ \< 9.2) and cyclical epidemic can emerge when *ϕ* ~*a*~ = *ϕ* ~*b*~ \> 9.2 (data not shown). Other model parameters are *β* ~*a*~ = 0.75,  *β* ~*b*~ = 0.74,  *β* ~*d*~ = 0.70, infectious period (1/*γ* ~*a*~ = 1/*γ* ~*b*~ = 1/*γ* ~*d*~) = 2 days, and life span (1/*μ*) = 55 years.](BMRI2014-375862.011){#fig11}

![Impact of simultaneous transmissibility on the interepidemic period and amplitude of infectious transmission dynamics. Two situations of transmissibility of single pathogen infection are considered: panels (a) and (b) *R* ~0~ ^*a*^ = *R* ~0~ ^*b*^ = 3.0 and immunity period = 3 years, (c) and (d) *R* ~0~ ^*a*^ = *R* ~0~ ^*b*^ = 17 and immunity period = 50 years. In panels (b) and (d), solid lines stand for the maxima while dotted ones lined the minima of the fraction of total infections. Other parameters are the efficiency of single transmission from coinfection (*ϕ*) = 2 and the infectious period (1/*γ* ~*a*~ = 1/*γ* ~*b*~ = 1/*γ* ~*d*~) in (a)-(b) which mimics influenza is 2 days while in (c)-(d) which mimics pertussis is 30 days. The first 100 years were discarded.](BMRI2014-375862.012){#fig12}

###### 

The dependence of the successful invasion by an indistinguishable new strain on simultaneous transmission rates of coinfections of multiplicity *m* = 2,3, 4 and the efficiency of partial transmission from coinfections (*ϕ*). We consider the condition where a new indistinguishable strain can invade a host population at equilibrium with *n* ~*e*~ endemic strains. After the challenge with a new identical strain, a new endemic will be established if condition ([5](#EEq5){ref-type="disp-formula"}) holds; otherwise, the invader is not successful and the old endemic crashes into one with a single strain. When the invasion is successful, the new endemic is established "invisibly" because the total prevalence, approximated by ([6](#EEq6){ref-type="disp-formula"}), remains unchanged. Note the prevalence of triple and quadruple infections is vanishingly small, implying their actual nonexistence. Other model parameter values are immunity period = 3 years, transmission rate of single infections (*β* ~1~) = 1.0, and infectious period = 2 days; hence *R* ~0~ = 2.0 for single infection.

  *n* ~*e*~   *β* ~2~          *β* ~3~          *ϕ*              *p*               *p* ~1~           *p* ~2~   *p* ~3~   *p* ~4~   *n* ~ec~\*
  ----------- ---------------- ---------------- ---------------- ----------------- ----------------- --------- --------- --------- ------------
  3           0.0              0.0              0.0--0.49        9.502 × 10^−4^    9.502 × 10^−4^    ---       ---       ---       1
  0.51        9.502 × 10^−4^   9.489 × 10^−4^   1.351 × 10^−6^   1.284 × 10^−9^    6.097 × 10^−13^   4                             
  0.5         0.0              0.0--0.24        9.502 × 10^−4^   9.502 × 10^−4^    ---               ---       ---       1         
  0.26        9.503 × 10^−4^   9.476 × 10^−4^   2.695 × 10^−6^   3.834 × 10^−9^    3.031 × 10^−12^   4                             
  0.4         0.0--0.24        9.502 × 10^−4^   9.502 × 10^−4^   ---               ---               ---       1                   
  0.26        9.503 × 10^−4^   9.476 × 10^−4^   2.697 × 10^−6^   6.398 × 10^−9^    5.463 × 10^−12^   4                             
  0.0--0.24   9.502 × 10^−4^   9.502 × 10^−4^   ---              ---               ---               1^&^                          
  0.26        9.503 × 10^−4^   9.476 × 10^−4^   2.697 × 10^−6^   6.400 × 10^−9^    9.109 × 10^−12^   4^&^                          
                                                                                                                                   
  2           0.0              0.0              0.0--0.49        9.502 × 10^−4^    9.502 × 10^−4^    ---       ---       ---       1
  0.51        9.506 × 10^−4^   9.494 × 10^−4^   1.201 × 10^−6^   7.612 × 10^−10^   ---               3                             
  0.5         0.0              0.0--0.24        9.502 × 10^−4^   9.502 × 10^−4^    ---               ---       ---       1         
  0.26        9.503 × 10^−4^   9.479 × 10^−4^   2.397 × 10^−6^   2.275 × 10^−9^    ---               3                             
  0.4         0.0--0.24        9.502 × 10^−4^   9.502 × 10^−4^   ---               ---               ---       1                   
  0.26        9.503 × 10^−4^   9.479 × 10^−4^   2.399 × 10^−6^   3.794 × 10^−9^    ---               3                             
                                                                                                                                   
  1           0.0              ---              0.0--0.50        9.502 × 10^−4^    9.502 × 10^−4^    ---       ---       ---       1
  ---         0.51             9.502 × 10^−4^   9.493 × 10^−4^   9.020 × 10^−7^    ---               ---       2                   
  ---         2                9.529 × 10^−4^   9.520 × 10^−4^   9.097 × 10^−7^    ---               ---       2                   
  0.5         ---              0.0--0.25        9.502 × 10^−4^   9.502 × 10^−4^    ---               ---       ---       1         
  ---         0.26             9.503 × 10^−4^   9.485 × 10^−4^   1.801 × 10^−7^    ---               ---       2                   

\**n* ~ec~ is the number of endemic strains after the challenge of the invader.

^&^In this situation we set *β* ~4~ = 0.4 while *β* ~4~ = 0 for all other situations.

###### 

Dependence of antigenic diversity structure in the steady coexistence on the transmission rates. Other parameter values are immunity period = 3 years, transmission rate of single infections (*β* ~1~) = 1.0, and infectious period = 2 days; hence *R* ~0~ = 2.0 for single infection; ten indistinguishable strains (*n* = 10) are assumed to cocirculate while three models are considered in which the most complex coinfections are of multiplicity *m* = 2,3, 4 strains, respectively.

  *m*    *β* ~2~   *β* ~3~   *β* ~4~          *ϕ*              *p*              *p* ~1~          *p* ~2~           *p* ~3~           *p* ~4~
  ------ --------- --------- ---------------- ---------------- ---------------- ---------------- ----------------- ----------------- -----------------
  4      0.5       0.0       0.0              0.75             9.535 × 10^−4^   9.502 × 10^−4^   3.263 × 10^−6^    7.469 × 10^−9^    1.664 × 10^−11^
  0.4    0.0       0.75      9.535 × 10^−4^   9.502 × 10^−4^   3.274 × 10^−6^   1.253 × 10^−8^   3.019 × 10^−11^                     
  0.4    0.4       0.75      9.535 × 10^−4^   9.502 × 10^−4^   3.274 × 10^−6^   1.256 × 10^−8^   5.038 × 10^−11^                     
  0.4    0.4       0.0^\$^   9.487 × 10^−4^   9.455 × 10^−4^   3.214 × 10^−6^   1.214 × 10^−8^   4.826 × 10^−11^                     
  0.75   0.0       0.0       0.38             9.535 × 10^−4^   9.470 × 10^−4^   6.444 × 10^−6^   1.714 × 10^−8^    6.170 × 10^−11^   
  0.6    0.0       0.38      9.536 × 10^−4^   9.470 × 10^−4^   6.529 × 10^−6^   4.349 × 10^−8^   1.328 × 10^−10^                     
  0.6    0.4       0.38      9.536 × 10^−4^   9.470 × 10^−4^   6.534 × 10^−6^   4.369 × 10^−8^   2.218 × 10^−10^                     
  0.6    0.4       0.0^\$^   9.487 × 10^−4^   9.423 × 10^−4^   6.386 × 10^−6^   4.209 × 10^−8^   2.117 × 10^−10^                     
                                                                                                                                     
  3      0.5       0.0       ---              0.75             9.535 × 10^−4^   9.502 × 10^−4^   3.263 × 10^−6^    7.478 × 10^−9^    ---
  0.4    ---       0.75      9.535 × 10^−4^   9.502 × 10^−4^   3.274 × 10^−6^   1.249 × 10^−8^   ---                                 
  0.4    ---       0.0^\$^   9.486 × 10^−4^   9.454 × 10^−4^   3.213 × 10^−6^   1.216 × 10^−8^   ---                                 
                                                                                                                                     
  2      0.5       ---       ---              0.75             9.535 × 10^−4^   9.502 × 10^−4^   3.256 × 10^−6^    ---               ---
  ---    ---       0.0^\$^   9.486 × 10^−4^   9.454 × 10^−4^   3.223 × 10^−6^   ---              ---                                 

^\$^These coexistence equilibria are unstable.

###### 

Four equilibria and their stability conditions of model system ([A.1](#EEqA.1){ref-type="disp-formula"}).

  ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  Equilibrium                        Solution                                                                                                                                                                                                                                          Stability conditions
  ---------------------------------- ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  ^&^Disease-free (*Q* ~0~\*)        $\overset{-}{S} = 1$,  ${\overset{-}{I}}_{a} = {\overset{-}{I}}_{b} = {\overset{-}{I}}_{d} = \overset{-}{R} = 0$                                                                                                                                  *R* ~0~ ^a^ \< 1, *R* ~0~ ^b^ \< 1, and *R* ~0~ ^d^ \< 1

                                                                                                                                                                                                                                                                                       

  ^\$^Strain *a* alone (*Q* ~a~\*)   $\overset{-}{S} = \frac{1}{R_{0}^{a}}$,  ${\overset{-}{I}}_{b} = {\overset{-}{I}}_{d} = 0$,\                                                                                                                                                      *R* ~0~ ^a^ \> 1, *R* ~0~ ^a^ \> *R* ~0~ ^b^, *R* ~0~ ^a^ \> *R* ~0~ ^d^, *R* ~0~ ^b^(*Q* ~a~\*) \< 1
                                     ${\overset{-}{I}}_{a} = \frac{\sigma + \mu}{\sigma + \mu + \gamma_{a}}\left( {1 - \frac{1}{R_{0}^{a}}} \right)$,  $\overset{-}{R} = \frac{\gamma_{a}}{\sigma + \mu + \gamma_{a}}\left( {1 - \frac{1}{R_{0}^{a}}} \right)$                         

                                                                                                                                                                                                                                                                                       

  Strain *b* alone (*Q* ~b~\*)       $\overset{-}{S} = \frac{1}{R_{0}^{b}}$,  ${\overset{-}{I}}_{a} = {\overset{-}{I}}_{d} = 0$, \                                                                                                                                                     *R* ~0~ ^b^ \> 1, *R* ~0~ ^b^ \> *R* ~0~ ^a^, *R* ~0~ ^b^ \> *R* ~0~ ^d^, *R* ~0~ ^a^(*Q* ~b~\*) \< 1
                                     ${\overset{-}{I}}_{b} = \frac{\sigma + \mu}{\sigma + \mu + \gamma_{b}}\left( {1 - \frac{1}{R_{0}^{b}}} \right)$,  $\overset{-}{R} = \frac{\gamma_{b}}{\sigma + \mu + \gamma_{b}}\left( {1 - \frac{1}{R_{0}^{b}}} \right)$                         

                                                                                                                                                                                                                                                                                       

  Coexistence (*Q* ~ab~\*)           Interior equilibrium $(\overset{-}{S},{\overset{-}{I}}_{a},{\overset{-}{I}}_{b},{\overset{-}{I}}_{d},\overset{-}{R})$ obtained numerically by solving equations ([A.1](#EEqA.1){ref-type="disp-formula"}) with the left sides being set to zero   *R* ~0~ ^a^(*Q* ~b~\*) \> 1, *R* ~0~ ^b^(*Q* ~a~\*) \> 1, and the real part of the dominant eigenvalue of Jacobian matrix of system ([A.1](#EEqA.1){ref-type="disp-formula"}) being negative (see [Appendix A.2](#secA.2){ref-type="sec"})
  ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

^&^ *R* ~0~ ^*i*^, *i* = *a*, *b*, *d*, are the basic reproduction number of strain *a*, strain *b*, and dual infection, respectively, and are given by *R* ~0~ ^i^ = β~i~/(γ~i~ + μ).

^\$^ *R* ~0~ ^b^(*Q* ~a~\*) represents the effective reproduction number of strain *b* upon the population in which strain *a* was established (*Q* ~a~\*) (see equation ([A.2](#EEqA.2){ref-type="disp-formula"})).

###### 

Impact of coinfection on invasion conditions.

  Contribution from coinfection to transmission                                                                                                                           Strain *b* can invade only if
  ----------------------------------------------------------------------------------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------------------------------------------------------------------------
  &*β* ~*d*~ = 0 and *ϕ* ~*b*~ = 0                                                                                                                                        $R_{0}^{b} > \left( {1 + \frac{\left( {\gamma_{a} + \mu} \right){\overset{-}{I}}_{a}R_{0}^{a}}{\left( {\gamma_{b} + \mu} \right)}} \right)R_{0}^{a} > R_{0}^{a}$
                                                                                                                                                                          
  ^\$^ *β* ~*d*~ \> 0 and *ϕ* ~*b*~ = 0                                                                                                                                   *R* ~0~ ^d^ \> *R* ~0~ ^a^
                                                                                                                                                                          
  *β* ~*d*~ = 0 and *ϕ* ~*b*~ \> 0                                                                                                                                        ^%^when *R* ~0~ ^b^ = *R* ~0~ ^a^, $\phi_{b} > \frac{1}{2}\frac{\gamma_{d} + \mu}{\gamma_{b} + \mu}$
  ^\#^when  *R* ~0~ ^b^ \< *R* ~0~ ^a^, $\phi_{b} > \frac{\gamma_{d} + \mu_{a}}{\beta_{a}{\overset{-}{I}}_{a}}\frac{1 - \Gamma}{1 + \Gamma}$                              
                                                                                                                                                                          
  *β* ~*d*~ \> 0 and *ϕ* ~*b*~ \> 0                                                                                                                                       ^%^when *R* ~0~ ^b^ = *R* ~0~ ^a^, $\frac{R_{0}^{d}}{R_{0}^{a}} + 2\phi_{b}\frac{\gamma_{b} + \mu}{\gamma_{d} + \mu} > 1$
  ^\#^when  *R* ~0~ ^b^ \< *R* ~0~ ^a^, $\frac{R_{0}^{d}}{R_{0}^{a}} + \phi_{b}\frac{\beta_{a}{\overset{-}{I}}_{a}}{\gamma_{d} + \mu}\frac{1 + \Gamma}{1 - \Gamma} > 1$   

We consider the successful invasion conditions by investigating the scenario where the host population in which strain *a* was already established is challenged by strain *b*. In principle, it is simply that the invasion reproductive number is larger than unit (i.e., *R* ~0~ ^b^(*Q* ~a~\*) \> 1). In [Table 4](#tab4){ref-type="table"} we explore the influence of the coinfection on the successful invasion by illustrating different special situations.

^&^ ${\overset{-}{I}}_{a}$ is the fraction of infected individuals at equilibrium *Q* ~a~\*.

^\$^Strain *a* will be replaced by the coinfection.

^\#^ $\Gamma \equiv {{({R_{0}^{b}/R_{0}^{a}})}{\lbrack{1 + {(({{\gamma_{a} + \mu)}/({\gamma_{b} + \mu)}})}{\overset{-}{I}}_{a}R_{0}^{a}}\rbrack}}^{- 1}$.

^%^The approximation under the condition $\gamma_{b} + \mu \gg \beta_{a}{\overset{-}{I}}_{a}$.

[^1]: Academic Editor: Laith Abu-Raddad
