The Impact of Coinfections and Their Simultaneous Transmission on Antigenic Diversity and Epidemic Cycling of Infectious Diseases

Abstract

Epidemic cycling in human infectious diseases is common; however, its underlying mechanisms have been poorly understood. Much effort has been made to search for external mechanisms. Multiple strains of an infectious agent were usually observed and coinfections were frequent; further, empirical evidence indicates the simultaneous transmission of coinfections. To explore intrinsic mechanisms for epidemic cycling, in this study we consider a multistrain Susceptible-Infected-Recovered-Susceptible epidemic model by including coinfections and simultaneous transmission. We show that coinfections and their simultaneous transmission widen the parameter range for coexistence and coinfections become popular when strains enhance each other and the immunity wanes quickly. However, the total prevalence is nearly independent of these characteristics and approximated by that of one-strain model. With sufficient simultaneous transmission and antigenic diversity, cyclical epidemics can be generated even when strains interfere with each other by reducing infectivity. This indicates that strain interactions within coinfections and cross-immunity during subsequent infection provide a possible intrinsic mechanism for epidemic cycling.

1. Introduction

Infectious diseases cause huge morbidity and mortality in human populations and remain as a dramatic threat to public health [1, 2]. 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 [3] to faecal-oral infections such as cholera and rotavirus and vector-borne diseases such as malaria and dengue fever [4, 5]. For example, seasonal flu develops as an epidemic during winter in temperate regions but it remains at very low levels during summer [6] while pertussis has 3-4-year cycles [7]. 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 [4, 6, 8]: 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”) [4, 5, 8–12]. However, these models cannot provide a reasonable explanation for, say, a quickly globally spreading pattern of influenza [13] and require unrealistically high reproduction rates for the observational seasonal patterns in influenza [14]. Moreover, Grassly et al. [15] 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 [16]. 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 [17], or models including an isolated class generate sustained oscillations albeit under unrealistic region of parameter space [18]. 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 [10, 19] or is impossible [20].

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 [21–26] and concurrent infections with multiple strains [27–31]. Polymorphic pathogens such as dengue virus [22], malaria parasite Plasmodium falciparum [21], influenza virus [23, 24], and Streptococcus pneumoniae [26] are responsible for enormous health burdens of human populations [26, 32, 33]. 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 [25]. In practice, this will also be supplied with the basis for the design of effective control strategies such as vaccination [34].

In general, once two different pathogen strains infect the same host individuals, they will interact with each other [32]. 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 [35]. The indirect interactions may be synergistic, neutral, or antagonistic in view of changes in transmissibility and virulence [2, 33, 35–37]. For example, epidemiological evidence suggests that primary infection with one strain of dengue virus can enhance transmission of subsequent infection with another strain [38–40]. Further, the strains present in coinfections can transmit simultaneously and separately [29, 41, 42]. 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. [29]. Cross-immunity between different subtypes of influenza has been observed [43–46]. 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 [47]. Cross-immunity was also found between different strains of other infectious agents [48, 49]. 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 [50]; further, coinfections may imply that the survivorship of strains is not necessarily maximized by maximizing the basic reproductive number [51].

Transmission consequences of coinfection and superinfection have received a great deal of interest over the last 30 years [37, 52–58]. Theoretical modelling and empirical analyses show that pathogen interactions could have profound implications for the pathogen communities [51, 59–61]. For example, it has been shown that cocirculation of different diseases and cross-immunity [59] and antibody-dependent enhancement in transmission [38, 39, 62] can induce cyclical and chaotic epidemics. Cross-immunity alone within a two-strain SIR or SIRS model cannot generate sustained oscillation [12, 63–65] unless other factors such as age-dependent survival [53, 57] or isolation [66] are included. Following the traditional assumption that individuals, while being infected with one strain, are temporarily not susceptible to the other [52, 53, 62, 67], 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 [12, 58, 61, 67, 68].

Recent empirical evidence has challenged the traditional assumption and indicated the concurrent infections during the infectious period of primary infection [29, 32, 41, 69, 70]. Due to technical difficulty and complexity involved in distinguishing multiple strains, multiple-strain infections are surely underestimated as argued by Balmer and Tanner [32]. The available empirical studies suggest that the occurrence of multiple strains was significantly more common than expected by chance. For example, Vaccarella et al. [70] 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 [42, 61]. 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 [71]. 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. [59, 72]). In view of evidence of strong cross-immunity among strains (e.g., [47]), 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. [38, 62, 67]). 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 [37, 42, 54]. 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

