PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
 
Epidemics. 2013 December; 5(4): 187–196.
PMCID: PMC3863957

Does homologous reinfection drive multiple-wave influenza outbreaks? Accounting for immunodynamics in epidemiological models[star]

Abstract

Epidemiological models of influenza transmission usually assume that recovered individuals instantly develop a fully protective immunity against the infecting strain. However, recent studies have highlighted host heterogeneity in the development of this immune response, characterized by delay and even absence of protection, that could lead to homologous reinfection (HR). Here, we investigate how these immunological mechanisms at the individual level shape the epidemiological dynamics at the population level. In particular, because HR was observed during the successive waves of past pandemics, we assess its role in driving multiple-wave influenza outbreaks. We develop a novel mechanistic model accounting for host heterogeneity in the immune response. Immunological parameters are inferred by fitting our dynamical model to a two-wave influenza epidemic that occurred on the remote island of Tristan da Cunha (TdC) in 1971, and during which HR occurred in 92 of 284 islanders. We then explore the dynamics predicted by our model for various population settings. We find that our model can explain HR over both short (e.g. week) and long (e.g. month) time-scales, as reported during past pandemics. In particular, our results reveal that the HR wave on TdC was a natural consequence of the exceptional contact configuration and high susceptibility of this small and isolated community. By contrast, in larger, less mixed and partially protected populations, HR alone cannot generate multiple-wave outbreaks. However, in the latter case, we find that a significant proportion of infected hosts would remain unprotected at the end of the pandemic season and should therefore benefit from vaccination. Crucially, we show that failing to account for these unprotected individuals can lead to large underestimation of the magnitude of the first post-pandemic season. These results are relevant in the context of the 2009 A/H1N1 influenza post-pandemic era.

Introduction

Mathematical models of infectious diseases often rely on a compartmental description in order to reduce the population diversity to a few key characteristics which are relevant to the infection under consideration. An extensively used model for influenza infection is of susceptible-exposed-infectious-removed (SEIR) form: after exposure to the virus, susceptible hosts (S) pass through an exposed state (E) of latent infection, become infectious (I) and are finally removed (R) from the infectious pool as they simultaneously recover (or die) and acquire permanent protection against the infecting strain. The SEIR model was particularly successful during the 2009 pandemic in estimating the key transmission parameters of the novel H1N1 virus (nH1N1) (Fraser et al., 2009) and assessing the effectiveness of alternative vaccination strategies (Baguelin et al., 2010).

Nevertheless, proper consideration of the primary immune response, which occurs on the first exposure to a novel influenza virus, motivates a more accurate description of the different stages from recovery to development of long-term protective immunity. Indeed, the primary immune response to influenza in humans operates on two different time scales. Usually, the viral load is cleared by the innate and cellular immune responses within a few days following infection (Woodland, 2003), thus leading to recovery of infected hosts. By contrast, the humoral (antibody-mediated) immune response, which provides long-term protection against the infecting strain as well as closely related strains (Fairlie-Clarke et al., 2008), takes several weeks to become efficient (Cox et al., 2004; Miller et al., 2010; Baguelin et al., 2011). Finally, at the population level, there is host heterogeneity in the development of this long-term protective immunity as some individuals show high antibody titres shortly after recovery whereas some other fail to reach a protective level (Cox et al., 2004; Miller et al., 2010; Chen et al., 2010; Hung et al., 2010; Chan et al., 2011).

In a recent study, Camacho et al. (2011) showed that a precise account of these host heterogenities was necessary to explain the reinfection episodes reported during the “natural experiment” of Tristan da Cunha (TdC), a remote island that underwent a two-wave A/H3N2 influenza epidemic in 1971 (Mantle and Tyrrell, 1973). More precisely, in the next few days that followed its introduction, the virus spread rapidly throughout the whole island population and after three weeks of propagation, 273 (96%) of 284 islanders had been infected. However, while the epidemic was nearing its end, several recovered islanders developed a second illness, thus initiating the second epidemic wave during which at least 92 (32%) islanders were reinfected (see section “Data” for more details). The main finding of Camacho et al. (2011) is that, among six biologically realistic reinfection mechanisms, only two could be retained: some hosts with either a delayed or deficient humoral immune response to the primary influenza infection were reinfected following rapid re-exposure to the same strain. This historical event illustrates that host heterogeneity at the individual level can not only lead to HR but also shape the epidemiological dynamics by triggering a second epidemic wave.

Historically, multiple-wave outbreaks and rapid reinfections have commonly been observed during influenza pandemics. The most striking example remains the “Spanish” influenza pandemic of 1918–1919 that occurred in three waves (Taubenberger and Morens, 2006) and during which several reinfection episodes were reported, sometimes in proportions similar to that of the 1971 TdC epidemic (Medical Department of the Local Government Board, 1919; Ministry Of Health, 1920; Barry et al., 2008). However, the three epidemic waves in 1918–1919 were spread out over 9 months (Taubenberger and Morens, 2006) whereas the two-wave epidemic on TdC lasted only 59 days (Mantle and Tyrrell, 1973). Accordingly, the time-scale between successive infections in the same individual was of the order of months during the pandemic whereas it was of the order of a few weeks for the TdC islanders, thus questioning their common underlying biological mechanisms. More recently, many populations experienced a spring and a fall waves during the 2009 pandemic and several cases of HR were virologically confirmed (Perez et al., 2010; Kim et al., 2010). Most of these HR episodes occurred within 2–3 weeks following recovery, a time-scale similar to that observed among the TdC islanders. However, both infection and HR occurred over the same epidemic wave in 2009 whereas they were separated across both waves in TdC, thus questioning the role of HR in driving multiple-wave outbreaks.

Overall, these observations call for clarification of the significance of HR and its role in driving multiple-wave outbreaks during pandemics. In particular, to what extent a better consideration of the immunological dynamics may be important in epidemiological models of influenza pandemics? In order to investigate these issues, we propose to explore and characterize the interplay between the immunological and epidemiological dynamics of a novel influenza virus. We start by defining an extended SEIR model accounting for the primary immune response to influenza and its inherent host heterogeneity. Using a maximum-likelihood (ML) approach, we confront our mechanistic model with the time-series of the daily incidence counts of the 1971 TdC epidemic and obtain ML estimates for the key immunological parameters. This analysis also reveals the exceptional setting of the TdC population and lead us to explore the impact of HR on the epidemiological dynamics for various population settings. We conclude with a discussion on the role of HR in the current post-pandemic era.