We consider a human population in which there are 𝑛 cocirculating strains and every individual is capable of being simultaneously infected by at most 𝑚(≤𝑛) strains. For convenience, we denote the coinfection of multiplicity 𝑙 as 𝑖₁𝑖₂⋯𝑖_𝑙, here 𝑖₁,𝑖₂,…,𝑖_𝑙∈𝐻≡{1,2,3,…,𝑛} and 𝑖₁<𝑖₂<⋯<𝑖_𝑙 representing the strains present in the coinfection, which is independent of order of such infections. Let 𝐼_{𝑙,𝑖1⋯𝑖𝑙} denote the fraction of individuals that carry the coinfection 𝑖₁𝑖₂⋯𝑖_𝑙. In this study we define the multiplicity of strains in a coinfection as the within-in host antigenic diversity and the number (𝑛) of cocirculating strains in the host population as the antigenic diversity of the population. We consider a SIRS epidemic model by assuming the coinfection 𝑖₁𝑖₂⋯𝑖_𝑙 can simultaneously transmit among individuals at rate 𝛽_{𝑙,𝑖1⋯𝑖𝑙}; the coinfection simultaneously progresses to the recovered at rate 𝛾_{𝑙,𝑖1⋯𝑖𝑙} (cf. [34, 73]); strain interactions within coinfection occur via the change in infectivity (cf. [34]). Infections are acquired through a mass action term with the force of infection for the coinfection 𝑖₁𝑖₂⋯𝑖_𝑙 given by Λ_{𝑙,𝑗1⋯𝑗𝑙}=𝛽_{𝑙,𝑗1⋯𝑗𝑙}𝐼_{𝑙,𝑗1⋯𝑗𝑙}+^𝑚_𝑣=𝑙+1𝜙_{𝑣,𝑙{𝑗1,…,𝑗𝑙}⊆{𝑘1,…,𝑘𝑣}}𝐼_{𝑣,𝑘1⋯𝑘𝑣}. (1) Here 𝜙_𝑣,𝑙 represents the factor that 𝑙 strains transmit simultaneously from a patient coinfected with 𝑣(>𝑙) strains, which reflects the strain interaction within coinfection. If strains interfere with each other, 𝜙_𝑣,𝑙<1; otherwise, 𝜙_𝑣,𝑙>1 implies enhancement in infectivity of strains. The factor 𝜙_𝑣,𝑙 here is distinct from that proposed in dengue virus models (e.g., [38, 62]) 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. [51, 61]). This assumption neglects the complicated history of infection and thus largely simplifies the system (cf. [51, 63, 67, 74]). 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 𝑑𝑆𝑑𝑡=−^𝑚_{𝑣=1{𝑘1,…,𝑘𝑣}}Λ_{𝑣,𝑘1⋯𝑘𝑣}𝑆+𝜇(1−𝑆)+𝜎1−𝑆−^𝑚_{𝑙=1𝑖1…𝑖𝑙∈𝐻}𝐼_{𝑙,𝑖1⋯𝑖𝑙},𝑑𝐼_{𝑙,𝑖1⋯𝑖𝑙}=𝑑𝑡_{𝑙𝑣=1{𝑗1,…,𝑗𝑣}∪{𝑗𝑣+1,…,𝑗𝑙}={𝑖1,…,𝑖𝑙}}Λ_{𝑣,𝑗1⋯𝑗𝑣}𝐼_{𝑙−𝑣,𝑗𝑣+1⋯𝑗𝑙}−^𝑚−𝑙_{𝑣=1{𝑗1,…,𝑗𝑣}∩{𝑖1⋯𝑖𝑙}=∅}Λ_{𝑣,𝑗1⋯𝑗𝑣}𝐼_{𝑙,𝑖1⋯𝑖𝑙}−𝛾_{𝑙,𝑖1⋯𝑖𝑙}+𝜇𝐼_{𝑙,𝑖1⋯𝑖𝑙}. (2) Here ∅ is used to stand for empty set. In the above equations, 𝑙=1,…,𝑚, and 𝑆(=𝐼₀) represents the fraction of the susceptible. For coinfection of multiplicity 𝑙, there are 𝑛_𝑙=𝑛!/𝑙!(𝑛−𝑙)! 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,𝑖=𝛽₁ and 𝛾_1,𝑖=𝛾₁,  𝑖=1,…,𝑛 (cf. [51, 67]). Further, strain interactions are identical and depend only on the multiplicity of coinfections: 𝛽_2,𝑖𝑗=𝛽₂,…,𝛽_{𝑚,𝑖1𝑖2⋯𝑖𝑚}=𝛽_𝑚,  𝜙_𝑗,𝑙=𝜙_𝑙 for 𝑗>𝑙 and 𝛾_2,𝑖𝑗=𝛾₂,…,𝛾_{𝑚,𝑖1𝑖2⋯𝑖𝑚}=𝛾_𝑚. At the equilibrium of such symmetric system the fraction of individuals that were coinfected with 𝑚 strains is the same among different combinations; for example, 𝐼_11,𝑖=𝐼₁,…,𝐼_{𝑚,𝑖1𝑖2⋯𝑖𝑚}=𝐼_𝑚,𝑖₁,…,𝑖_𝑚∈𝐻. The model system reduces to that of 1+𝑚 equations 𝑑𝑆𝑑𝑡=−^𝑚_𝑖=1𝑛𝑖Λ_𝑖𝑆+𝜇(1−𝑆)+𝜎1−𝑆−^𝑚_𝑙=1𝑛_𝑙𝐼_𝑙,𝑑𝐼_𝑙=𝑑𝑡^𝑙_𝑖=1𝑙𝑖Λ_𝑖𝐼_𝑙−𝑖−^𝑚−𝑙_𝑗=1𝑗Λ𝑛−𝑙_𝑗𝐼_𝑙−𝛾_𝑙+𝜇𝐼_𝑙,𝑙=1,…,𝑚, (3) with the force of infection, Λ_𝑙=𝛽_𝑙𝐼_𝑙+𝜙_𝑙^𝑚_𝑗=𝑙+1𝐼𝑛−𝑙𝑗−𝑙^𝑗,𝑙=1,…,𝑚. (4) The prevalence of coinfections of multiplicity 𝑙 is 𝑝_𝑙=𝑛_𝑙𝐼_𝑙 and the total prevalence is 𝑝=∑^𝑚_𝑙=1𝑝_𝑙. For further simplicity we assume that all the efficiencies of partial transmission from coinfections are the same (i.e., 𝜙₁=𝜙₂=⋯=𝜙_𝑚=𝜙) and the infectious periods of all infections are the same (i.e., 1/𝛾₁=1/𝛾₂=⋯=1/𝛾_𝑚=1/𝛾).

3. Results

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) and (4). 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

Consider a situation where a host population was already at equilibrium with 𝑛_𝑒 endemic strains, each of which has a basic reproductive number 𝑅₀=𝛽₁/(𝜇+𝛾). What is the condition for a new indistinguishable strain to invade into the population? In the case of 𝑛_𝑒=1, as shown in Appendix, it is that the efficiency of partial transmission from coinfections (𝜙) and the simultaneous transmission rate of double infections (𝛽₂) satisfy the inequality 𝜙≥𝜙_𝑠𝛽≡0.5−²2𝛽₁. (5) 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) was analytically obtained for the system of two strains, numerical simulations confirm it as a good approximation for multistrain systems (see Table 1). Results show that when the efficiency of partial transmission from coinfections is less than the critical value (𝜙_𝑠), 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) is satisfied, arbitrary number of strains can be maintained. Table 1 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., 𝜙_𝑠 is much smaller than the unit) (cf. [73]). The critical value (𝜙_𝑠) is 50% in the absence of simultaneous transmission of double infections; with simultaneous transmission (𝛽₂), the critical value (𝜙_𝑠) further reduces as predicted in (5) (Table 1). 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. [67]). Hence it is the strain interaction within coinfection (i.e., 𝛽₂ and 𝜙) that changes the traditional competitive exclusion [67, 73] and facilitates the increase in antigenic diversity.

Table 1 

The dependence of the successful invasion by an indistinguishable new strain on simultaneous transmission rates of coinfections of multiplicity 𝑚=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 𝑛_𝑒 endemic strains. After the challenge with a new identical strain, a new endemic will be established if condition (5) 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), 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.0, and infectious period = 2 days; hence 𝑅₀=2.0 for single infection.

3.2. Structure of Antigenic Diversity

Under the traditional circumstance where only one strain can transmit (i.e., 𝛽₁>0 but 𝛽_𝑙=0 for 𝑙>1), the antigenic diversity structure is shown in Figure 1. It shows that as the antigenic diversity (𝑙) of coinfections increases, the proportion of coinfections of the same 𝑙 strains (i.e., 𝐼_𝑙) decrease as 1/𝑛_𝑙; however, their prevalence (i.e., 𝑝_𝑙) remains roughly the same. The reduction in 𝐼_𝑙 with multiplicity 𝑙 is so remarkable that single infections are the most common infections (Figures 1(a) and 1(b)). When transmissibility is relatively low, triple and quadruple infections are prohibited (Figures 1(a) and 1(c)). However, with a high transmissibility, high diverse infections can safely survive (Figures 1(b) and 1(d)). 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 𝑝=^𝑚_𝑙=1𝑝_𝑙≈𝜎+𝜇1𝜎+𝜇+𝛾1−𝑅₀ (6) (see Figure 1(e)). This approximation can be simply confirmed from (3) and (4) under the situation where coinfections are rare. Ignoring coinfections, (3) and (4) reduce to that of one-strain SIR model with the prevalence of one strain replaced by total prevalence of 𝑛 strains (cf. [67]). The derivation of approximation (6) for a special situation of two strains is given in Appendix.

(a) 𝑅0=2 and 1/𝜎 = 10 yrs

(b) 𝑅0=20 and 1/𝜎 = 10 yrs

(c) 𝑅0=2 and 𝑛=50

(d) 𝑅0=20 and 𝑛=50

(e) 𝑛=50 and 1/𝜎 = 10 yrs

Figure 1 

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 𝛽𝑙=0,  𝑙≥2,  𝜙=1, and life span (1/𝜇) = 70 years. In panels (a) and (c) infectious period = 2 days and 𝑅0=2.0; in panels (b) and (d) infectious period = 20 days and 𝑅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 𝑝 is well approximated by (6). Note that in panel (d) the prevalence 𝑝4 exceeds 𝑝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.

The influences of simultaneous transmission of coinfections are shown in Figures 2 and 3 under two different immunity periods. Increasing 𝛽₂ increases the prevalence of coinfections but reduces single infections when strains interfere with each other (𝜙<1) (Figures 2(a) and 3(a)). However, if strains interact to enhance the infectivity within coinfections (𝜙>1), cyclical epidemics may emerge when immunity lasts longer (Figure 2(b)) or the prevalence of highly diverse coinfections increases substantially while single infections decrease greatly (Figure 3(b)). However, the impact of simultaneous transmission of more complicated coinfections is much weaker such that 𝛽₄ hardly have any effect (Figures 2(e), 2(f), 3(e), and 3(f)). Under the circumstance of 𝛽₂>𝛽₃>𝛽₄, 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)) and 𝜙<1 (Figures 2(a), 2(c), 2(e), and 3(a)). However, when strains facilitate each other to enhance their infectivity (𝜙>1), single infections cannot dominate (Figures 2(b), 2(d), and 2(f)). 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), 3(d), and 3(f)). This indicates that decreasing the duration of immunity can increase coinfection and stabilize the coexistence (cf. Figures 2(b) and 3(b); see also Figure 5(a)). When the coexistence is steady, the total prevalence (𝑝) slightly increases with these characteristics. Given all others being the same, increasing 𝛽₂ from 0 to (𝛽₁=)1 only increases 𝑝 less than 3% while increasing 𝛽₃ from 0 to 0.75 and increasing 𝛽₄ from 0 to 0.6 only increase 𝑝 by 0.9% and 0.3%, respectively.

(a) 𝜙=0.5, 𝛽3=0.0, 𝛽4=0.0

(b) 𝜙=2, 𝛽3=0.0, 𝛽4=0.0