Materials and methods

The primary immune response to influenza infection in humans

A multi-pronged innate (McGill et al., 2009) and adaptive (Brown et al., 2004) immune response has been described for clearing influenza infection. The innate response is the first to be activated and plays a key role through its ability to control early viral replication and to promote and regulate the virus-specific adaptive immune response (McGill et al., 2009). The adaptive response itself may be broken into two critical sub-components: (i) the cellular immune response by which antigen-specific cytotoxic T lymphocytes (CTLs) eliminate infected cells and thus prevent viral release; and (ii) the humoral immune response by which serum and mucosal antibodies efficiently neutralize the virus (as explained in Text S1 the separation between serum and mucosal antibodies is not necessary for our study). Antibodies can remain detectable for years after infection and prevent reinfection by the same strain as well as by sufficiently cross-reactive variants (Fairlie-Clarke et al., 2008). Genetic variation in any of these immune components might determine whether or how rapidly an individual develops protective immunity following influenza infection.

As schematized in Fig. 1A, it is important to note that, during a primary influenza infection, the innate and cellular responses (blue curve) play the key role in viral clearance whereas neutralizing antibodies (green curve) are generated later and do not play a significant role unless the viral load is high and sustained (Woodland, 2003). The primary CTL response is detectable in blood after 6–14 days whereas the neutralizing antibody response peaks at 4–6 weeks (Cox et al., 2004). Critically, the CTL response is down-regulated after viral clearance (Woodland, 2003), disappears by day 21 post-infection (Cox et al., 2004) and is followed by a state of immunological “memory” with antigen-specific T cells. The memory cells cannot prevent HR as well as specific antibodies could, but they can reduce the severity of the disease (Woodland, 2003). Finally, it has been reported that a serum or mucosal antibody response cannot be detected in approximately 10 to 20% of subjects after natural influenza infection (Cox et al., 2004; Tamura and Kurata, 2004; Miller et al., 2010; Chen et al., 2010; Hung et al., 2010; Chan et al., 2011).

Fig. 1
Mechanistic modelling of the primary immune response to influenza. (A) Schematized dynamics of the viral load as well as the innate and adaptive immune responses, as described in section “The primary immune response to influenza infection in humans”. ...

Mechanistic modelling

Fig. 1B shows the SEICWH model which extends the classical SEIR model to account for the dynamics and host heterogeneity of the primary immune response to influenza in humans. Following recovery, hosts remain temporarily protected against HR thanks to the cellular response. Accordingly, they enter the C stage (cellular protection). Then, following down-regulation of the CTL response, the humoral response has a probability α to reach a level sufficient to protect against HR. In this case, recovered hosts enter the H stage (humoral protection) but otherwise they remain unprotected and re-enter the susceptible pool (S). Finally, in order to account for potential delay between completion of CTL contraction and full development of the neutralizing antibody response, recovered hosts pass through a time window of susceptibility (W) before entering the H stage. Crucially, while in the W stage, individuals can be reinfected following re-exposure to the same strain

In order to account for host heterogeneity in the development of the immune response, we use a stochastic framework to simulate the durations of the successive immunological stages. Defining τE, τI, τC and τW as the times spent by an infected host in the indexed immunological stages, we assume that these four random variables follow independent Erlang distributions with shapes kE, kI, kC and kW and means ϵ−1, ν−1, γ−1 and ω−1, respectively. Erlang distribution with mean m and shape k is modelled by k consecutive sub-stages, each being exponentially distributed with mean m/k. As illustrated in Fig. S1 the flexibility of the Erlang distribution ranges from the exponential (k = 1) to Gaussian-like (k [dbl greater-than sign] 1) distributions. In particular, whenever k > 1 the memory-less property of the exponential distribution is lost, thus providing more realistic distribution for biological processes with delays such as recovery or contraction of the CTL response (Wearing et al., 2005). In Fig. 1B, as well as in the rest of the paper, we use kE = kI = 2, kW = 1 and kC = 5. Justification of these values is provided in section “Parameter inference via maximum likelihood”.

Finally, regarding disease transmission, we make the standard assumptions of a well mixed, isolated and constant size (= Ω) population, as well as a constant contact rate (= β) among individuals. These simplifications permit us (i) to focus on the direct impacts of the immunological mechanisms on the epidemic dynamics; and (ii) to rapidly assess these impacts for various population settings, i.e. contact rates. We refer to the last section of this paper for a discussion on the inclusion of further refinements in the transmission mechanisms.

Data

The data counts are clinical records based on symptom observation and were drawn from the notes taken during the regular work of the local practice who visited all but three houses during the course of the epidemic in TdC (Mantle and Tyrrell, 1973). In addition, blood sample of 11 individuals provided serological confirmation of the circulation of A/H3N2 on the island, a subtype to which the TdC population had never been exposed before 1971. We refer to the paper of Mantle and Tyrrell (1973) for a detailed description of the 1971 data set as well as to the paper of Camacho et al. (2011) for a summary of the influenza experiences in TdC before 1971.

We note that the data are not available at the individual level. However, since the data set consists of 312 cases for 284 islanders some individuals must appear twice in the data counts. More precisely, in their original article, Mantle and Tyrrell (1973) states that among the 284 islanders 273 were infected at least once whereas 92 were reinfected. Unfortunately, only 312 of the 365 total cases were recorded with daily precision and were included in their data set (Mantle and Tyrrell, 1973) which is reproduced in Fig. 2 (black dots). As such, we can only conclude that at least 49 HR cases appear in the data.

Fig. 2
Detailed analysis of the 105 realizations of the stochastic SEICWH model for the 1971 TdC epidemic using Gillespie's algorithm (Gillespie, 1977). Upper panel: original incidence data (black dots) and expected incidence (red line) conditioned on non-extinction ...

Simulation