(c) 𝜙=0.5, 𝛽2=0.75, 𝛽4=0.0

(d) 𝜙=2, 𝛽2=0.75, 𝛽4=0.0

(e) 𝜙=0.5, 𝛽2=0.75, 𝛽3=0.6

(f) 𝜙=2, 𝛽2=0.75, 𝛽3=0.6

Figure 2 

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 (𝛽𝑙) of coinfections of multiplicity are shown: panels (a) and (b) 𝑙=2; (c) and (d) 𝑙=3; and (e) and (f) 𝑙=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 𝑅0=20; life span (1/𝜇) = 70 years, 𝑛=10, and 𝑚=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 𝑀) and minimum (with subscript 𝑚) prevalence are shown.

(a) 𝜙=0.5,  𝛽3=0.0,  𝛽4=0.0

(b) 𝜙=2,  𝛽3=0.0,  𝛽4=0.0

(c) 𝜙=0.5,  𝛽2=0.75,  𝛽4=0.0

(d) 𝜙=2,  𝛽2=0.75,  𝛽4=0.0

(e) 𝜙=0.5,  𝛽2=0.75,  𝛽3=0.6

(f) 𝜙=2,  𝛽2=0.75,  𝛽3=0.6

Figure 3 

As in Figure 2 except immunity period = 3 years. Under this short immunity period, no cyclical epidemics appear within the ranges of parameters investigated.

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). When 𝜙 is strong, given all others being the same, coinfections such as double and triple coinfections become common and predominate (Figure 3) and may also induce cyclical epidemics (Figure 2(b)). Nevertheless, as 𝜙 increases from 0 to 2, coinfections predominate but the total prevalence 𝑝 only increases about 5%. Overall, the results cannot exceed 106% of (6) for all the steady scenarios examined in Figures 2 and 3. These observations are also supported by the results listed in Table 2. This indicates that (6) is a good approximation for the total prevalence of steady coexistence. Thus strain interactions (i.e., 𝛽₂ and 𝜙) determine the structure of antigenic diversity but cannot much alter the total prevalence.

Table 2 

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.0, and infectious period = 2 days; hence 𝑅₀=2.0 for single infection; ten indistinguishable strains (𝑛=10) are assumed to cocirculate while three models are considered in which the most complex coinfections are of multiplicity 𝑚=2,3,4 strains, respectively.

3.3. Cyclic Epidemics

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 shows the impact of simultaneous transmission of coinfections and antigenic diversity on the critical value of efficiency of partial transmission from coinfections (𝜙_𝑐) for emergence of cyclical epidemics. In the absence of simultaneous transmission of coinfections, 𝜙_𝑐>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 [38, 62, 75]. Two previous multistrain models [59, 63] 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 (𝛽₂) and high antigenic diversity, 𝜙_𝑐 reduces to below unit. For example, when 𝛽₂=𝛽₁/2=0.5,  𝜙_𝑐=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) and 4(c)). These indicate that 𝛽₂ has a predominant effect on the critical value 𝜙_𝑐. 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.

(a)

(b) 𝑚=𝑛 when 𝑛≤4; 𝑚=4 when 𝑛>4

(c) 𝑛=10, 𝑚=4

(d) 𝜙=𝜙𝑐

Figure 4 

Impact of simultaneous transmission of coinfections and antigenic diversity on the critical value of the efficiency of partial transmission from coinfections (𝜙𝑐) 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 𝑛=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 𝑅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 (𝑅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).

(a)

(b) 𝜙=𝜙𝑐

(c)

Figure 5 

Impact of the basic reproductive number and immunity period on the critical value of the efficiency of partial transmission from coinfections (𝜙𝑐) 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 𝜙=𝜙𝑐. In panel (c) the interepidemic periods in years are shown against 𝜙 and 𝑅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, 𝑛=10, and 𝑚=4. The basic reproductive number of single infection (𝑅0) varies by altering the infectious period (1/𝛾).

Figure 5 shows that cyclic epidemics are also controlled by other model parameters. For a given immunity period, Figure 5(a) shows that cyclic epidemics become impossible for a too small or too large reproductive number 𝑅₀. 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) and 5(c) and Figure 6). However, infectious period can only slightly increase the interepidemic period (data not shown). Figure 6 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) and the fraction of complicated coinfections can become popular as seen in steady coexistence shown in Figure 1(d).

(a) 𝜙=0.6

(b) 𝜙=1.0

(c) 𝜙=2.0

Figure 6 

Example time series for the proportion of infected individuals. Other model parameter values are 𝛽1=1.0 day−1, 𝑅0=15,  1/𝜇=70 years, 𝑛=10, and 𝑚=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) 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.

4. Discussion

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

In contrast to other coinfection models [51], 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 and 8; [40]) and simultaneous transmission enables a less transmissible strain to successfully invade (Table 4). 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; [67, 73]). This is comparable with multistrain models for dengue virus which include cross-reactive interaction during subsequent infection by other strains [76, 77]. Kooi et al. [76] 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. [77] 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. [51], 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, 2, and 3). 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 [25].

Figure 7 

Impact of strain interactions on coexistence and sustained oscillation in incidence of infections. The host population was initially endemic with strain 𝑎 and then challenged with strain 𝑏. Two strains are assumed to be equally transmissible with 𝛽𝑎=𝛽𝑏=0.75 day−1 and infectious period (1/𝛾𝑎=1/𝛾𝑏) = 2 days (i.e., 𝑅𝑎0=𝑅𝑏0=1.5). Coinfection can simultaneously transmit at a rate (𝛽𝑑)=0.375 day−1 and with an infectious period (1/𝛾𝑑)=2 days. Other parameters are the immunity period (1/𝜎)=3 years and life span (1/𝜇)=55 years.

(a) 𝑅𝑑0=0 and 𝜙≤0.5 or 𝑅𝑑0=𝑅𝑎0/2 and 𝜙≤0.25

(b) 𝑅𝑑0=0 and 𝜙=5

(c) 𝑅𝑑0=𝑅𝑎0/2 and 𝜙=5

Figure 8 

Impact of efficiency of single transmission (𝜙) and simultaneous transmissibility (𝑅𝑑0) from dual infection on the distribution of equilibria. The four equilibria on the (𝑅𝑎0,𝑅𝑏0) plane are separated by four curves: 𝑅𝑎0=1, 𝑅𝑏0=1, 𝑅𝑏0(𝑄∗𝑎)=1, and 𝑅𝑎0(𝑄∗𝑏)=1. In panel (a) where 𝑅𝑑0=0 and 𝜙≤0.5 or 𝑅𝑑0=𝑅𝑎0/2 and 𝜙≤0.25,  𝑅𝑏0>𝑅𝑎0 is required to maintain 𝑅𝑏0(𝑄∗𝑎)≥1 (and similarly requiring 𝑅𝑎0>𝑅𝑏0 to maintain 𝑅𝑎0(𝑄∗𝑏)≥1); hence no coexistence is possible. In panels (b) and (c) we consider the symmetric situation in which 𝜙𝑎=𝜙𝑏=𝜙=5 under different simultaneous transmissibility. The other parameters are infectious period (1/𝛾𝑎=1/𝛾𝑏=1/𝛾𝑑)=2 days, immunity period = 3 years, and life span = 55 years.

The observed levels of antigenic diversity may be large but surely are not infinite, say, influenza A [78]. So what limits the antigenic diversity? As the total prevalence is roughly constrained by (6) 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 [25, 67, 74]. Further, when infections experience oscillations, trough extinction (Figure 6 and Figure 12; [65]) 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 𝛽₂ 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

Though coinfection enhances coexistence of pathogen strains it alone cannot induce cyclical epidemics (Figure 5(a)). For example, Zhang et al. [36] 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 [79] and Martcheva and Pilyugin [57] 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 [2, 3] and is an important factor that leads to sustained oscillation [11, 15, 20]. Empirical analyses show that, for example, vaccination increases the degree of seasonal variation in the incidence of rotavirus [5]. However, theoretical investigations show that cross-immunity alone within a two-strain SIR model cannot support cyclical epidemics [63, 66, 80]. 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. [2, 38]). This is comparable with the assumption of transient strain-transcending immunity [74, 78] 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); [59]). Figure 11 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. [59, 63]). In the presence of simultaneous transmission, strain interactions within coinfections interplay with cross-immunity to generate cyclical epidemics (Figures 4 and 5). The prevalent view of intrinsic mechanisms for cyclical epidemics in the situation of strong cross-immunity is the enhancement in sequential infections between strains [38, 59, 62, 63, 75]. 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 and 5).

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 𝑅₀ gives rise to cyclical epidemics of short periods (Figure 5). This prediction is consistent with the results of Ferguson et al. [38] 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 [12, 81] and 𝑅₀ about 2 [82] was bound to annual cycles [24] while pertussis of a longer immunity period of about 30 years [83] leads to 3-4-year cycles [7]. Moreover, the interepidemic period initially increases with efficiencies of partial transmission from coinfections but quickly saturates (Figure 5(c) and Figure 6). It first increases with simultaneous transmission rate of double infections (𝛽₂) and then decreases (Figure 12) while it is insensitive to simultaneous transmission of higher diverse coinfections when epidemic cycle first emerges (Figure 4(d)). Unfortunately, we currently do not have empirical data to test these relationships.

4.3. Comparison to Other Models of Epidemic Cycling

Oscillations in incidence were predicted in many previous epidemic models with different underlying mechanisms, for example, in one-strain SIR models [4, 8–11], two-strain SIR models [68], and two-strain SIRS models [12] where the seasonal forcing was thought to be the mechanism. Castillo-Chavez et al. [53] established that age-structure and cross-immunity generate sustained oscillations in a two-strain SIR model while Nuño et al. [66] showed that host isolation and cross-immunity lead to periodic epidemic outbreaks. Gupta et al. [59] 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 𝛽_𝑙=0 with 𝑙>1 within our model. Gupta et al. [59] 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): cyclical and chaotic epidemics are prohibited when simultaneous transmission of and partial transmission from coinfections are absent. Ferguson et al. [38] and others [39, 62, 75] 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 [84] and Bianco et al. [85] 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. [38] 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. [66]: 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. [61] 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 [53, 57] 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

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 [1, 73, 86, 87] 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 [15, 88]. 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 [4, 8, 19]. 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 [13] and especially for the cyclical dynamics of syphilis [15]. It has been suggested that the mechanisms for the common oscillation of infectious diseases must be multiple [6, 10, 19]. The true mechanism might lie in the interplay of the intrinsic (as suggested by this investigation and Grassly et al. [15]) and external periodic forcing [4, 8].

For simplicity, we assume complete cross-immunity during subsequent infection. In reality, cross-immunity is more likely to be partial [43–45, 89] 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.

Appendix

A. Two-Strain SIRS Model