Our aim is (i) to fit the SEICWH model to the 1971 TdC epidemic in order to estimate the immunological parameters (ϵ, ν, γ, ω, α) and (ii) to explore the model dynamics for various population settings. For this purpose, we used both stochastic and deterministic simulations of the SEICWH model, as we now explain.

In a previous study, Camacho et al. (2011) showed that, given the small population of TdC, demographic stochasticity should be taken into account when fitting a mechanistic model to the 1971 TdC epidemic. This is because the risk of epidemic fade-out during the trough between waves depends on the model parameters and must therefore be accounted for when maximizing the likelihood over the parameter space. Accordingly, we exclusively resorted to stochastic simulations to fit the 1971 TdC epidemic. In the stochastic SEICWH model, the number of individuals in each immunological (sub-)stage is a discrete random variable and a possible state of the population at time t is defined by the random vector

equation M1

where C1:5 [equivalent] C1, …, C5 (and similarly for I1:2 and E1:2). The state of the population at time t is therefore a realization

equation M2

of X(t). The time course of x(t) is led by the possible transitions described in Table 1 and was simulated using Gillespie's exact algorithm (Gillespie, 1977).

Table 1
Transitions between classes in the stochastic SEICWH model. The notation A → B means that when the event occurs one individual is transferred from compartment A to compartment B.

Our second aim was to explore the dynamics of the SEICWH model over a wide range of parameter values. In this context, stochastic simulations become computationally intensive and one is tempted to resort on deterministic simulations for the sake of efficiency. However, this approximation is acceptable only as one controls that the stochastic effects should remain negligible. Accordingly, we assessed that the inter-wave extinction probability remains negligible (p < 10−3) in the parameter range explored (see Text S3), thus justifying the use of a deterministic approximation. In the deterministic SEICWH model, the state of the population x(t) becomes a continuous variable governed by a set of ordinary differential equations (given in Text S3) that can be obtained by the large population limit (Ω→ ∞) of the stochastic process (Kurtz, 1971). These equations were numerically integrated using the fourth-order Runge-Kutta routine of the GSL library (Galassi et al., 2010).

Parameter inference via maximum likelihood

For a time series y1:T of T successive observations and a state-space model with parameter vector θ, the likelihood is given by equation M3. For our stochastic model, the likelihood is analytically intractable and we resorted to an iterated filtering procedure which converges to the ML parameter estimate (θML) to the incidence data (Ionides et al., 2006) (code available upon request). In short, this inference framework only requires: (i) an algorithm to simulate the stochastic model; and (ii) an observation process to link the model-predicted incidence (i.e. the daily number of new hosts entering the infectious class I1) to the daily incidence counts reported in the data set. Following Camacho et al. (2011), we used Gillespie's algorithm for model simulation and a Poisson process, whose reporting rate ρ was also inferred, for observation. We performed log-likelihood profiles in order to check convergence to the maximum likelihood and to calculate 95% confidence intervals (CI95%) for parameter estimates.

A detailed description of the inference procedure can be found in the Supplementary Material of Camacho et al. (2011). In particular, it is shown that, in contrast to the other parameters, inference of the shape parameters kE, kI, kC and kW is computationally too expensive. To tackle this issue, we followed Wearing et al. (2005) who fitted a SEIR model to an influenza epidemic in a boarding school and obtained the best fit for kE = kI = 2. Then, we performed sensitivity analyses on kC and kW and found that, whatever the value of kC, the likelihood was maximized when kW = 1 (results not shown). Finally, we performed a log-likelihood profile on kC and found the maximum at kC = 5. In the remainder of this paper we therefore fix kE = kI = 2, kW = 1 and kC = 5 and obtain the model of Fig. 1B.

Quantities of epidemiological interest

Our detailed description of the different immunological stages allows us to derive several quantities of epidemiological interest. We denote by τ the time elapsed since the start of the infectious period and define the following probabilities:

equation M4, the probability that, at time τ, the host is still infectious,
equation M5, the probability that, at time τ, the host is still temporarily protected against HR thanks to the innate and cellular immunity,
equation M6, the probability that, at time τ, the host is already protected on the long-term against HR thanks to the humoral immunity,
P4(τ) = 1 − P2(τ) − P3(τ), the probability that, at time τ, the host is unprotected and can potentially suffer from HR if re-exposed,

where fIC is the probability density function of τI + τC and similarly for fICW. The probability distributions P1–4 were computed via Monte-Carlo integration.

These probability distributions can be compared with empirical ones obtained from volunteer challenge studies (Carrat et al., 2008) or population surveys during natural infections (Miller et al., 2010; Hung et al., 2010; Chan et al., 2011). In a recent study, Baguelin et al. (2011) fitted a Weibull distribution to a data set consisting of 115 time intervals to seroconversion that were obtained from a serological survey during the second wave of the 2009 A/H1N1 pandemic in the UK (Miller et al., 2010). More precisely, the authors defined the seroconversion interval of each individual as the time taken since symptom onset to reach an hemagglutination-inhibition (HI) titre of ≥32 (Baguelin et al., 2011). Historically, HI assay has been considered to be the gold standard for evaluation of the humoral serum response, with an HI titre of ≥32 considered as a surrogate marker for recent infection during the 2009 pandemic (Hardelid et al., 2010). We investigated the link between seroconversion and efficient protection by comparing the Weibull distribution obtained by Baguelin et al. (2011) to the distribution P3(τ) under our ML parameter estimates. For this purpose we implicitly assume that the time of influenza symptom onset coincides with the start of the infectious period.

In order to assess the impact of HR on the epidemiological dynamics, we define the fraction of the population that was (i) infected at least once (FI.), (ii) reinfected (FII.) and iii) infected only once but remains unprotected (FIS) at the end of the epidemic. We note however that the model presented in Fig. 1B does not allow us to compute these fractions since it only tracks the epidemiological status of individuals (e.g. susceptible, infectious, protected) rather than their infection and immunological histories. For instance, it is impossible to distinguish between those infected only once and those reinfected (they all pass through the same stages), nor between those susceptible who escaped infection and those unprotected who escaped reinfection (they all end in the S class). To tackle this issue, we developed a version of the SEICWH model that allows us to track both the infection and the immunological histories of individuals and used this model to compute FI., FII. and FIS. Since this model is of much higher dimension than that of Fig. 1B, we refer to Text S4 for a more detailed description.

Finally, we define the daily inflow of unprotected hosts equation M7 by counting the number of recovered hosts who lose their temporary immunity conferred by the cellular protection (i.e. leave the C5 compartment) during day d, independently of the outcome of their humoral response.

Exploration

Human communities differ from each other in their contact structure. We seek to characterize the interplay between the immunological and epidemiological dynamics of the SEICWH model for various contact rates (β) ranging from highly to less mixed populations. We keep the immunological parameters constant and equal to those inferred from the 1971 TdC epidemic. We can thus express changes in the contact rate in terms of the more meaningful basic reproduction number (R0 = β/ν).

We then focus on the first post-pandemic influenza seasons and seek to determine under which conditions a subsequent variant to the pandemic strain can break population herd immunity, thus leading to a typical seasonal epidemic. Rapid evolution of the pandemic strain through mutations mainly results in changes of its antigenic properties and/or its transmissibility. Other properties such as the duration of the infectious period could also evolve but we keep them constant for simplicity. As such, evolution of the transmissibility translates into a difference ΔR0 between the basic reproduction number of the mutated variant and that of the pandemic strain. On the other hand, antigenic evolution is modelled through an immune escape factor σ [set membership] [0, 1] which corresponds to the proportion of antigenic properties of the mutated variant that differs from the pandemic strain. For instance, σ = 0 means that both the pandemic and the mutant strain share the same antigenic properties. Finally, we assume that immune escape translates into cross-immunity by reducing the susceptibility against infection by the new variant by a multiplicative factor 1 − σ. As such, individuals who have developed a protective humoral response to the pandemic strain are partially protected against infection by the mutant strain whereas those who escaped infection or remained unprotected at the end of the pandemic season are fully susceptible. As detailed in Text S5, previous empirical and modelling studies suggest a relative increase of the transmissibility ΔR0/R0 [set membership] [0, 1] as well as an immune escape σ [set membership] [0, 0.5] for a post-pandemic variant. In the following, we explore these parameter ranges.

Results

Parameter inference

ML estimates and CI95% for the parameter set are presented in Table 2. Our estimates reveal the exceptional epidemiological context of the 1971 epidemic, in the small and fully isolated TdC community, characterized by a high contact rate among the islanders (R0 = 11.78, CI95% = [7.7–25.5]), as well as a very low level of pre-existing immunity at the beginning of the epidemic (S0/Ω ≈ 98 %, [97–99]), the origin of which we speculate upon in Text S6. ML estimates of the generation time (average time between primary and secondary cases: 3.34 days [2.53–4.7]) and of the reporting rate in the data counts (due to asymptomatic infections and observation errors: 71%, [62–82]) are in good agreement with those previously published (Carrat et al., 2008). Similarly, we find that 17%, [0–51] of the infected islanders did not mount an efficient humoral immune response, which is in the range of the estimates available in the literature (Cox et al., 2004; Miller et al., 2010; Chen et al., 2010). Finally, ML estimates of the duration of the short-term protection (13.37 days [10.37–16.31]) that follows recovery and of the window of susceptibility (2.75 days [0–6.03]) that precedes the establishment of a long-term humoral protection are both in good agreement with the timings of the completion of the CTLs contraction (Cowling et al., 2012) and the peak of neutralizing antibodies (see section “The primary immune response to influenza infection in humans”).

Table 2
Results of the maximum likelihood statistical inference for the 1971 TdC epidemic. ML estimates and 95% confidence intervals for the SEICWH model parameters.

We note that the CI95% of the mean window of susceptibility (ω−1) contains 0 which could suggest a more parsimonious model. However, as we show in Text S7, the broad CI95% of ω−1 is due to a strong correlation with the parameter α (the probability to mount an efficient humoral response). In particular, the lower bound ω−1 = 0 corresponds to values of α~ 50 % that are far below empirical estimates found in the literature (~80–90%). Conversely, we show that fixing α~80–90% leads to a much tighter CI95% for ω−1 that excludes 0. As such, we conclude that despite its broad CI95% the window of susceptibility is justified for the sake of biological realism.

Immunodynamics

The immunodynamics under our ML estimates is summarized in Fig. 3. In particular, we find that the waiting time in the window of susceptibility is exponentially distributed, thus revealing a high level of host heterogeneity in the development of the humoral response. Indeed, among those who do mount a protective humoral immune response (83% of the population), 30% will stay in the window of susceptibility less than one day and 20% more than four days (see Fig. S1, green curve). It also shows that, at the population level, the probability P4(τ) to sample an unprotected individual rapidly peaks to 0.25 three weeks after the date of symptom onset, owing to the window of susceptibility, and then decreases to (1 − α) = 0.17 on a longer time scale due to the lack of humoral protection.

Fig. 3
Dynamics of the immune response, under the SEICWH model, inferred from the 1971 TdC epidemic. At the population level, our framework allows us to reconstruct the proportion of individuals that are infectious (dotted black line), short-term protected thanks ...

Finally, we investigate the link between seroconversion and efficient protection by comparing the cumulative distribution obtained by Baguelin et al. (2011) with our probability distribution P3. Fig. 3 suggests that seroconversion occurs faster (~1 week) and in a slightly greater proportion of infected hosts (87% vs. 83%) than efficient protection. We discuss this discrepancy in section “Immunodynamics model”.

Model fit

In order to better characterize the 1971 two-wave epidemic on TdC, we simulated 105 realizations of the stochastic SEICWH model (Fig. 1B) under the ML estimates (Table 2) and computed the 50 and 95 percentile intervals (PI50,95%) of the distribution conditioned on non extinction. Fig. 2 (upper panel) demonstrates the goodness of fit of the SEICWH model since the data lie in the PI95% while the shape and dynamics of the epidemic are closely captured by the mean predicted incidence with PI50% envelope (we refer to Text S8 for a similar analysis using parameter sets sampled from the CI95%).