In this appendix we consider a special scenario with 𝑛=𝑚=2 but allow two strains to differ in epidemiological properties (i.e., so-called asymmetric two-strain model). From the model systems (1) and (2), the two-strain SIRS model is given by 𝑑𝑆𝛽𝑑𝑡=𝜇−_𝑎𝐼_𝑎+𝜙_𝑎𝐼_𝑑+𝛽_𝑏𝐼_𝑏+𝜙_𝑏𝐼_𝑑+𝛽_𝑑𝐼_𝑑+𝜇𝑆+𝜎1−𝑆−𝐼_𝑎−𝐼_𝑏−𝐼_𝑑,𝑑𝐼_𝑎𝑑𝑡=𝛽_𝑎𝐼_𝑎+𝜙_𝑎𝐼_𝑑𝑆−𝛽_𝑏𝐼_𝑏+𝜙_𝑏𝐼_𝑑𝐼_𝑎−𝛾_𝑎+𝜇𝐼_𝑎,𝑑𝐼_𝑏𝑑𝑡=𝛽_𝑏𝐼_𝑏+𝜙_𝑏𝐼_𝑑𝑆−𝛽_𝑎𝐼_𝑎+𝜙_𝑎𝐼_𝑑𝐼_𝑏−𝛾_𝑏+𝜇𝐼_𝑏,𝑑𝐼_𝑑𝑑𝑡=𝛽_𝑏𝐼_𝑏+𝜙_𝑏𝐼_𝑑𝐼_𝑎+𝛽_𝑎𝐼_𝑎+𝜙_𝑎𝐼_𝑑𝐼_𝑏+𝛽_𝑑𝐼_𝑑𝑆−𝛾_𝑑+𝜇𝐼_𝑑, (A.1) where 𝐼_𝑎,  𝐼_𝑏, and 𝐼_𝑑 are the proportions in infections with strain 𝑎 alone, strain 𝑏 alone, and both. 𝛽_𝑎,  𝛽_𝑏, and 𝛽_𝑑 are transmission rate of strain 𝑎, strain 𝑏, and both, respectively, and 𝛾_𝑎,  𝛾_𝑏, and 𝛾_𝑑 are their corresponding recovery rates.

There are four equilibria (see Table 3); 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.

Table 3 

Four equilibria and their stability conditions of model system (A.1).

Table 4 

Impact of coinfection on invasion conditions.

A.1. Invasion Conditions and Stability Analysis of Equilibria

The threshold condition for pathogen strain 𝑎 to successfully invade a naive host population was that its basic reproductive number (𝑅^𝑎₀) must be larger than one where 𝑅^𝑎₀ 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 𝑏 to invade and persist in a host population in which strain 𝑎 was established at equilibrium (𝑄^∗_𝑎) is its invasion (or effective) reproductive number being larger than unit; that is, 𝑅^𝑏₀(𝑄^∗_𝑎)>1 [57]. Applying the method of van den Driessche and Watmough [90], we obtain 𝑅^𝑏₀𝑄^∗_𝑎=12Φ+Φ²+4Δ, (A.2) where 𝛽Φ≡^𝑏𝑅^𝑎₀1𝛾_𝑏+𝜇+𝛽_𝑎𝐼_𝑎+𝛽_𝑏𝜙_𝑏𝐼_𝑎+𝛽_𝑑/𝑅^𝑎₀𝛾_𝑑,+𝜇𝛽Δ≡^𝑏𝑅^𝑎₀𝜙_𝑏𝛽_𝑎𝐼_𝑎−𝛽_𝑑/𝑅^𝑎₀𝛾_𝑑+𝜇𝛾_𝑏+𝜇+𝛽_𝑎𝐼_𝑎. (A.3)(A.4) In expressions (A.3)-(A.4), 𝐼_𝑎 is the fraction of infective individuals with strain 𝑎 at equilibrium 𝑄^∗_𝑎;  𝑅^𝑑₀ represents the basic reproductive number of double infection. It is worth mentioning that 𝑅^𝑏₀(𝑄^∗_𝑎)=1 when Φ+Δ=1; thus, when Φ+Δ>1,  𝑅^𝑏₀(𝑄^∗_𝑎)>1. To demonstrate how the interaction between strains determines the successful invasion and coexistence, we rearrange the inequality Φ+Δ>1 into 𝛽_𝑑>𝛾_𝑑+𝜇−𝜙_𝑏𝛽_𝑎𝐼_𝑎1+Γ𝑅1−Γ^𝑎0 (A.5) with Γ≡(𝑅^𝑏₀/𝑅^𝑎₀)(1+(𝐼_𝑎𝛽_𝑎/(𝛾_𝑏+𝜇)))⁻¹. Under the situation where 𝑅^𝑎₀=𝑅^𝑏₀ and 𝛾_𝑏+𝜇≫𝛽_𝑎𝐼_𝑎 which equivalently requires that the infectious period is much shorter than the duration of immunity and life span, we have Γ≈1−𝐼_𝑎𝛽_𝑎/(𝛾_𝑏+𝜇), and the above inequality is approximated as follows: 𝛽_𝑑>𝛾_𝑑+𝜇−2𝜙_𝑏𝛾_𝑏+𝜇𝑅^𝑎₀. (A.6)

By the symmetry between strain 𝑎 and strain 𝑏, 𝑅^𝑎₀(𝑄^∗_𝑏) can be obtained by swapping the indexes 𝑎 and 𝑏 in the above equations. When 𝑅^𝑎₀(𝑄^∗_𝑏)>1 and 𝑅^𝑏₀(𝑄^∗_𝑎)>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) and the other for distribution of four equilibria and the effect of simultaneous transmission (Figure 8).

The stability conditions for equilibria can be obtained by analysing the eigenvalues of Jacobian matrix of the system (A.1) at equilibria. Consider the stability condition of equilibrium 𝑄^∗_𝑎 of single infection with strain 𝑎; its characteristic equation can be reduced to ||=𝐌−𝜆𝐈−𝛽_𝑎𝐼_𝑎−𝜇−𝜎−𝜆−𝛽_𝑎𝛽𝑆−𝜎_𝑎𝐼_𝑎×𝛽−𝜆_𝑏𝑆−𝛽_𝑎𝐼_𝑎−𝛾_𝑏−𝜇−𝜆𝜙_𝑏𝛽_𝑏𝑆𝛽_𝑎+𝛽_𝑏𝐼_𝑎𝜙_𝑏𝛽_𝑏𝐼_𝑎+𝛽_𝑑𝑆−𝛾_𝑑≡|𝐂|×|𝐃|−𝜇−𝜆=0. (A.7) Here 𝑆 and 𝐼_𝑎 are given in Table 3. Matrix C always has two negative eigenvalues while matrix D has negative eigenvalues only if the following inequalities hold: 𝛽_𝑏1𝑅^𝑏₀−1𝑅^𝑎₀+𝛽_𝑑1𝑅^𝑑₀−1𝑅^𝑎₀>−𝛽_𝑎−𝛽_𝑏𝜙_𝑏𝐼_𝑎,𝛽1>^𝑏𝑆𝛾_𝑏+𝜇+𝛽_𝑎𝐼_𝑎+𝛽_𝑏𝜙_𝑏𝐼_𝑎+𝛽_𝑑𝑆𝛾_𝑑+𝛽+𝜇_𝑎𝜙_𝑏𝐼_𝑎−𝛽_𝑑𝑆𝛾_𝑏+𝜇+𝛽_𝑎𝐼_𝑎𝛾_𝑑𝛽+𝜇𝑅^𝑏𝑎₀≡Φ+Δ. (A.8)(A.9) Inequality (A.9) is equivalent to 𝑅^𝑏₀(𝑄^∗_𝑎)<1 and it can also be expressed as follows: 𝑅^𝑎₀>𝑅^𝑏₀𝑅1−^𝑑0𝑅^𝑎₀𝛾_𝑏+𝜇𝛾_𝑏+𝜇+𝛽_𝑎𝐼_𝑎+𝑅^𝑑₀+𝑅¹₀𝛽_𝑏𝜙_𝑏𝐼_𝑎𝛾_𝑑𝛾+𝜇1+^𝑎+𝜇𝛾_𝑏+𝜇+𝛽_𝑎𝐼_𝑎. (A.10) Combining this with inequality (A.8), we have 𝑅^𝑎₀>1,  𝑅^𝑎₀>𝑅^𝑏₀, and 𝑅^𝑎₀>𝑅^𝑑₀. These are stability conditions listed in Table 3. Similarly we can obtain the stability conditions of equilibria 𝑄^∗_𝑏 and 𝑄^∗₀. Due to the complexity, the stability condition for coexistence equilibrium 𝑄^∗_𝑎𝑏 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 𝑅^𝑎₀<1 and 𝑅^𝑏₀<1, conditions that 𝑅^𝑑₀>1 plus 𝜙_𝑎>0 and 𝜙_𝑏>0 can guarantee the coexistence of both strains (cf. [90]).

Some special situations are listed in Table 4 to illustrate the effect of coinfections on the successful invasion. Without simultaneous transmission of coinfections, strain 𝑏 can invade only when it has a higher basic reproductive number 𝑅^𝑏₀ (see also Figure 8). However, strain interactions within coinfections change this outcome. If 𝑅^𝑏₀=𝑅^𝑎₀, strain 𝑏 can invade if the efficiency of transmitting strain 𝑏 from coinfection (𝜙_𝑏) is larger than (1/2)(1/𝛾_𝑏)/(1/𝛾_𝑑), suggesting that if coinfections have doubled infectious periods, strain 𝑏 can invade if 𝜙_𝑏>25%. On the other extreme scenario where 𝜙_𝑏 is nil, strain 𝑏 can invade only when 𝑅^𝑑₀>𝑅^𝑎₀ and the previous single infection will be replaced by coinfections. In particular, it is possible for a less transmissible strain 𝑏 (i.e., 𝑅^𝑏₀<𝑅^𝑎₀) to invade if both 𝛽_𝑑 and 𝜙_𝑏 are strong enough.

Next we consider a symmetric situation where two strains share the same properties: 𝛽≡𝛽_𝑎=𝛽_𝑏,  𝛾≡𝛾_𝑎=𝛾_𝑏=𝛾_𝑑 and 𝜙≡𝜙_𝑎=𝜙_𝑏 so that from (A.1) we have at equilibrium 𝐼₁=𝐼₂ and 𝐼₁+𝐼₂+𝐼_𝑑=𝜎+𝜇𝜎+𝜇+𝛾1−𝑆. (A.11) When 𝜙 is not very strong and 𝛽>𝛽_𝑑 so that only the steady coexistence is possible, numerical results show that 𝑆≫𝐼₁=𝐼₂≫𝐼_𝑑 (Tables 1 and 2); from the second and third equations of (A.1), we obtain the approximation 𝑆≈𝜇+𝛾𝛽=1𝑅^𝑎₀=1𝑅^𝑏₀. (A.12) Putting (A.12) into (A.11) gives rise to approximation (6) in the text for the situation of two symmetric strains.

A.2. Analysis of Hopf Bifurcation