In addition, Fig. 2 (lower panel) reveals that although the extinction probability increases at the beginning of the epidemic and during the inter-wave period (i.e. when the transmission chain can be broken due to the low number of infectious hosts), the risk of disease fadeout remains below 5%. By contrast, this risk rapidly increases during the downturn of the second epidemic wave due to depletion of the susceptible pool (i.e. most HR individuals have gained long-term protection). We can thus conclude that, despite the small community size of TdC, the infection and HR dynamics were robust to demographic stochasticity during the 1971 epidemic. Put another way, given the population settings of TdC, the HR wave was not a twist of fate but did have a high probability to occur.

Interplay between the immunological and epidemiological dynamics

Three typical epidemic profiles

The first striking result of the exploration is that three typical epidemic profiles can be distinguished, depending on the value of R0. When the contact rate is high (R0 [greater than or approximate] 5), as on TdC, the epidemic is composed of two waves with two distinct peaks (Fig. 4C). By contrast, at intermediate contact rates (R0 [set membership] [2 − 5]), the epidemic is composed of a single wave with a “long tail” end (Fig. 4B). Finally, when the contact rate is low (R0 [less-than or approximate] 2), the tail disappears so that the epidemic becomes bell shaped (Fig. 4A). These three epidemic profiles arise from the interplay between the immunological and epidemiological dynamics that modulates the effective reproduction number Re(t) throughout the epidemic. For our model, Re(t) = R0(S(t) + W(t))/Ω and as long as Re(t) > 1 the epidemic is increasing.

Fig. 4
Example of the three typical epidemic profiles generated by the SEICWH model, depending on the value of R0. A: R0 = 1.4, B: R0 = 4 and C: R0 = 10. These values reflect the tendency of the contact rate among ...

In the parameter region of high R0 (Fig. 4C), the disease spreads so rapidly that almost the entire population is infected over a short time interval which is similar to the duration of the cellular protection. As a result, many recovered hosts lose their cellular protection simultaneously, leading to an important inflow of unprotected individuals equation M8 shortly after the end of the infection wave. Accordingly, Re rapidly increases so that, in the event that the chain of transmission is maintained until the threshold Re = 1 is reached, HR becomes sustained and a second epidemic wave is observed.

By contrast, in the parameter region of low R0 (Fig. 4A), the disease spreads over a much longer time scale than the immune response so that equation M9 peaks during the downturn of the infection wave. This timing, together with the low R0, help to explain why Re remains below one after the first epidemic peak. In this case, the reinfection wave is not sustained but mainly driven by the infection wave.

Finally, in the parameter region of intermediate R0 (Fig. 4B), the disease spreads slowly enough that the reinfection wave is initially driven by the infection wave while equation M10 is sufficient to maintain the chain of transmission after the end of the infection wave. However, in contrast to the high R0 case, Re remains below one so that the epidemic does not peak again, but subsides in a tail of reinfection.

The HR threshold

As shown in Fig. 5 (upper panel), the fraction of individuals infected at least once during the epidemic (FI.) increases rapidly with R0. In particular, the value of FI. is greater than expected with a SEIR model since the latter does not account for new cases resulting from contact with reinfected hosts (the difference is plotted as a dotted-line).