Because of the complexity in model system (A.1), 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 (𝑄^∗_𝑎𝑏) is as follows: ||=𝑚𝐌−𝜆𝐈¹¹−𝜆𝑚₂₁𝑚₃₁𝑚₄₁𝑚₁₂𝑚₂₂−𝜆𝑚₃₂𝑚₄₂𝑚₁₃𝑚₂₃𝑚₃₃−𝜆𝑚₄₃𝑚₁₄𝑚₂₄𝑚₃₄𝑚₄₄−𝜆=0, (A.13) where elements 𝑚_𝑖𝑗 are functions of the interior equilibrium (𝑆,𝐼_𝑎,𝐼_𝑏,𝐼_𝑑,𝑅) and model parameters. The characteristic equation can be written as follows: 𝜆⁴+𝑎₃𝜆³+𝑎₂𝜆²+𝑎₁𝜆+𝑎₀=0. (A.14) The coefficients 𝑎_𝑘,  𝑘=0,1,…,3, are the functions of matrix elements 𝑚_𝑖𝑗,  𝑖,𝑗=1,…,4.

The quartic equation (A.14) 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 (𝜙_𝑎,𝜙_𝑏) of single transmission from the coinfection. In general, the four eigenvalues can be expressed as follows: 𝜆_𝑗𝜙_𝑎,𝜙_𝑏=𝜂_𝑗𝜙_𝑎,𝜙_𝑏+𝑖𝜔_𝑗𝜙_𝑎,𝜙_𝑏,𝑗=1,…,4. (A.15) When the real parts 𝜂_𝑗 are negative, the steady coexistence equilibrium is stable. However, when values of (𝜙_𝑎,𝜙_𝑏) increase to some critical values (𝜙_𝑎0,𝜙_𝑏0), real part 𝜂_𝑗 of some eigenvalues will increase to zero and the equilibrium loses its stability and a Hopf bifurcation occurs [91], from which oscillation in incidence will be induced. The imaginary part of the dominant eigenvalue, 𝜔(𝜙_𝑎,𝜙_𝑏), is the frequency of the oscillation when the real part 𝜂(𝜙_𝑎,𝜙_𝑏) passes through zero to become positive and the interepidemic period can be calculated as follows: 𝑇=2𝜋𝜔𝜙_𝑎,𝜙_𝑏. (A.16) One numerical example was shown in Figure 9. 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). Similarly, this can also be tested against other model parameters.

(a)

(b)

(c)

Figure 9 

Numerical integration of (A.1) with the efficiency of single transmission from dual infection 𝜙(=𝜙𝑎=𝜙𝑏)=1.75, 1.76, and 1.77. The real parts of two identical eigenvalues of characteristic equation (A.6) 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, 𝑇=410 days from (A.16), 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 𝛽𝑎=𝛽𝑏=0.75 day−1, 𝛽𝑑=0.375 day−1, infectious period (1/𝛾𝑎=1/𝛾𝑏=1/𝛾𝑑)=2 days, immunity period = 3 years, and life span = 55 years.

(a)

(b) 𝜙𝑎=𝜙𝑏=1.5

(c) 𝜙𝑎=𝜙𝑏=3.9

(d) 𝜙𝑎=𝜙𝑏=5.0

Figure 10 

(a) Bifurcation diagram plotting the local maxima of 𝐼𝑎 (the fraction of the infective individuals with strain 𝑎) against efficiency (𝜙=𝜙𝑎=𝜙𝑏) of single transmission from coinfection. Numerical integration of the equations is illustrated at different values of 𝜙𝑎=𝜙𝑏: (b) 1.5 (period-2), (c) 3.9 (period-4), and (d) 5.0 (chaotic) for 𝐼𝑎. In panel (b) the fraction of infective individuals with strain 𝑏 and that of coinfections are also included. The other parameters are immunity period = 5 years, 𝛽𝑎=0.75 day−1, 𝛽𝑏=0.74 day−1, 𝛽𝑑=0.7 day−1, and infectious period (1/𝛾𝑎=1/𝛾𝑏=1/𝛾𝑑)=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 (𝜙).

Figure 11 

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<𝜙𝑎=𝜙𝑏<9.2) and cyclical epidemic can emerge when 𝜙𝑎=𝜙𝑏>9.2 (data not shown). Other model parameters are 𝛽𝑎=0.75,  𝛽𝑏=0.74,  𝛽𝑑=0.70, infectious period (1/𝛾𝑎=1/𝛾𝑏=1/𝛾𝑑)=2 days, and life span (1/𝜇)=55 years.

(a) 𝑅𝑎0=𝑅𝑏0=3, immunity period = 3 years

(b) 𝑅𝑎0=𝑅𝑏0=3, immunity period = 3 years

(c) 𝑅𝑎0=𝑅𝑏0=17, immunity period = 50 years

(d) 𝑅𝑎0=𝑅𝑏0=17, immunity period = 50 years

Figure 12 

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) 𝑅𝑎0=𝑅𝑏0=3.0 and immunity period = 3 years, (c) and (d) 𝑅𝑎0=𝑅𝑏0=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/𝛾𝑎=1/𝛾𝑏=1/𝛾𝑑) 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.

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), the peak fraction of single infections with dominant strain 𝑎 is 0.24% while it is 0.04% at equilibrium if it alone infects the host population (from Table 1 and see also Figure 10(a)). Therefore we should be cautious when explaining the meanings and estimating values of the basic reproductive number (𝑅₀) 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), 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). Given low reproductive numbers for two different strains but relatively high transmissibility for double infection as in Figure 11, the critical efficiency of single transmission from double infection for cyclical epidemics is 𝜙_𝑎0=𝜙_𝑏0=4.0 for immunity that lasts for one month or 𝜙_𝑎0=𝜙_𝑏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) and 10(c)). While it first increases with the simultaneous transmissibility (𝑅^𝑑₀), after reaching peak at a critical value, it decreases (Figure 12).

Conflict of Interests

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

Acknowledgments

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.