Fig. 5
Change in the fractions of infected (FI., upper panel), reinfected (FII., middle panel) and unprotected (FIS, lower panel) individuals at the end of the epidemic as a function of R0. Each colour refers to an epidemic profile of Fig. 4 (bell: green, tail ...

Regarding the fraction of the population that is reinfected (FII., middle panel) or that remains unprotected (FIS, lower panel) at the end of the epidemic, we note a qualitative change at intermediate values of R0. Specifically, we define the HR threshold by equation M11. When equation M12, most hosts have closed their window of susceptibility before re-exposure to the virus whereas most of those with deficient humoral response are likely to escape HR until the end of the epidemic and remains unprotected. By contrast, when equation M13, most hosts are likely to be re-exposed to the virus before their window of susceptibility closes whereas most unprotected hosts gain long-term immunity via HR during the second epidemic wave. We numerically calculated equation M14 under the ML estimates of the immunological parameters.

Finally, we contend that this HR threshold should not be confused with the reinfection threshold introduced by Gomes et al. (2004). Although both thresholds indicate important qualitative change of the epidemic dynamics, we show in Text S9 that they have different epidemiological interpretations as well as different dynamical implications.

Implications for the first post-pandemic season

Fig. 6 (upper panels) shows the expected fraction of individuals infected at least once by a mutant during the first post-pandemic season (FI. post-pdm) as a function of the immune escape σ and the relative increase in transmissibility ΔR0/R0, for five different values of R0 in agreement with pandemic scenarios in large populations (Lessler et al., 2007). As R0 increases the pandemic becomes more and more severe so that the expected fraction of protected individuals at the beginning of the post-pandemic season (equation M15) increases. Accordingly, from an evolutionary point of view, it becomes more and more efficient for the mutant to increase its immune escape than its transmissibility in order to invade the population. By contrast, when R0 is close to 1, a mutant antigenically similar to the pandemic strain can invade the population following moderate increase in transmissibility.

Fig. 6
Implication of immunodynamics for the first post-pandemic season. Upper panels: expected fraction of individuals infected at least once by a mutant strain during the post-pandemic season (FI. post-pdm, colour coded), as a function of the relative increase ...

This pattern can be compared to that predicted by a SEIR model, i.e. assuming that all the individuals infected during the pandemic develop an efficient humoral response and are therefore partially protected against the mutant strain. One can show that, at the beginning of the post-pandemic season, the SEIR model underestimates the value of Re for the mutant strain by equation M16, where equation M17 is the fraction of unprotected hosts at the end of pandemic season in the SEICWH model. Fig. 6 (lower panels) reveals a parameter region, below the invasion threshold (Re = 1) of the SEIR model, where the SEICWH model predicts epidemics involving up to 25% of the population. Furthermore, even above this invasion threshold, the epidemic sizes differ by the same order of magnitude as a typical seasonal influenza epidemic (ΔFI. ≈ 5–20%).

Discussion

Immunodynamics model

Our study supports the view that host heterogeneity in the timely development of a protective immunity can explain HR. More precisely, although short lived (innate and cellular) immunity should prevent HR within 2–3 weeks following the primary infection (Cowling et al., 2012), incomplete immune formation and non seroconversion can lead to HR following re-exposure to the same strain on an intermediate and a long time-scale, respectively. These mechanisms provide an explanation to the HR cases reported over 2–5 weeks (Perez et al., 2010; Kim et al., 2010) as well as over several months (Trakulsrichai et al., 2012) during the 2009 A/H1N1 pandemic.

To our knowledge, the present statistical analysis is the first one that attempts to provide joint estimates for the duration of the short-term protection, the duration of the window of susceptibility and the probability to develop a long-term protection. Although strong correlations between ω−1 and α prevent us from identifying these key quantities with tight CI95%, our ML estimate of α is in very good agreement with empirical estimates in the literature. Moreover, we show in Text S7 that fixing α around these empirical estimates leads to much tighter CI95% for ω−1. We have so far assumed that hosts are fully susceptibility to HR while in the window of susceptibility. In Text S10 we show that partial susceptibility can be modelled by means of an extra parameter and has the effect to lengthen the window of susceptibility. However, we also found that this extra parameter suffers from serious identifiability issue and choose not to include it in the present SEICWH model.

Finally, comparison of the immunological dynamics under the best fit model with empirical estimates from the 2009 pandemic in the UK (Fig. 3) suggests that either (i) seroconversion occurred more rapidly during the 2009 A/H1N1 pandemic in the UK than during the 1972 A/H3N2 epidemic on TdC or (ii) it should take a higher HI titre than 32 for efficient protection against HR. The first explanation could be justified by different immunogenic properties between A/H1N1 and A/H3N2 as well as different immuno-genetic background between the UK and TdC populations. However, we note that ML estimates of the immunological parameters for the TdC epidemic are in close agreement with empirical literature (see section “The primary immune response to influenza infection in humans”). Accordingly, we believe that our estimates can reasonably be extended to other human populations and influenza viruses. On the other hand, the second explanation is in good agreement with the results of a recent meta-analysis showing that a HI titre of 32 corresponds to less than 50% reduction in the risk of contracting influenza whereas it takes a titre of ≥100 to decrease this risk to 10% (Coudeville et al., 2010). Similarly, we note that most labs outside UK fix the protective threshold at 40 instead of 32, thus increasing the time to seroconversion while reducing the proportion of seroconverted.

Does HR drive multiple-wave influenza outbreaks?

Our study indicates that HR could drive multiple-wave influenza outbreaks in communities with exceptional contact configurations like schools or isolated settlements. However, we have assumed so far a constant contact rate between infected and susceptible (or unprotected) individuals as well as no prior immunity to the new virus. These assumptions seem justified for the isolated and close-knit TdC community. Indeed, the high attack rate (96%) during the 1971 epidemic (Mantle and Tyrrell, 1973) suggests that those who were infected at the beginning of the infection wave have rapidly been re-exposed while caring for the sick during the inter-wave period, thus initiating the HR wave. By contrast, we expect that in less isolated and better prepared communities, past influenza exposures and vaccination should reduce the number of susceptible individuals at the beginning of the epidemic, thus increasing the HR threshold equation M18. On the other hand, distancing or containment measures should rapidly be implemented as the epidemic progresses, thus considerably mitigating the risk of re-exposure. For instance, school closure could rapidly drive the epidemic to extinction, whereas rapid isolation of suspected cases could efficiently reduce the contact rate of infected host, thus preventing the HR wave.

In line with these scenarios, no reinfection was reported during an influenza outbreak that occurred in a boarding school of 763 boys in 1978 (one year after the re-emergence of A/H1N1) despite the high attack rate (67%) (Anonymous, 1978). In this case, infectious boys were confined to bed and cared by 130 adults, presumably already immune as only one of them reported symptoms. As a result, the epidemic died out after 13 days, while most of the recovered boys would still have benefited from a cellular protection (see Fig. S11). By contrast, in 1924, only 121 (13%) of 904 boys of the Royal Navy School of Greenwich had already been infected after 23 days of disease propagation (presumably because of the high level of prior-immunity acquired since the 1918 pandemic), when 40 new boys were distributed indiscriminately throughout the school (Dudley, 1926). On the 26th day, two of this batch developed influenza, and within a week nine new boys had been infected. Meanwhile the incidence among the old boys, which had previously been on the wane, rose again, and 12 reinfections were reported. Interestingly, a further batch of 40 new boys were introduced just at the end of the epidemic, when the chances of infection must have considerably diminished; six of these were ultimately infected during the period in which 14 cases of infection and 4 cases of reinfection were reported among the old boys (see Fig. S12). This epidemic pattern can simply be explained by an increase of Re due to the simultaneous effect of the replenishment of the susceptible pool and the inflow of unprotected old boys.

On the other hand, our results clearly indicate that HR is not sufficient in itself to generate the multiple-wave outbreak patterns observed during past pandemics in large populations. Indeed, R0 has been estimated around 2 (Lessler et al., 2007) which is below the HR threshold equation M19. In these cases, HR would have only increased the force of infection, and thus the number of infected hosts, by a few percent. Once again, we contend that our simple transmission model willingly ignores many known mechanisms at work in larger and more structured populations such as age-dependency in the contact rate (Mossong et al., 2008) (i.e. heterogeneous mixing) and behavioural changes (Funk et al., 2009). Furthermore, propagation of a new influenza virus over several months must depend on seasonal variations in transmissibility (i.e. change in absolute humidity (Shaman and Kohn, 2009)) and contact rate (i.e. school closing and reopening (Hens et al., 2009)) as well as meta-population coupling (Balcan et al., 2009). Indeed, recent studies suggest that the timing of the first and second waves during the 2009 pandemic influenza was controlled by a combination of these mechanisms (Chao et al., 2010; Shaman et al., 2010). Nevertheless, despite the simplicity of our transmission model, we believe that our qualitative conclusions on the role of HR in large populations (i.e. below the HR threshold) remain valid even including these additional mechanisms.

The necessity to account for immunodynamics and host heterogeneity in epidemiological models

In communities with exceptional contact rates, HR can rapidly drive a second epidemic wave involving not only individuals with deficient humoral response but also those who are re-exposed before their humoral protection becomes efficient. Accordingly, the risk of a wave of reinfection cannot be anticipated without a precise description of the immunodynamics that follows recovery from influenza infection as in the SEICWH model. By contrast, in larger and less mixed populations, HR does not significantly alter the epidemiological dynamics so that a simple SEIR framework should be sufficient to predict or infer (retrospectively) the impact of a new virus in these populations.

On the other hand, one should also bear in mind that the SEIR model overestimates the level of population immunity at the end of the pandemic by assuming that all infected hosts develop a protective humoral response. Accordingly, our results indicate that consideration of host heterogeneity in the humoral response is essential in order to anticipate the impact of a mutant in the post-pandemic era. In addition, this would permit to quantify the fraction of infected hosts that should remain unprotected after a pandemic and could therefore benefit from vaccination in order to boost their humoral response. Although it seems difficult to separate protected from unprotected hosts without individual serological tests, we can nonetheless assume that random vaccination of symptomatic cases should, in principle, increase population herd immunity through cross-immunity to subsequent antigenic variants.

The SEICWH model represents a step forward in the consideration of the immune response, and its heterogeneity among individuals, in epidemiological models. However, further research and refinements could be envisaged to improve its realism. First, reinfected hosts may benefit from T-cell “memory” and be less infectious than infected hosts, and even more often asymptomatic, thus reducing their risk to transmit the disease. Second, host heterogeneity in the development of a protective humoral response could vary depending on the immunogenic properties of each influenza virus and the population under study. For instance, it was recently reported that although 90% of the infected hosts in the age range 16–29 seroconverted during the 2009 pandemic, this proportion decreased to 70% for those aged 50 years and over (Hung et al., 2010). Finally, although we have assumed a life-long humoral protection once in the H stage, the same study conducted during the 2009 pandemic also revealed that 7 and 16% of patients who seroconverted had a decline of antibody titre of 4- and 2-fold, respectively, after one year (Hung et al., 2010). As for the lack of immune response, this additional mechanism could have significant implications for the current post-pandemic era by increasing the effective reproduction number of subsequent nH1N1 variants.

Acknowledgments

We sincerely thank Pr John Edmunds for his assistance in editing the manuscript. AC is supported by the Medical Research Council (fellowship MR/J01432X/1). The computational resources used in this work were funded by an investment grant from the Région Île-de-France through the scientific program DIM MALINF.

Footnotes

[star]This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Appendix A. Supplementary data

The following are the supplementary data to this article:

References

Anonymous Influenza in a boarding school. Br. Med. J. 1978;1:587.
Baguelin M., Hoek A.J.V., Jit M., Flasche S., White P.J., Edmunds W.J. Vaccination against pandemic influenza A/H1N1v in England: a real-time economic evaluation. Vaccine. 2010;28:2370–2384. [PubMed]
Baguelin M., Hoschler K., Stanford E., Waight P., Hardelid P., Andrews N., Miller E. Age-specific incidence of A/H1N1 2009 influenza infection in England from sequential antibody prevalence data using likelihood-based estimation. PLoS ONE. 2011;6:e17074. [PubMed]
Balcan D., Colizza V., Gonçalves B., Hu H., Ramasco J.J., Vespignani A. Multiscale mobility networks and the spatial spreading of infectious diseases. Proc. Natl. Acad. Sci. U.S.A. 2009;106:21484–21489. [PubMed]
Barry J.M., Viboud C., Simonsen L. Cross?Protection between successive waves of the 1918–1919 influenza pandemic: epidemiological evidence from US Army Camps and from Britain. J. Infect. Dis. 2008;198:1427–1434. [PubMed]
Brown D.M., Román E., Swain S.L. CD4 T cell responses to influenza infection. Semin. Immunol. 2004;16:171–177. [PubMed]
Camacho A., Ballesteros S., Graham A.L., Carrat F., ratmann O., Cazelles B. Explaining rapid reinfections in multiple-wave influenza outbreaks: Tristan da Cunha 1971 epidemic as a case study. Proc. R. Soc. B. 2011;278:3635–3643. [PMC free article] [PubMed]
Carrat F., Vergu E., Ferguson N.M., Lemaitre M., Cauchemez S., Leach S., Valleron A.J. Time lines of infection and disease in human influenza: a review of volunteer challenge studies. Am. J. Epidemiol. 2008;167:775–785. [PubMed]
Chan K., To K.K.W., Hung I.F.N., Zhang A.J.X., Chan J.F.W., Cheng V.C.C., Tse H., Che X.Y., Chen H., Yuen K.Y. Differences in antibody responses of individuals with natural infection and those vaccinated against pandemic H1N1 2009 influenza. Clin. Vacc. Immunol. 2011;18:867–873. [PMC free article] [PubMed]
Chao D.L., Elizabeth Halloran M., Longini I.M., Jr. School opening dates predict pandemic influenza A(H1N1) outbreaks in the United States. J. Infect. Dis. 2010;202:877–880. [PubMed]
Chen M.I., Barr I.G., Koh G.C.H., Lee V.J., Lee C.P.S., Shaw R., Lin C., Yap J., Cook A.R., Tan B.H., Loh J.P., Barkham T., Chow V.T.K., Lin R.T.P., Leo Y.S. Serological response in RT-PCR confirmed H1N1-2009 influenza A by hemagglutination inhibition and virus neutralization assays: an observational study. PLoS ONE. 2010;5:e12474. [PubMed]
Coudeville L., Bailleux F., Riche B., Megas F., Andre P., Ecochard R. Relationship between haemagglutination-inhibiting antibody titres and clinical protection against influenza: development and application of a bayesian random-effects model. BMC Med. Res. Methodol. 2010:1–11. [PubMed]
Cowling B.J., Fang V., Nishiura H., Chan K., Ng S., Ip D.K.M., Chiu S.S., Leung G.M., Peiris J.S. Increased risk of noninfluenza respiratory virus infections associated with receipt of inactivated influenza vaccine. Clin. Infect. Dis. 2012;54:1778–1783. [PubMed]
Cox R.J., Brokstad K.A., Ogra P. Influenza virus: immunity and vaccination strategies. Comparison of the immune response to inactivated and live, attenuated influenza vaccines. Scand. J. Immunol. 2004;59:1–15. [PubMed]
Dudley S. 1926. The spread of “Droplet Infection” in Semi-isolateed Communities. Technical Report 111, Oxford.
Fairlie-Clarke K.J., Shuker D.M., Graham A.L. Perspective article: Why do adaptive immune responses cross-react? Evol. Appl. 2008;2:122–131.
Fraser C., Donnelly C.A., Cauchemez S., Hanage W.P., Van Kerkhove M.D., Hollingsworth T.D., Griffin J., Baggaley R.F., Jenkins H.E., Lyons E.J., Jombart T., Hinsley W.R., Grassly N.C., Balloux F., Ghani A.C., Ferguson N.M., Rambaut A., Pybus O.G., Lopez-Gatell H., Alpuche-Aranda C.M., Chapela I.B., Zavala E.P., Guevara D.M.E., Checchi F., Garcia E., Hugonnet S., Roth C., The WHO Rapid Pandemic Assessment Collaboration Pandemic potential of a strain of influenza A (H1N1): early findings. Science. 2009;324:1557–1561. [PubMed]
Funk S., Gilad E., Watkins C., Jansen V.A.A. The spread of awareness and its impact on epidemic outbreaks. Proc. Natl. Acad. Sci. U.S.A. 2009;106:6872–6877. [PubMed]
Galassi M., Theiler J., Jungman G., Gough B., Davies J., Priedhorsky R., Booth M., Rossi F. 2010. GNU Scientific Library Reference Manual.
Gillespie D.T. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 1977;81:2340–2361.
Gomes M.G.M., White L.J., Medley G.F. Infection, reinfection, and vaccination under suboptimal immune protection: epidemiological perspectives. J. Theor. Biol. 2004;228:539–549. [PubMed]
Hardelid P., Andrews N.J., Hoschler K., Stanford E., Baguelin M., Waight P.A., ZAMBON M., Miller E. Assessment of baseline age-specific antibody prevalence and incidence of infection to novel influenza A/H1N1 2009. Health Technol. Assess. 2010;14:115–192. [PubMed]
Hens N., Ayele G., Goeyvaerts N., Aerts M., Mossong J., Edmunds J.W., Beutels P. Estimating the impact of school closure on social mixing behaviour and the transmission of close contact infections in eight European countries. BMC Infect. Dis. 2009;9:187. [PubMed]
Hung I.F.N., To K.K.W., Lee C.K., Lin C.K., Chan J.F.W., Tse H., Cheng V.C.C., Chen H., Ho P.L., Tse C.W.S., Ng T.K., Que T.L., Chan K.H., Yuen K.Y. Effect of clinical and virological parameters on the level of neutralizing antibody against pandemic influenza A virus H1N1 2009. Clin. Infect. Dis. 2010;51:274–279. [PubMed]
Ionides E.L., Bretó C., King A.A. Inference for nonlinear dynamical systems. Proc. Natl. Acad. Sci. U.S.A. 2006;103:18438–18443. [PubMed]
Kim T.S., Ho K.M., Yim K.R., Oh W.S., Chon S.B., Ryu S.W., Yie K., Lee S.J. Three reinfection cases of the pandemic influenza (H1N1 2009) Infect. Chemother. 2010;42:257.
Kurtz T.G. Limit theorems for sequences of jump Markov processes approximating ordinary differential processes. J. Appl. Probabil. 1971:344–356.
Lessler J., Cummings D.A.T., Fishman S., Vora A., Burke D.S. Transmissibility of swine flu at Fort Dix, 1976. J. R. Soc. Interface. 2007;4:755–762. [PubMed]
Mantle J., Tyrrell D.A. An epidemic of influenza on Tristan da Cunha. J. Hyg. 1973;71:89–95. [PubMed]
McGill J., Heusel J.W., Legge K.L. Innate immune control and regulation of influenza virus infections. J. Leukoc. Biol. 2009;86:803–812. [PubMed]
Medical Department of the Local Government Board 1919. Report of the Medical Department of the Local Government board for 1918–19. Technical Report, London.
Miller E., Hoschler K., Hardelid P., Stanford E., Andrews N., Zambon M. Incidence of 2009 pandemic influenza A H1N1 infection in England: a cross-sectional serological study. Lancet. 2010;375:1100–1108. [PubMed]
Ministry of Health 1920. Reports on Public Health and Medical Subjects. Technical Report 4, London.
Mossong J., Hens N., Jit M., Beutels P., Auranen K., Mikolajczyk R., Massari M., Salmaso S., Tomba G.S., Wallinga J., Heijne J., Sadkowska-Todys M., Rosinska M., Edmunds W.J. Social contacts and mixing patterns relevant to the spread of infectious diseases. PLoS Med. 2008;5:e74. [PubMed]
Perez C.M., Ferres M., Labarca J.A. Pandemic (H1N1) 2009 reinfection, Chile. Emerging Infect. Dis. 2010;16:156–157. [PubMed]
Shaman J., Goldstein E., Lipsitch M. Absolute humidity and pandemic versus epidemic influenza. Am. J. Epidemiol. 2010;173:127–135. [PubMed]
Shaman J., Kohn M. Absolute humidity modulates influenza survival, transmission, and seasonality. Proc. Natl. Acad. Sci. U.S.A. 2009;106:3243–3248. [PubMed]
Tamura S.i., Kurata T. Defense mechanisms against influenza virus infection in the respiratory tract mucosa. Jpn. J. Infect. Dis. 2004;57:236–247. [PubMed]
Taubenberger J.K., Morens D.M. 1918 Influenza: the mother of all pandemics. Emerging Infect. Dis. 2006;12:15–22. [PubMed]
Trakulsrichai S., Watcharananan S.P., Chantratita W. Influenza A (H1N1) 2009 reinfection in Thailand. J. Infect. Public Health. 2012;5:211–214. [PubMed]
Wearing H.J., Rohani P., Keeling M.J. Appropriate models for the management of infectious diseases. PLoS Med. 2005;2:e174. [PubMed]
Woodland D.L. Cell-mediated immunity to respiratory virus infections. Curr. Opin. Immunol. 2003;15:430–435. [PubMed]