About this Journal Submit a Manuscript Table of Contents
ISRN Biomathematics
Volume 2012 (2012), Article ID 471653, 23 pages
http://dx.doi.org/10.5402/2012/471653
Research Article

Dynamics of Single-City Influenza with Seasonal Forcing: From Regularity to Chaos

Centre for Nutrition Modelling, Department of Animal and Poultry Science, University of Guelph, Guelph, ON, Canada N1G 2W1

Received 7 September 2011; Accepted 17 October 2011

Academic Editors: C. B-Rao and I. Rogozin

Copyright © 2012 John H. M. Thornley and James France. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Seasonal and epidemic influenza continue to cause concern, reinforced by connections between human and avian influenza, and H1N1 swine influenza. Models summarize ideas about disease mechanisms, help understand contributions of different processes, and explore interventions. A compartment model of single-city influenza is developed. It is mechanism-based on lower-level studies, rather than focussing on predictions. It is deterministic, without non-disease-status stratification. Categories represented are susceptible, infected, sick, hospitalized, asymptomatic, dead from flu, recovered, and one in which recovered individuals lose immunity. Most categories are represented with sequential pools with first-order kinetics, giving gamma-function progressions with realistic dynamics. A virus compartment allows representation of environmental effects on virus lifetime, thence affecting reproductive ratio. The model's behaviour is explored. It is validated without significant tuning against data on a school outbreak. Seasonal forcing causes a variety of regular and chaotic behaviours, some being typical of seasonal and epidemic flu. It is suggested that models use sequential stages for appropriate disease categories because this is biologically realistic, and authentic dynamics is required if predictions are to be credible. Seasonality is important indicating that control measures might usefully take account of expected weather.

1. Introduction

1.1. Background and Objectives

Both seasonal influenza and influenza epidemics continue to cause, despite prophylactic and therapeutic efforts, considerable morbidity, mortality, and financial loss across the world (e.g., [1, 2]). Further concern is generated by connections between human influenza with avian influenza ([3]; see below), which, in its own right, gives rise to economic damage and distress [4]. Models provide a method of summarizing current ideas about mechanisms of disease development and propagation, understanding the contributions from the different processes, and exploring possible consequences and interventions [57]. Because of the world-wide relevance of influenza, human and avian, there continue to be many publications on the topic. The recent outbreak of H1N1 swine influenza, and especially the varied and conflicting prognoses from experts, reveals how thinly based much of our knowledge still is.

Our paper is partly a review, partly sets out in detail a particular approach to influenza modelling, and partly has some original content on seasonality, forcing, and the resulting predictions. The particular objectives are to construct a model with multiple sequential pools in various disease categories (Figure 1) so that the underlying biological dynamics is realistic, to validate the model by applying it to data describing an influenza epidemic (Figure 5), to suggest a possible mechanism for the direct effects of environment (season) on influenza dynamics (Figure 6), and last to apply various levels of environmental forcing to demonstrate how the model, without reparameterization, can give rise to a wide range of dynamics, ranging from regular (e.g., twice yearly, annual, biennial, etc.) epidemics to chaos, where sometimes predictions resemble some of the great influenza pandemics (Figures 7 and 8).

471653.fig.001
Figure 1: Model scheme. The 32 state variables are denoted by n + subscript. All pools other than the virus pool (𝑛vir) and the dead-from-influenza pool (𝑛ded) suffer from a natural mortality rate of μ (Table 1, where parameters are listed). Many categories (e.g., “Infected”) consist of sequential pools. The four infected pools all have (incubation) rate constants of 𝑘inf. Similarly, 𝑘sic applies to the four nonhospitalized clinically sick pools, 𝑘hos to the three hospitalized clinically sick pools, 𝑘asy to the seven asymptomatic pools, and 𝑘del to the ten pools in the loss-of-immunity delay pipe. 𝑓asy and 𝑓hos are the fractions of the fluxes between 𝑛inf1 and 𝑛inf2, and 𝑛sic1 and 𝑛sic2, which are asymptomatic and hospitalized. 𝜇sic and 𝜇hos are the death rates (caused by flu) of the clinically sick and hospitalized pools.

Before proceeding with this agenda, a brief review of some of the issues and modelling approaches relevant to influenza is given, and our contribution is discussed in this wider context.

1.2. Brief Review of Influenza Issues and Modelling Approaches
1.2.1. Avian Influenza

One of the reasons for the ongoing concern with influenza is the continuing threat from avian flu. Avian influenza is endemic in many parts of the world and, from time to time, there are outbreaks of highly pathogenic strains, such as the recent outbreak of the H5N1 strain of the virus [8]. The virus can be transferred to other animals, including horses, pigs, and humans, frequently with fatal consequences [9]. Because of antigenic drift (a high rate of immunologically significant mutations) and shift (reassortment between different strains of influenza within a single host), coupled to other factors such as loss of immunity, there is a continuing threat to the human population [10, 11]. This threat depends on possible changes in the to-date limited capacity of the highly pathogenic avian strains for human-to-human transmission [8, 9]. This risk makes it important that key factors determining virus evolution, epidemic occurrence, and cross-species transmission are well understood, so that effective strategies for containment and control might be designed [1214].

There are also models dealing with the avian-human influenza nexus. For example, [14] presents an ordinary differential equation model which combines an SI (susceptible, infected) avian model with an SIR (susceptible, infected, resistant) human model. They suggest that measures such as both extermination (of avians) and quarantine (for humans) could be needed to avoid a pandemic of influenza.

1.2.2. Two-Strain Influenza and Influenza in a Single Person

An interesting paper [15] extends what might be called the complex models of simple influenza by presenting a simple deterministic model of a more complex situation: they treat the dynamics of two-strain influenza, focussing on competition and cross-immunity. Isolation period and crossimmunity are critical parameters. Some of their results are similar to those reported here (with variable interepidemic periods from 2 to 10–13 years), although the model and mechanisms are quite different.

At a more detailed level, the cell level in a single person, a model of the immune response to the influenza virus which treats innate and adaptive immunity has been proposed [16]. The model has 10 ordinary differential equations, representing interferon, T-cells, killer cells, antibodies, and other states. They explore the impact of initial viral load on disease progression. When this is small, the disease progresses asymptomatically. The model builds on [17].

1.2.3. Network Models and Stochastic Simulations

Some recent influenza models are based on networks and stochastic simulations [11, 1822]. The model of [23] is of considerable interest: it includes an individual level (age, treatment, vaccination status) as well as a community level (household, workplace, supermarkets, schools, etc.). Within a city, most contacts occur in a few locations. Interventions at these locations can be expected to be more effective than less-targeted interventions. Isolation and quarantine (e.g., [24]) are possible treatments, as is antiviral use for both prophylaxis and therapy. A stochastic calculation with a half-day time step is applied. Their model gives valuable quantitative indications of how epidemics/pandemics may be prevented or controlled. These highly articulated models may be well suited to predictions of outcomes from interventions for particular situations. However, stochastic models are highly demanding, of computing resources, high-resolution data on populations, contacts, transmission, age-specific characteristics, and so forth. Moreover, their predictions are often highly sensitive to initial values and settings.

1.2.4. Deterministic Models

Deterministic models of influenza are also numerous, and the model presented here is of this type. While biological variability is a reality, there is a feeling amongst some modellers, especially those with a background in engineering or physics or applied mathematics, that a deterministic approach gets closer to the real science, to understand what is going on, than a stochastic approach. Reflecting this are many quotations, such as “God does not play dice,” and “if you need statistics, then do a better experiment.” Deterministic models are easier to build, easier to understand and use, easier to falsify, and easier for tracing cause and effect. With a deterministic model, a clearly labelled box diagram tells the reader most of what he needs to know. In the next two paragraphs, we refer to two examples of the genre. Although many models have the potential to represent directly seasonal effects of environment, we have not found one that actually does this. Some models represent places where people assemble and interact such as schools and clinics which affect contact rates: indirectly, these can represent a seasonality. But here, effects of weather variables on model parameters are represented directly.

The model of [25] is a typical example of a simple deterministic influenza model. The model has seven significant state variables and ordinary differential equations. From our point of view, their use of a single stage to represent each category (e.g., latent, infectious, asymptomatic, hospitalized, etc.) gives a biologically unreasonable representation of the transit time distribution through that category. Since here we represent day-to-day effects of weather on influenza dynamics, it is important that the model should represent the dynamics of influenza realistically. This means using several serial stages for progress through each disease category.

A typical large deterministic influenza model is in [26]. This comprises over 1000 differential equations and allows for many demographic and clinical parameters (such as risk, age, four levels of sickness, treated or untreated at home, and treated or untreated in hospital) so that it is useful in planning. Their model does employ multiple stages for the different disease categories, and therefore the model is able to be much more dynamically realistic than that of [25]. The model has been used to explore the consequences of pharmaceutical and nonpharmaceutical interventions by [27]. The interested reader might find it useful to start with [27, Figure 7] and proceed to [26, Tables 2, 3, and 4]. The model has not been applied to direct seasonal effects, although effects such as school closures giving a decrease in contact rates are included. There is no mention of chaos resulting from such forcing.

[28] addresses generally for SIR and SEIR models some of the issues covered here: the inclusion of more realistic distributions; the destabilizing consequences of this so that lower levels of forcing are required to give chaos; the conclusion that the assumptions made in formulating the model have a major impact on its dynamical properties.

1.2.5. Simplicity versus Complexity in Influenza Modelling

There are numerous other examples of models of the deterministic ordinary differential equation type. Many parameters of these models are uncertain; predictions can be highly sensitive to initial conditions. The model may be applied as a large regression equation, irrespective of the known biological lacunae and the well-documented dangers of such procedures [29]. Pertinent to this dilemma, it has been remarked that “phenomenological approaches are deficient in their lack of attention to underlying processes; individual-based models, on the other hand, may obscure the essential interactions in a sea of detail [30]. The challenge then is to find ways to bridge these levels of description, ….” Others (e.g., [31]) have argued the importance of modelling the dynamics of influenza surveillance data, in order to provide early predictions of epidemic events; to this end, they apply purely statistical methods. Cogently, it has been suggested that “as a general policy in preparing for an outbreak of a disease whose parameters are not yet known, it would be better to use a general compartment model using relatively few parameters and not depending critically on the particular as yet unknown setting” [32]. We concur with this view. Such models are easier to construct and explore. They are better suited for elucidating general properties of these systems, as is done here.

1.2.6. Seasonality

Seasonality has long been implicated in influenza incidence and severity, although the basis for this is not understood [33]. Also, it has arguably been given little detailed attention by epidemiologists. Seasonality is a significant factor in mortality from several causes including influenza in temperate countries, with more people dying in the winter (approximately, November to March in the northern hemisphere) than in the summer [34, 35]. The contribution of influenza to these excess deaths is disputed [35, 36], as vaccination against influenza protects against deaths from other conditions [37]. However, [35] summarizes with “our findings are compatible with the hypothesis that the cause of winter-season excess mortality is singular and is most likely to be influenza.” This conclusion agrees with an earlier European study [38]. Any model which covers many community levels (e.g., household, workplace, supermarkets, schools, as in [23]), offers many possibilities for applying seasonality by altering mixing patterns. Often seasonal forcing is represented empirically, by adding a sinusoid to the infection rate parameter [39, 40]. This approach to forcing leads to annual, biennial, and multiple cycles including chaos. It has been suggested that large seasonal oscillations in incidence can result from an amplification of very small seasonal changes in influenza transmission [41]. Large amplification occurs when the driving frequency is close to the natural frequency of the unforced system. A two-state variable SIRS model (S + I + R = constant) is applied. A two compartment model with linear transfers (giving rise to negative exponentials) can only give a limited representation of the biological dynamics. They do not explore their model other than to support their suggestion of dynamical resonance and make no mention of chaos.

1.3. Current Model and Its Contributions

Here a simple model for epidemic influenza in a single city with seasonal forcing is constructed and evaluated. We are not aware of existing work which treats directly seasonal effects (but see [42], where the contact rate is reduced by a factor of ten for the 6-month nonepidemic season). The focus is on the essential biology of the problem using traditional scientific reductionism. The model is of the compartmental deterministic type with homogeneous mixing (but see [30]), and is without age or any other non-disease-status stratification. Various categories are represented, but, recognizing the importance of mechanistically realistic dynamics, and at variance with usual practice in SIR models, each category is represented by three or more stages. The reason for using several sequential stages (spelt out in Appendix A) to represent a given category is that this allows a more credible gamma-function progression through the category [43]. With a single-stage category, the most probable time a person spends in that category is zero, which is hardly biologically defensible although it is widely applied. Because dynamics is so important when considering influenza and especially its interaction with seasonal forcing, it is necessary to use the more realistic multistage categories. This adds to the size of the model, but importantly, it does qualitatively change the dynamics (Appendix A). Essentially, our model is of the SIRS (susceptible, infected, resistant, susceptible) type but with most categories represented by sequential stages: for example, the I (infected) category in Figure 1 is broken down into four sequential stages.

There are many ways in which seasonality could impinge on influenza dynamics. We choose one of the simplest. A virus compartment allows effects of environment to be represented on virus lifetime, which might be an important environmental forcing mechanism. Some simulations are presented for the unrealistic situation in which recovered individuals are immune for life, in order to illustrate the basic characteristics of the model. Also, simulations are given for the more realistic situation where immunity is lost over a few years. For this latter situation, seasonal forcing gives rise to an unexpectedly wide range of pertinent dynamics, including regularly spaced epidemics from two per year through one per year to one epidemic every several years (two or more), sometimes with slightly chaotic spacing but sometimes regularly spaced but with chaotic amplitudes, and sometimes quite chaotic in both spacing and amplitude.

Our main objective is to show how, in a mechanistically-oriented model with credible biological assumptions and minimal parameter adjustment to obtain specific outcomes, seasonal forcing functions of different magnitudes can constrain, entrain, and amplify the natural rhythms of influenza, giving rise to a wide range of epidemic/pandemic patterns, from biennial, annual, at intervals of several years, and chaotic. Our novel contributions are first to suggest an explicit direct mechanism for the effects of weather on influenza dynamics and then give simulations that show that such mechanisms can have a profound but realistic effect on the dynamics of influenza epidemics. This suggests that the approach could be an important (but hitherto neglected) part of influenza models and planning tools.

2. Methods and Modelling

2.1. Model Scheme

The scheme is drawn in Figure 1. State variables are denoted by n + subscript. There is no age stratification. There are eight categories of persons: susceptible (sus); infected, which is considered as four sequential stages (inf𝑗,𝑗=1,,4); clinically sick, also considered as four stages of sickness (sic𝑗,𝑗=1,,4); recovered or immune (rec); asymptomatic (asy), which branches off after the first infected stage and has seven stages before recovery is achieved; hospitalized (hos) or isolated, which branches off after the first clinically sick stage and has three stages before recovery is achieved; dead from influenza (ded). Note that clinically sick persons, whether hospitalized or not, may die from influenza, or recover. There is a delay 𝜏imm (day) during which recovered persons can lose immunity and return to the susceptible compartment—persons in this delay pipe, represented by ten sequential stages, are denoted by 𝑛del1,𝑛del2,, to 𝑛del10. Infectious persons give rise to virus particles (𝑛vir), which, while they may be inactivated or killed, can give rise to further infection events.

We choose, although this is not the usual procedure in influenza models, to have a virus pool (𝑛vir). The reader may wonder why? First, we are fairly sure that such a pool must exist. Second, with a virus pool it becomes easier to think in concrete terms of the effects of weather variables (air temperature, relative humidity, wind speed, and radiation) on the virus, for example, its viability and longevity ((15), (18)). An alternative to having a virus pool (𝑛vir) would be to assume that weather affects transmission rate β directly. In the equations to follow, transmission rate β and virus pool 𝑛vir always appear multiplicatively ((2), first equation of (3), (20)). Thus, it perhaps makes little difference whether we have a virus pool which can be modified by weather, or no virus pool and simply assume that weather affects transmission coefficient. We take the view that transmission is a multistage process and that the components of weather may impinge on different parts of this process. It may then be helpful to have an explicit virus pool. Note that our transmission rate β has units of virus−1 day−1 (rather than the customary day−1) (Table 1). Also, in the expression for the basic reproductive ratio (20), β is divided by the virus mortality rate. In fact, the first term on the right side of (20) can be viewed as a traditional transmission rate with units of day−1, and the modulation of virus mortality ((15), (18)) may be considered as modulating the traditional transmission rate.

tab1
Table 1: State variables, parameters and definitions. Initial values of state variables, principal parameter values and definitions of principal variables. See Figure 1 for state variables. Relevant equation numbers are given.

Initial values and parameters are listed in Table 1, although some parameters which are only used once are defined in or after the equation where they appear.

All routes from susceptible to recovered (Figure 1) pass through eight pools. We have assigned the same value of 2 day−1 to the four rate constants 𝑘inf,𝑘asy,𝑘hos, and 𝑘sic (Table 1). This gives a mean transit time from susceptible to recovered of 4 days (8/2) with simple gamma-type distributions of transit times applying to the whole path and its components (last paragraph of Methods section; Appendix A; e.g. [44], pp. 818–822).

Because of the importance of the assumption of sequential pools for giving biologically realistic dynamics (whereas the traditional assumption of a single pool gives biologically unacceptable dynamics [43]), a discussion of sequential pools in relation to the gamma function is given in Appendix A. In Appendix A, it is shown that the use of two sequential pools give qualitatively different dynamics than a single pool also that three sequential pools gives again a qualitatively different result than two pools; with three or more sequential pools, the dynamics only changes quantitatively. Although observation and data clearly support the existence of a minimum time span being required to traverse a given clinical category, which can be represented by sequential pools as is done here, measurements in this area are extremely difficult. Where the data do not speak clearly, for example, as to whether one should use 3 or 4 pools, or 7 or 8, we have made simplifying and convenient assumptions in order that we could proceed with the calculations (see Section 2.3).

[25] uses a model of similar type but with a simpler structure. The authors employ least squares to estimate most of their model's parameters for the spring and autumn waves of the 1918 influenza epidemic (their Table 1). They find substantial differences between some of the parameters for the spring and autumn epidemics, which may raise questions about what such parameters describe. Our parameters (Table 1) are generally compatible with those obtained and applied in [25], allowing for differences in model structure.

All state variables (Figure 1) and all terms in the differential equations scale linearly with population size and the size of the virus pool, so that solutions (expressed proportionately) are independent of population (𝑛city, Table 1).

2.1.1. Some Definitions

First, define some totals in terms of the state variables (Figure 1):𝑛inf=𝑛inf1+𝑛inf2+𝑛inf3+𝑛inf4,𝑛sic=𝑛sic1+𝑛sic2+𝑛sic3+𝑛sic4,𝑛hos=𝑛hos2+𝑛hos3+𝑛hos4,𝑛asy=𝑛asy2+𝑛asy3++𝑛asy7+𝑛asy8,𝑛del=𝑛del1+𝑛del2++𝑛del9+𝑛del10,𝑛liv=𝑛sus+𝑛inf+𝑛sic+𝑛hos+𝑛asy+𝑛rec+𝑛del,𝑛fer=𝑛liv=𝑛city,𝑛iash=𝑛inf+𝑛asy+𝑛sic+𝑛hos.(1) The notation of the first five equations is obvious. Total live population 𝑛liv is given by the sixth equation and is assumed equal to the fertile population 𝑛fer (seventh equation); both of these are equal to the city population, 𝑛city ((2) and following paragraph). Finally, total infected population 𝑛iash is obtained by adding together the inf, asy, sic, and hos categories. All the variables in (1) can be turned into fractions by division by 𝑛city.

State variables are next treated by categories and pools.

2.1.2. Susceptibles

The differential equation for 𝑛sus isd𝑛sus(𝑡)d𝑡=𝜇𝑛fer+𝑘del𝑛del10+𝐼ded𝜇𝑛sus𝛽𝑛sus𝑛vir𝑛liv.(2) The two principal inputs are first from births with the number of fertile persons 𝑛fer given by (1) and second from the output of compartment 10 of the delay pipe (Figure 1) where people recovered from influenza are slowly losing their immunity. μ is the natural birth (and death) rate (Table 1). It is assumed that all births are free of infection and are without immunity. A small additional input to the 𝑛sus pool is included, 𝐼ded. This is the death rate caused by influenza (Figure 1). 𝐼ded is given by (16). This is done so that the live person number, 𝑛liv (1), remains constant at value 𝑛city (Table 1). This makes the results easier to check as true steady states can be obtained. It is of no biological significance as over an epidemic, deaths from influenza are typically less than 1% of those from natural mortality (the −𝜇𝑛sus term in (2)). Also, the instantaneous death rate from influenza, 𝐼ded (16), is transient, and even at its maximum value, is usually less that the death rate from natural mortality. This assumption of equal birth and death rates was also made in [28].

The two outputs (negative terms) are natural death and infection, the latter giving a transfer to the first infected pool. The natural death rate μ is assumed to be the same as the birth rate (Table 1). β is the infection transmission rate (Table 1). This is multiplied by the number of susceptibles 𝑛sus, the virus quantity 𝑛vir and is divided by the number of live persons 𝑛liv (1). All terms in (2) scale equally with population size.

2.1.3. Infecteds

The differential equations for the four sequential infected pools ared𝑛inf1=d𝑡𝛽𝑛sus𝑛vir𝑛liv𝑘inf𝑛inf1𝜇𝑛inf1,d𝑛inf2=d𝑡1𝑓asy𝑘inf𝑛inf1𝑘inf𝑛inf2𝜇𝑛inf2,d𝑛inf3d𝑡=𝑘inf𝑛inf2𝑘inf𝑛inf3𝜇𝑛inf3,d𝑛inf4d𝑡=𝑘inf𝑛inf3𝑘inf𝑛inf4𝜇𝑛inf4.(3) In the first equation, the input of infecteds is the last term in (2). All pools suffer equally from natural death at rate μ. The output term 𝑘inf𝑛inf1 from the inf1 pool is partly asymptomatic with fraction 𝑓asy (Figure 1, Table 1; (7)), the remainder entering the inf2 pool (second equation, first term). The third and fourth equations are straightforward. Rate constant 𝑘inf is 2 day−1, such that with four pools the mean time from infection to clinical sickness is 4/𝑘inf=2 days [45]. This gives a gamma-distributed lag for the overall exit time from the fourth pool (Appendix A, (A.4), Figure 9.)

2.1.4. Nonhospitalized Clinically Sick

Persons who become clinically sick enter the sic1 pool (Figure 1). A fraction 𝑓hos (Table 1) of these may be hospitalized or isolated, the remainder continuing to recover from the illness at home. Some of the clinically sick (whether in hospital or not) will die from influenza. The differential equations ared𝑛sic1d𝑡=𝑘inf𝑛inf4𝑘sic𝑛sic1𝜇+𝜇sic𝑛sic1,d𝑛sic2=d𝑡1𝑓hos𝑘sic𝑛sic1𝑘sic𝑛sic2𝜇+𝜇sic𝑛sic2,d𝑛sic3d𝑡=𝑘sic𝑛sic2𝑘sic𝑛sic3𝜇+𝜇sic𝑛sic3,d𝑛sic4d𝑡=𝑘sic𝑛sic3𝑘sic𝑛sic4𝜇+𝜇sic𝑛sic4.(4) These equations are similar to (3), but the natural death rate μ is augmented by flu-induced deaths at rate 𝜇sic (Table 1). The time between the onset of sickness and recovery is similarly (to the 𝑛inf pools above) gamma-distributed (Appendix A, Figure 9, (A.4)).

The rate at which individuals are admitted to hospital, 𝐼hos2 (Figure 1), is𝐼hos2=𝑓hos𝑘sic𝑛sic1,(5) with units of persons day1. This is often a recorded statistic.

The total flu-related death flux is (input from all four sick pools to the dead-from-influenza 𝑛ded box of Figure 1; units: persons day1)𝐼sicded=𝜇sic𝑛sic1+𝑛sic2+𝑛sic3+𝑛sic4.(6)

2.1.5. Asymptomatic Infecteds

The differential equations for these seven pools (Figure 1) ared𝑛asy2d𝑡=𝑓asy𝑘inf𝑛inf1𝑘asy𝑛asy2𝜇𝑛asy2,d𝑛asy3d𝑡=𝑘asy𝑛asy2𝑘asy𝑛asy3𝜇𝑛asy3,d𝑛asy8d𝑡=𝑘asy𝑛asy7𝑘asy𝑛asy8𝜇𝑛asy8.(7) Asymptomatic infecteds are fed from (3) (first equation, 2nd term on the right side). Seven pools are employed so that the path from susceptibles to recovered (Figure 1) traverses eight pools in total.

2.1.6. Hospitalized Clinically Sick

These represent hospitalized or isolated or specially treated clinically sick. The differential equations ared𝑛hos2d𝑡=𝑓hos𝑘sic𝑛sic1𝑘hos𝑛hos2𝜇+𝜇hos𝑛hos2,d𝑛hos3d𝑡=𝑘hos𝑛hos2𝑘hos𝑛hos3𝜇+𝜇hos𝑛hos3,d𝑛hos4d𝑡=𝑘hos𝑛hos3𝑘hos𝑛hos4𝜇+𝜇hos𝑛hos4.(8) The input term 𝑓hos𝑘sic𝑛sic1 represents hospital admissions (5), a recorded statistic, which may be useful for comparing with data. The flu-induced death rate 𝜇hos (Table 1) is assumed to be the same as that for the nonhospitalized clinically sick, 𝜇sic (4). This may be justified in that the fraction 𝑓hos taken into hospital is more ill, but they then receive better care.

Total flu-related death rate from the three hospitalized clinically sick pools is𝐼hosded=𝜇hos𝑛hos2+𝑛hos3+𝑛hos4.(9)

2.1.7. Recovered and Immune

The state variable 𝑛rec is governed by the equationd𝑛recd𝑡=𝑘asy𝑛asy8+𝑘sic𝑛sic4+𝑘hos𝑛hos4𝜇𝑛rec𝑘rec𝑛rec.(10) The positive (input) terms are from the last equations of (7), (4), and (8). Death occurs at its natural rate (μ). The rate constant 𝑘rec is set to zero if it is assumed that immunity is not lost; otherwise it is set to a nonzero value (e.g., 1 day−1, Table 1, the precise value is unimportant as long as it is of order of 1 day−1 or more). A non-zero 𝑘rec causes rapid entry into the delay sequence of pools (Figure 1) which leads to loss of immunity.

2.1.8. Delay Pipe Representing Loss of Immunity

Loss of immunity might arise from loss of immunological memory, or from drift and shift in the antigenic character of the virus [11]. This process is represented by ten sequential compartments giving a gamma function delay (Appendix A; also e.g., [44], pp. 818–822). If 𝜏imm (day) is the time period during which immunity is lost (three years is assumed; Table 1), and 𝑘del (day−1) is the rate constant out of each of the ten compartments, then𝑘del=10𝜏imm.(11) The standard deviation in the exit times is (10)/𝑘del and the coefficient of variation is 1/10 (=standarddeviation/𝜏imm) (A.7). The differential equations for the pools ared𝑛del1d𝑡=𝑘rec𝑛rec𝑘del𝑛del1𝜇𝑛del1,d𝑛del2d𝑡=𝑘del𝑛del1𝑘del𝑛del2𝜇𝑛del2,d𝑛del9d𝑡=𝑘del𝑛del8𝑘del𝑛del9𝜇𝑛del9,d𝑛del10d𝑡=𝑘del𝑛del9𝑘del𝑛del10𝜇𝑛del10;(12) The first term on the right of the 1st equation (𝑘rec) is the last term in (10); the second term (𝑘del) is transfer to the next compartment; the last term (μ) is the natural death rate. The 2nd term on the right of the last equation (𝑘del) is transfer into the susceptible pool (2).

2.1.9. Virus Pool

It is assumed that virus production, 𝑠vir, which is a surrogate for “infectious strength”, is calculated with𝑠vir(𝑡)=𝑘vir𝑤inf1𝑛inf1+𝑤inf2𝑛inf2+𝑤inf3𝑛inf3+𝑤inf4𝑛inf4+𝑤sic1𝑛sic1+𝑤sic2𝑛sic2+𝑤sic3𝑛sic3+𝑤sic4𝑛sic4+𝑤asy2𝑛asy2+𝑤asy3𝑛asy3++𝑤asy8𝑛asy8+𝑤hos2𝑛hos2+𝑤hos3𝑛hos3+𝑤hos4𝑛hos4.kvir=200virusunits(weightedinfectedperson)1day1.(13)𝑘vir is a production rate constant which is applied to weighted infected persons. The value above gives a basic reproductive ratio 𝑅0 (20) of 4.78 (see Section 2.3). The weighting factors 𝑤 of (13) are (0𝑤1)𝑤inf1=0,𝑤inf2=0.5,𝑤inf3=1,𝑤inf4𝑤=1;sic1=1,𝑤sic2=1,𝑤sic3=1,𝑤sic4𝑤=0.5;asy2=0.2,𝑤asy3=0.4,𝑤asy4𝑤=0.4,asy5=0.4,𝑤asy6=0.4,𝑤asy7𝑤=0.2,asy8𝑤=0;hos2=0.2,𝑤hos3=0.2,𝑤hos4=0.1.(14) The precise values do matter but they are not important, so long as the values are reasonable. In the usual SIR model there would not be multiple stages as here and any weighting factor would be implicitly unity. In (14) infectivity increases a half day after infection, reaches a maximum, and then decreases as recovery takes place. It is stated in [47], with reference to [46], that “the standard pattern of an influenza A virus in adults is characterized by an exponential growth of virus titre, which peaks 2 to 3 days postinfection (DPI), followed by an exponential decrease until it is undetectable after 6 to 8 DPI.” 𝑠vir of (13) is the input to the virus pool (15).

The outputs from the virus pool (15) are death by natural mortality (𝑘vir,dnm, Table 1) and death induced by a suboptimal environment, represented by rate constant 𝑘vir,env (18). The value of 𝑘vir,dnm corresponds to a half-life of free virus of about 3 h [47]. No extra virus death process is ascribed to the contact/infection process with persons (whether susceptible or immune). Thus, when the basic reproductive ratio is calculated in (20), this is proportional to β, and there is no β in the denominator. Such virus death processes are subsumed in the natural mortality term 𝑘vir,dnm. The differential equation for the virus pool 𝑛vir isd𝑛vird𝑡=𝑠vir𝑘vir,dnm+𝑘vir,env𝑛vir.(15)

2.1.10. Flu-Related Dead Pool

The inputs to this pool are from the sick pools with (6) and (9), so the total input to the 𝑛ded pool (Figure 1) is𝐼ded=𝐼sicded+𝐼hosded.(16) There are no outputs. Therefore,d𝑛dedd𝑡=𝐼ded,𝑡=0,𝑛ded=0.(17) Note that, in order to maintain the total population constant as a mathematical convenience, the birth rate of susceptibles is augmented by the death-from-influenza rate (2). As explained after (2), this has a negligible effect on the performance of the model.

2.1.11. Environmentally induced Virus Mortality 𝑘vir,env and Seasonal Forcing

The environmentally dependent function 𝑘vir,env in (15) (Figure 1) is assumed here to depend only on daily mean values of air temperature, 𝑇air. Possible influences of relative humidity [33], radiation, or wind speed are ignored but could be similarly treated. In a study of relative humidity and temperature on virus transmission, [48] remark that “although the seasonal epidemiology is well characterized, the underlying reasons for predominant wintertime spread are not clear.” The rate constant for environmentally induced death, 𝑘vir,eny (d1; Figure 1) is written as𝑘vir,env=𝑠env𝑘Tair𝑇air𝑇air0𝑇+ABSair𝑇air02𝑞𝑇,𝑠env=0(or1),𝑘Tair=0.1C2d1,𝑇air0=13C,𝑞𝑇=2.(18)𝑠env is used to switch environmental effects off (0) or on (1). Air temperature 𝑇air above a threshold 𝑇air0 increases virus mortality according to power 𝑞𝑇. Parameter 𝑞𝑇 is assumed equal to two giving a quadratic dependence of temperature above the threshold 𝑇air0 (Table 1). It is usual for the biological effects of temperature to be nonlinear, sometimes approximating to exponential, as in the use of a 𝑄10 factor for the consequences of a 10°C temperature rise on chemical reaction rate, or the application of the Arrhenius equation for chemical reactions (e.g., [44], pp. 103–105). 𝑘𝑇air is a rate parameter. We did not find controlled-environment studies on the effect of temperature on virus longevity which we could use, and therefore the values assigned the parameters are estimates.

In southern Britain, daily mean air temperature 𝑇air varies from c. 3 (January) to 17°C (July) ([44], p. 270; [49] p. 142). It is assumed that 𝑇air varies sinusoidally, with𝑇air=𝑇air,mean+𝑇air,varsin2𝜋𝑖365Julian𝜏Tair;𝑇air,mean=10C(varied),𝑇air,var=7𝜏C,𝑇air=115day.(19) Annual mean and seasonal variation are 𝑇air,mean and 𝑇air,vr. 𝑖Julian is the Julian day number (1 on 1 January, leap years are ignored). 𝜏Tair (d) is the phase of the sinusoid, which is maximum on 25 July. Combining equations (19) with (18) modifies environmentally induced virus death rate, 𝑘vir,env, and hence reproductive ratio (20). The same seasonal pattern is applied every year. These equations for temperature are a good approximation to long-term weather means in the UK (Note: the study reported in [34] suggests that, on an annual timescale, cold weather does not predict winter deaths, but, on a shorter term timescale, cold weather could be a significant trigger).

2.2. Basic Reproductive Ratio, R0, and the Disease Generation Time, 𝜏gen (Day)

Two important epidemiological parameters are basic reproductive ratio, 𝑅0, and the disease generation time, 𝜏gen (day). These are both derived from the basic parameters of the model. 𝑅0 is defined as the number of infections directly caused by a single infected individual during its infectious period when in a population of susceptibles. 𝜏gen is the average time it takes the direct infections which contribute to 𝑅0 to arise. 𝜏gen is also called the “serial interval”, the average time between the primary case and secondary cases.

In the absence of forcing (see previous section, in which case 𝑅0 can only be calculated numerically), 𝑅0 can be calculated analytically [6] by travelling round the infectious loops in Figure 1 and adding the terms together. This leads to (see Appendix B for an alternative equivalent statement of 𝑅0) 𝑅0=𝛽𝑘vir,nm+𝑘vir,env𝑘vir𝑘inf𝑘inf×𝑤+𝜇inf1𝑘inf+1𝑓asy×𝑘inf𝑘inf𝑤+𝜇inf2𝑘inf+𝑘inf𝑘inf+𝜇2𝑤inf3𝑘inf+𝑘inf𝑘inf+𝜇3×𝑤inf4𝑘inf+𝑘sic𝑘sic+𝜇+𝜇sic×𝑤sic1𝑘sic+1𝑓hos𝑘sic𝑘sic+𝜇+𝜇sic𝑤sic2𝑘sic+𝑘sic𝑘sic+𝜇+𝜇sic2𝑤sic3𝑘sic+𝑘sic𝑘sic+𝜇+𝜇sic3𝑤sic4𝑘sic+𝑓hos𝑘hos𝑘hos+𝜇+𝜇hos𝑤hos2𝑘hos+𝑘hos𝑘hos+𝜇+𝜇hos2𝑤hos3𝑘hos+𝑘hos𝑘hos+𝜇+𝜇hos3𝑤hos4𝑘hos+𝑓asy𝑘asy𝑘asy𝑤+𝜇asy2𝑘asy+𝑘asy𝑘asy𝑤+𝜇asy3𝑘asy+𝑘asy𝑘asy+𝜇2𝑤asy4𝑘++asy𝑘asy+𝜇6𝑤asy8𝑘asy.(20)If (20) is multiplied out, each term corresponds to one of the 18 weighting factors in (14) (4(inf) + 4(sic) + 3(hos) + 7(asy) = 18) and represents one complete infective loop passing through the virus pool in Figure 1 (see Appendix B, (B.1)– (B.4)). The first term is the simplest loop, passing from 𝑛vir to 𝑛inf1 and back to 𝑛vir (Figure 1). The first term is𝛽𝑘vir,nm+𝑘vir,env𝑤inf1𝑘vir𝑘inf𝑘inf𝑘inf+𝜇.(21) Start with a single virus particle in the 𝑛vir box (we could equally well start with a single person in the inf1 first-infected compartment). The first factor of (21), with units of persons per virus particle, is the probable number of persons infected by a virus particle during its life; the mean lifetime of a virus particle is 1/(𝑘vir,dnm+𝑘vir,env). It is assumed that there is no additional virus death process due to exposure. 𝑤inf1 is a dimensionless weighting factor (14). The third factor (𝑘vir/𝑘inf) is the number of virus particles produced per person in the inf1 state during his life: 1/𝑘inf is the average lifetime of a person in the inf1 state. The last factor, with 𝑘inf+𝜇 in the denominator, is the probability that this lifetime (of 1/𝑘inf) is actually achieved. We can continue in this way via all the compartments which can give rise to infection (i.e., produce virus particles (13)), adding up the terms. (20) for 𝑅0 can be written out as a sum of the 18 (potentially) contributing terms (Appendix B)—an equivalent formulation which is sometimes useful.

Equation (20) can be (and was) checked numerically by placing a single infected individual into the 𝑛inf1 box and diverting the (primary) infected individuals into an accumulator, rather than allowing them to enter the 𝑛inf1 box, where they can lead to secondary infections. The algebraic and numerical methods agree to six decimal digits, giving a basic reproductive ratio 𝑅0=4.78 for the default parameters without forcing.

A “dynamic” reproductive ratio, 𝑅dyn, can be calculated when the total population is not entirely susceptible by means of𝑅dyn=𝑅0𝑛sus𝑛liv.(22) The live population, 𝑛liv, is given by (1). This allows approximately for a slowly changing fraction of susceptible individuals. However, 𝑅0 may itself be dynamic on a shorter time scale due to seasonal effects on virus death rate 𝑘vir,env (see previous section).

The disease generation time, 𝜏gen (day), can be computed analytically by (B.7), giving 𝜏gen=2.51 day. A numerical calculation using Runge-Kutta integration gives also 𝜏gen=2.51 day agreeing with the analytical result to nine decimal digits, although Euler integration gives 𝜏gen=2.48 day (in each case the integration interval is 1/32 day; see first paragraph of Section 3). Due to overlapping generations, 𝜏gen is not equal to the time constant at which the total infecteds increase. Assume that total infecteds 𝑛iash (1) increase initially at an exponential rate according to𝑛iashe𝑡/𝜏,(23) where 𝜏 (day) is the growth time constant (the proportional growth rate of infecteds is 1/𝜏 day−1). This depends on both 𝑅0 and 𝜏gen and the structure of the model [6, Section  1.2.3]. This can be extracted numerically (there is a period of constant exponential growth that lasts for some 10 days) to give τ = 1.29 day. This is considerably less than the generation time of 2.51 day.

2.3. Parameterization

Parameters have been introduced while developing the model. Their values are listed for reference in Table 1. Here we summarize the evidential basis for the parameter values used.

As a preliminary, statements from two physicians are quoted. In [50] a general practitioner in the Doncaster (UK) area describes his study of the 1969-1970 pandemic as it affected his urban practice. He said “the true incidence of influenza during an epidemic is probably impossible to assess.” This is perhaps equally true today and sets the scene for the significance of processes such as “validation” and data fitting (see Figure 5 and Section 3.4). Another general practitioner, this time in Kent (UK) [51], states that “an epidemic of influenza tends to last in this area for between two and three months. Beginning slowly the epidemic reaches its peak in four to five weeks and then subsides slowly. The extent and severity of any attack will depend on factors such as the strain of influenza virus, on the state of the host-immunity of the population and on the timing of the epidemic; the fatality and complication rates are always higher during the cold and foggy winter months.” His statement agrees with many of our seasonal-forcing simulations (Figures 6 and 7).

Five key epidemiological quantities for influenza are (a) the time which elapses after the infection event until the subject becomes infectious to others, denoted by the 𝜏inf0 (day); (b) the latent period, namely the time which elapses between the infection event and the appearance of clinical sickness, denoted by 𝜏lat (day); (c) the infectious period 𝜏inf (day), which is the time during which the subject is infectious (Figure 1—producing virus); (d) the period of immunity 𝜏imm (day)—the time period after recovery from influenza during which the subject has immunity, before gradually losing it and returning to the susceptible pool (Figure 1, 𝜏imm is represented by the box at the bottom of the diagram with 10 sequential pools); (e) the basic reproductive ratio 𝑅0 (dimensionless) (20).

Addressing these quantities, first consider the time period between the infection event and becoming infectious, 𝜏inf0 (day). With the infectivity weighting factors in (14), the first infected pool 𝑛inf1 with a mean lifetime of 0.5 day is not infectious but the second infected pool 𝑛inf2 is, and therefore𝜏inf0=1𝑘inf=2day1=0.5day.(24) Next the latent period 𝜏lat (day) is given by (there are four sequential pools in the infected category of Figure 1)𝜏lat=4𝑘inf=2day1=2day.(25) The infectious period 𝜏inf (day) can be estimated as follows. Note that (a) all paths from the susceptible category to the recovered category in Figure 1 pass through eight pools with the same outgoing rate constant (natural mortality excluded); (b) it is assumed that the disease-related rate constants 𝑘inf,𝑘asy,𝑘sic, and 𝑘hos are all equal (Table 1), say to 𝑘dis; (c) the first infected pool (𝑛inf1) is assumed not infectious and all sick and hospitalized pools are infectious (14). We therefore write the infectious period 𝜏inf as𝜏inf=7𝑘dis=2day1=3.5day.(26) Loss of immunity occurs during a time period of 𝜏imm, assumed to be 3 years. This delay is represented by (Figure 1) 10 sequential pools each having an outgoing rate constant of 𝑘del (day−1) (ignoring natural mortality (12)). This gives a gamma-distributed delay (Appendix A, (A.5) to (A.8) with 𝑚=10). 𝜏imm and 𝑘del are related by𝑘del=10𝜏immwith𝜏imm=1095day=3year.(27) Next compare the values in (24)–(27) with the literature. Our value of 𝜏inf0=0.5 day (24) can be compared with that of [25], who referring to [45], use a rather different value of 1.9 day. However, the value given in [45] is based on fitting a homogenous-mixing deterministic SEIR (susceptible, exposed (meaning infected but latent), infected (meaning clinical), resistant) model to the excess pneumonia and influenza deaths in 45 cities during the 1918 pandemic. [23] use a value of 0.5 day in their model, without citing a specific source. Our value could be doubled to 1 day by taking 𝑤inf2=𝑤asy2=0 in (14), a relatively small change to the model. Both 0.5 and 1 day are compatible with the clinical evidence, although not with the 1.9 day value of [45].

Our latent period of 𝜏lat = 2 day (25) is both clinically acceptable and is not very different from [25]'s (fitted) values of 2.4 and 3.6 days for the spring and autumn waves of the 1918 epidemic. Note [25] uses the term “latent” to describe the period between the infection event and the time at which the person becomes infectious, rather than our use (25) which refers to the period between the infection event and the time at which the person becomes clinically sick.

[41] report that the “infectious period” 𝜏inf lies in the range of 6 to 10 days, but then their single infected pool represents everything in our Figure 1 between the susceptible and recovered categories which makes meaningful comparison difficult. [41] states that “for simplicity, we do not explicitly model the exposed population but instead include people infected but not yet infectious in the “I” box. Including an exposed class yields similar results.” The last sentence is not supported by simulations. Students of dynamical systems may find this statement surprising. However, their range of 6 to 10 days is very different from our 3.5 day (26) which is compatible with the clinical evidence and also the 2 or 3 day period obtained by [25] when fitting the spring and autumn waves of the 1918 pandemic.

The loss of immunity of recovered patients is not a feature of the model in [25]. [41] uses a range of 4 to 8 years: we were not able to chase their values to earth. We employ 3 years, as in (27) for 𝜏imm, which, while clinically reasonable, should be regarded as a guess at what is probably a rather variable quantity.

The basic reproductive ratio 𝑅0 (20) is hardly a clinical or easily observed quantity, but it is much loved by modellers and epidemiologists, not without reason. It is an “emergent” parameter of the model, as its value depends on underlying parameters. Because influenza is highly seasonal, it can be concluded that the environment is important, yet values of 𝑅0 given in the literature rarely say anything about the environment. [41] gives its range as 4 to 16. [45] says that “estimates vary widely, varying from 1.68 to 20.” After fitting the Geneva data for the 1918 pandemic, [25] gives values of 1.49 and 3.75 for the spring and fall waves. [27] states “Assuming a basic reproduction number of 𝑅0 = 2.5 and using the standard parameter set of InfluSim” [26]. We assume that this value of 2.5 is the outcome of their standard parameter set. Other values given are 1.2 to 2.4 [7, Table 2] and 2.07 ([23], after their Figure 3). Since many of these models use biologically inappropriate assumptions [43], it is relevant to quote the statement from [43] that “ignoring the latent period or assuming exponential distributions will lead to an underestimate of 𝑅0 and therefore will underestimate the level of global control measures … that will be needed ….” [52] reviews 𝑅0 values of several flu epidemics and pandemics, focussing on the possible control of an H1N1 epidemic. Our standard parameter set without any environmental forcing gives a value of 𝑅0 = 4.78 (20). This can be regarded as an upper baseline applying to an optimum environment (one maximizing 𝑅0), and changing the environment can only decrease this number (Figure 6).

The disease generation time used in [23] is 2.44 day, which is close to our value of 2.51 given following (23).

Finally, we wish to comment on the number of pools used to represent the infected and sick categories in Figure 1, where there are four in each. Using trajectory matching on an influenza outbreak at an English boarding school [43, 53] suggests that two pools are appropriate for each category. Unfortunately they give few details of their procedure. We found a good fit (Figure 5) using the current model to fit the same data. We therefore continued to use four pools, although arguably, whether or not two, or three or four pools, or even a nonintegral number of pools are applied is perhaps less important than the principle of applying two or more pools.

3. Results

This section reports simulations of the model described above. First we describe the general technical aspects applied in the simulations.

3.1. Numerical Methods

The model was programmed in ACSL (Advanced Continuous Simulation Language, Aegis Research, Huntsville, AL, USA; version 11.2.2 for DOS), an ordinary differential equation solver. In all simulations, equations are integrated using Euler's method, a fixed integration interval of Δt = 0.03125 = 1/32 day (45 minutes), and results were communicated for plotting at half-day or daily intervals. There were no difficulties with model implementation. Not unexpectedly, in some of the chaotic simulations, different (but still chaotic) results were obtained if different integration methods and intervals or different Fortran compilers were used (the results shown used the Watcom compiler).

The simulations focus on the daily hospital admission rate, 𝐼hos2 (persons day−1); Figure 1), as this is often a recorded statistic. An “epidemic” is deemed to have occurred if 𝐼hos2 shows a maximum (in steps of Δt) and if at the maximum 𝐼hos20.1 persons day1.

State variable initial values and parameter values are as listed in Table 1 (unless stated otherwise). In general, the parameters have not been tuned for any particular performance (but see Figure 5) and as far as possible have been estimated mechanistically (see Section 2.3). Some of the simulations (e.g., Figures 2, 3, and 4) are intended to illustrate important characteristics of the model—they are not intended to be compared with actual epidemics/pandemics. Other simulations (e.g., Figure 5, parts of Figures 7 and 8) are intended to demonstrate at least a partial realism and to provide credibility to the model.

fig2
Figure 2: Dynamics without loss of individual immunity or seasonal forcing. There is no loss of individual immunity (𝑘rec=0; (10), Figure 1) nor seasonal forcing ((18); 𝑠env=0). (a) Short-term dynamics. (b) Long-term dynamics. Parameters and initial conditions have the values in Table 1. Successive epidemics are caused by loss of immunity at the population level from the natural birth and death processes (μ, (2), Table 1).
fig3
Figure 3: Effect of varying key parameters. Parameters and initial conditions have the values in Table 1. There is no loss of individual immunity (𝑘rec=0; (10), Figure 1) and no seasonal forcing ((18); 𝑠env=0). In each case, the value of the altered parameter, basic reproductive ratio 𝑅0 (20), and total hospital admissions (HA, the integral of 𝐼hos2, (5)) are given; the curve for the default parameters is thickened and indicated. (a) Changes to infection parameter β. (b) Changes to 𝑘inf(=𝑘sic=𝑘hos=𝑘asy; Figure 1). (c) Changes to the initial (time t = 0) immune fraction in the recovered category, 𝑓rec0. This is achieved by altering the time t = 0 value of 𝑛rec state variable (Figure 1) to 𝑛rec(𝑡=0)=𝑓rec0𝑛city (Table 1, (22)).
fig4
Figure 4: Dynamics with loss of individual immunity. Individuals now lose immunity (𝑘rec=1; Figure 1). Parameters and initial conditions have the values in Table 1. (a) Short-term dynamics: susceptible 𝑛sus and delayed 𝑛del numbers are calculated with (2) and (1), the dynamic reproductive ratio 𝑅dyn with (22), and hospital admission rate 𝐼hos2 with (5). (b) Longer-term dynamics.
471653.fig.005
Figure 5: Model validation. The model is validated by comparing predictions (solid line) with data (●) from [53]. See text for details.
fig6
Figure 6: Environmental effects on basic reproductive ratio, 𝑅0. (a) Seasonal change in virus inactivation rate, 𝑘vir,env (Figure 1; (15), (18)), giving virus inactivation by temperature alone. 𝑇air is mean daily air temperature (appropriate to southern Britain) (19); 𝑇air0 is an assumed temperature threshold for virus inactivation, here equal to 13°C ((18), Table 1)); 𝑅𝟎 is the seasonally modulated basic reproductive ratio (20), calculated assuming temperature is constant. (b) For different values of mean annual air temperature, 𝑇air,mean (19), basic reproductive ratio 𝑅0 (20) is seasonally modulated as in (a) via the seasonal modulation of 𝑘vir,env, for a constant threshold temperature for virus inactivation of 𝑇air0 = 13°C (18), giving for 𝑅0 a maximum, minimum, annual mean, and amplitude (maximum-minimum). 𝑁fdy denotes the number of forcing days per year (a “forcing day” is a day on which air temperature is sufficiently high as to increase virus mortality).
fig7
Figure 7: Weak environmental forcing. Seasonal forcing (18) combined with loss of individual immunity (Figure 4; Figure 1, 𝑘rec=1 day−1). The system is initially in a steady state without forcing (Figure 4(b)); forcing is applied after 1 year. The effects of seasonal forcing on the basic reproductive rate 𝑅0 (20) are shown in parts (a)–(c) (right-hand 𝑦-axis). Increases in mean annual air temperature 𝑇air,mean increase the level of forcing ((18), (19), Figure 6(a)). 𝑅0 is the mean annual value of 𝑅0. Parts (c) and (d) are from a single simulation and have the same level of forcing.
fig8
Figure 8: Moderate environmental forcing. Seasonal forcing (18) combined with loss of individual immunity (Figure 4; Figure 1, 𝑘rec=1 day−1). The system is initially in a steady state without forcing (Figure 4(b)); forcing is applied after 1 year. The effects of seasonal forcing on the basic reproductive rate 𝑅0 (20) are shown (right-hand 𝑦-axis). Increases in mean annual air temperature 𝑇air,mean increase the level of forcing ((18), (19), Figure 6(a)). 𝑅0 is the mean annual value of 𝑅0.
471653.fig.009
Figure 9: Gamma function basics. Illustration of consequences of using several sequential pools in representing a clinical category. Results for a single pool (a), two pools (𝑏2), three pools (𝑐3), four pools (𝑑4), and eight pools (8) are given, showing in each case the value of the state variable for the final pool in the sequence. For 𝑚 pools, the rate constant 𝑘m has been set to give the same mean overall transit time of 𝑚/𝑘m=4 day (indicated by ■). For each curve, the maximum value (●) denotes the most probable overall transit time.
3.2. Dynamics without Intrinsic Loss of Immunity and without Forcing

Here the model is exercised with parameter 𝑘rec=0 ((10), Figure 1, Table 1), so that there is no loss of immunity as represented by the delay box in Figure 1, and without environmental forcing (18). In this case, the natural death and birth rates (μ of Table 1 and (2) to (12)) lead to a slow loss of immunity at the population level as births are assumed to be susceptible.

3.2.1. Short-Term Dynamics

These are illustrated in Figure 2(a) and are unexceptionable. From initial infection, it is two to three weeks before the disease is visible, and the epidemic is over within a further two weeks. Infected number (𝑛inf) peaks a few days before hospitalized numbers (𝑛hos), both then falling to very low values. Most (99%) of the population joins the recovered (immune) box (𝑛rec), with 1% escaping infection altogether. 90% of those infected travel via the asymptomatic route (Figure 1; 𝑓asy, Table 1, (7)). The epidemic in Figure 2(a) is short compared to UK experience, but Figure 2(a) is for the no-seasonal-forcing situation where initially the population is 100% susceptible, which is not comparable with actual epidemics where seasonality is always a factor as is partial immunity. Later (Figure 5), it is shown how the model, without significant “tuning”, is able to fit data on a UK epidemic. Also, seasonal forcing lengthens the period of the epidemic and decreases the fraction of the susceptible population which becomes infected (Figures 7 and 8).

3.2.2. Long-Term Dynamics

These are illustrated over 100 years in Figure 2(b). On this time scale, the first epidemic, shown in Figure 2(a), occurs at zero time and rapidly dies out with the dynamic reproductive ratio 𝑅dyn (22) (not plotted) falling to near zero, as the fraction of susceptibles (𝑛sus/𝑛liv) becomes small. In the first epidemic, at time 𝑡20 day, the initial peak in hospital admissions 𝐼hos2=949 (Figure 1, (5)), decreasing to less than 100 at the second epidemic. After the first epidemic, 𝑅dyn (22) slowly recovers as immune numbers (𝑛rec) decrease and susceptibles (𝑛sus) increase due to the natural birth and death processes. The second smaller epidemic occurs after a further 28 years, followed by epidemics of decreasing amplitude and increasing frequency until a steady state is reached with fractions of susceptibles of 0.209, of all four infected categories together (infected, asymptomatic, sick, hospitalized, Figure 1) of 0.00014 and of recovered (immune) of 0.791; there are 0.18 hospital admissions per day, and 𝑅dyn (22) = 1.

3.2.3. Responses to Key Parameters

Figure 3 illustrates the effects of changing three key parameters: infectivity parameter β (Figure 1; (2); Table 1), rate parameter 𝑘inf (with 𝑘inf=𝑘sic=𝑘hos=𝑘asy; Figure 1; Table 1; (3)), and the initial (time zero) immune fraction, 𝑓rec0. Figure 3(a) shows how the duration, initial proportional growth rate, and severity (indicated by hospital admissions, HA) of a modelled epidemic are influenced by the value of infectivity parameter β (or equally 𝑘vir; both are linear factors of 𝑅0 (20)). Indeed, the effects of increasing β are monotonic, moving the epidemic towards shorter time scales, increasing initial proportional growth rate, total hospital admissions (HA) towards an asymptote of 5000, and narrowing the width of the epidemic. Also, increasing β causes the long-term steady state to be more quickly attained—the spikes of 𝐼hos2 (e.g., Figure 2(b)) becomes closer together.

Figure 3(b), where the rates of transit of infected persons through the system are increased (𝑘inf=𝑘sic=𝑘hos=𝑘asy=𝑘xxx, Figure 1), is less straightforward. Here a low value for the 𝑘xxx (e.g., 0.25 day−1) gives a high basic reproductive ratio 𝑅0 (20) (Appendix B) an epidemic which is slow to take hold but (surprisingly) moves more rapidly towards a long-term steady state (i.e., the 𝐼hos2 spikes as in Figure 2(b) are closer together). As the 𝑘xxx are increased basic reproductive ratio 𝑅0 decreases (20) first the epidemic becomes narrower and faster (e.g., 𝑘xxx=2 day−1), but further increases in 𝑘xxx (e.g., 𝑘xxx=5 and 8 day−1) cause the epidemic to become slower to take hold and less peaked, with fewer hospital admissions (HA). Increasing 𝑘xxx always causes the 𝐼hos2 spikes (Figure 2(b)) to move further apart and lengthens the time taken to reach a steady state and decreases the number of infected persons in the steady state (𝑛iash,(1)).

Finally, Figure 3(c) illustrates how increasing the initial immune fraction 𝑓reco has a similar effect to that of decreasing β: delaying the onset of the epidemic and increasing its width, decreasing the initial proportional growth rate, and decreasing severity (number of hospital admission, HA).

3.3. Dynamics with Intrinsic Loss of Immunity and without Forcing

Now the discrete delay box of Figure 1 giving loss of immunity is switched on by making𝑘rec=1 day−1, which causes recovered individuals to be moved quite rapidly into the delay sequence of ten compartments where immunity is lost after three years (𝜏imm=3×365days; Table 1, (11)). Immunity is also being lost due to the natural birth and death processes, as newborns are assumed to be susceptible. The responses of Figure 2(a) are little changed by taking 𝑘rec=1 day−1 (instead of 0) if the 𝑛rec variable of Figure 2(a) is replaced by the 𝑛del variable (1). However, over a longer time period Figure 4(a) illustrates that there is now a switching of susceptibles between a low (fractional) value and a high (fractional) value (fractions of 𝑛city, Table 1, (1) and following paragraph), with an inverse switching of numbers in the delay compartment, 𝑛del (Figure 1). Figure 4(b) shows that the oscillations, as in Figure 2(b), decrease in amplitude and increase in frequency until a steady state is attained. In the steady state, there are 3.7 hospital admissions per day, a dynamic reproductive ratio 𝑅dyn = 1 (22), and the fractions in the three principal categories (Figure 1) of susceptible; all infected categories lumped together (infected, asymptomatic, sick, hospitalized, 𝑛iash (1)) and the delayed categories are 0.209, 0.0030, and 0.787. Comparing these values with the situation in which there is no intrinsic loss of immunity (Figure 2(b)), it is seen that hospitalizations have increased by a factor of 20 as has the fraction in all infected categories, although the fractions of susceptibles and recovered or immunes have barely changed.

3.4. Validation

“Validation” is often a misused and misunderstood concept and is perhaps better described as an evaluation of applicability. Validity is not a property of a model alone, neither is it a “zero or one” concept. It describes the relationship between model predictions and a set of data obtained under prescribed conditions. In this section, the model of Figure 1 is “validated” by fitting the predictions of the model, with minimal parameter adjustment (a perfectly formulated mechanistic model would permit no adjustment of parameters or initial values) to data on an influenza epidemic [53]. Success in this endeavour gives the model some credibility, although it does not make the model valid for general use.

The data in [53] relate to an influenza outbreak at an English boarding school, which provides a simple situation that seems comparable to the single-city homogeneous-mixing model of Figure 1, remembering that the model does not have the stratification which might be needed if a larger region was considered.

Minimal tuning is applied. One infected person is introduced to the school at noon on 18 January, otherwise all are susceptible. Table 1 parameters are altered for the Figure 5 simulation as follows: infection parameter β is changed to 0.07 (virus units)−1 day−1; all birth and death rates are set to zero (𝜇,𝜇hos,𝜇sic) for such young persons; total population 𝑛city=763; and the asymptomatic fraction, 𝑓asy, is one-third, to reflect the finding that only two-thirds of the boys became sick, and the community is assumed to be “well-mixed”. There is no initial immunity. With these values, the basic reproductive ratio 𝑅0 (20) is 6.44, the mean generation time is 2.65 day (A.7), and the initial proportional growth rate of total infecteds, 𝑛iash(1), is 1.015 day−1.

For comparison with the data in [53], we define the number of persons confined to bed, 𝑛ctb, as𝑛ctb(𝑡)=𝑛sic1(𝑡)+𝑛sic2(𝑡)+nsic3(𝑡)+𝑛sic4(𝑡)+𝑛hos2(𝑡)+𝑛hos3(𝑡)+𝑛hos4(𝑡).(28) Since this definition is somewhat arbitrary, in comparing the predicted 𝑛ctb from the model with the data from Anon (1978), 𝑛ctb(𝑡) is scaled with an adjustable factor so that the comparison line drawn in Figure 5 is1.35𝑛ctb(𝑡).(29) Fitting was done by eye, as this can produce (see below) a better focus on the biological significance of the parameter being adjusted and possible limitations in the biological data that may be obtained with more automated methods. See [29] for possible problems arising from formulaic parameter adjustment in mechanistic models.

The degree of fit in Figure 5 is satisfactory. Apart from the two outliers on 26th and 27th January, the fit is good. In an actual epidemic, there may be underreporting during the early states and overreporting later, as the performance of those handling the epidemic changes. Overall, we believe it is reasonable to assume that the model has some credibility as a result of this validation.

Finally, a comment on whether the number of pools used in the infected and sick categories in Figure 1 is appropriate. [43] fits an SEIR model to the observed data given in [53] and shown here in Figure 5. They minimize the sum of the squared errors, arguably this underweights the skirts of the distribution, which are sensitive to the numbers of sequential pools assumed (Figure 9). They assume 𝑚 pools for the E category and 𝑛 pools for the I category. They find a best fit with 𝑚=2 and 𝑛=2. We were not able to discover the details of their general parameterization or indeed how they define a confined-to-bed number from the categories and pools of the SEIR model. They remark that there is a sensitivity to the number of points used to obtain the fit and that basic reproductive ratio 𝑅0 can change substantially. More notably, it can be seen in their Figure 3(c) that the model fits the last three data points as the epidemic is subsiding rather poorly. This suggests (see Figure 9) that a higher number of pools is required than their best values of 2 for the exposed and infected categories. In view of these difficulties, and the comparison shown in Figure 5, we consider that the number of pools per category suggested in Figure 1 and used throughout this paper is reasonable. Although our particular choice cannot rigorously be defended, it seems to be “good enough” at the present time.

3.5. Dynamics with Intrinsic Loss of Immunity and Seasonal Forcing

Now we add direct seasonal forcing by weather to the simulations illustrated in Figure 4 ((18), (19)). Four weather factors which could impinge on virus longevity are air temperature (𝑇air), relative humidity (𝑅H), radiation (possibly multicomponent), and wind speed. The effects can be highly complex: for example, [48] which examines the effects of 𝑇air and 𝑅H on virus transmission, finding that cold dry conditions favour transmission (but see their Figure 6). The topic is far from being well understood, but since our concern here is to represent broadly weather forcing within the model, we make the simplifying assumption that 𝑇air alone is operative.

Figure 6(a) illustrates the seasonal variation of mean daily air temperature 𝑇air in the southern UK. Using (19) with (18), and 𝑠env=1, this affects environmentally induced virus death rate 𝑘vir,env (Figure 1, (15)), and thereby basic reproductive ratio 𝑅0 (20). Mean daily air temperature 𝑇air varies between 3°C (Jan 24) and 17°C (Jul 25). When 𝑇air crosses the temperature threshold of 𝑇air0 = 13°C in May, this increases environmentally induced virus mortality, 𝑘vir,env (18) and decreases basic reproductive ratio 𝑅0 (20). With this formulation, influenza is most likely during the months of October to April when 𝑅0 is highest.

Figure 6(b) shows how varying (e.g., increasing) mean annual air temperature, 𝑇air,mean (19), (or equivalently, decreasing threshold temperature 𝑇air0, (18)) varies the duration and intensity of seasonal forcing of 𝑅0. Changing mean annual air temperature, 𝑇air,mean (19), changes the average annual value of 𝑅0, 𝑅0, as well as its maximum, minimum and amplitude (maximum-minimum). The dependence of these quantities on 𝑇ait,mean is also shown, together with the number of forcing days per year (𝑁fdy, a forcing day is one on which mean air temperature is above threshold temperature and virus mortality is reduced). With mean annual air temperature 𝑇air,mean6C, there is no forcing at all (𝑁fdy=0): air temperature 𝑇air (19) never exceeds threshold temperature 𝑇air0 = 13°C (18); 𝑅0 is invariant at its maximum value. Environmentally induced virus death rate 𝑘vir,env ((15), (18)) stays at zero; the system remains in the steady state shown in Figure 4(b) and 𝑅0=𝑅0 (this situation is equivalent to 𝑇air,mean=10C and 𝑇air0 = 17°C). As mean annual air temperature 𝑇air,mean increases above 6°C, there are more days in the summer months when mean daily temperature 𝑇air> threshold temperature 𝑇air0, number of forcing days per year 𝑁fdy increases, and 𝑅0<𝑅0. The amplitude increases to a maximum when𝑇air,mean=21C, before decreasing. With 𝑇air,mean>28.5C, influenza epidemics do not occur, because the annual average of 𝑅0 is less than unity.

In the simulations presented in Figures 7 and 8, various values of mean annual air temperature 𝑇air,mean are taken between −4.5 and 16°C, assuming always an annual variation 𝑇air,var of ±7°C (19) as in the UK. Otherwise, parameters have the values in Table 1. In these simulations, the long-term steady state reached in Figure 4, with loss of immunity (𝑘rec=1; Figure 1; (10)) but without seasonal forcing, is used for initial values. Forcing is applied after one calendar year. The aim is to illustrate the wide variety in dynamic behaviour which results from this type of seasonal forcing (which decreases the reproductive ratio). In each case, the mean value, its maximum, minimum and amplitude (maximum-minimum) of basic reproductive ratio can be read off Figure 6(b). In the UK the influenza season is considered to be over by May, when the mean daily temperature ranges from 10.7 (1 May) to 14.1°C (31 May) (Figure 6). Therefore, of the simulations described in Figures 7 and 8, those given in Figures 7(b), 7(c),7(d), 8(a), 8(b), and 8(c) are more relevant to the UK.

Figure 7 illustrates the effects of low levels of seasonal forcing on influenza hospital admissions, 𝐼hos2 (Figure 1,(5)). Forcing is seen as a downward modulation of the basic reproductive ratio 𝑅0, which decreases the annual average, 𝑅0.

The lower graph in Figure 7(a) shows that 𝑅0 is slightly decreased beginning on 24 June, from 4.78 in the steady state to 4.69 on 25 July. Influenza hospital admissions (the upper graph in Figure 7(a)) oscillate twice about the steady state value (Figure 3(b); 3.7 admissions day1) in a 12-month period with the two steady-state maxima 168 days apart on 27 September and 14 March. A regular variation is quickly established. Note that the natural response time of the system as indicated by the upper graph in Figure 7(a) is not commensurate with the annual cycle imposed by the environment (shown by the lower graph in Figure 7(a)). This sets the scene for potential chaos.

In Figure 7(b) the level of forcing is increased (lower graph). This results in a more complex (but still regular) schedule of hospital admissions with a biennial pattern superposed on a twice-yearly variation. In [38] monthly data are presented in their Figures 2 and 3; some of these are suggestive of a twice-yearly pattern.

A further increase in forcing (Figures 7(c) and 7(d)) results in chaos, with sometimes one, two or even three peaks in hospital admissions occurring within a twelve-month period. However, there is a tendency towards annual epidemics, with 387 epidemics occurring in 250 years, and the epidemics (10% points) lasting about seven weeks (cf. [51] which gives a duration of two to three months; also cf. Figure 8(a) with annual late spring epidemics of 5-week duration). The susceptible fraction of the population varies from c. 30% before an epidemic to 15% just after each epidemic, so that one half of the susceptibles becomes infected.

Further increases in forcing are shown in Figure 8. First, in Figure 8(a) with mean annual air temperature 𝑇air,mean=11C (Figure 6), there is a transition to an annual epidemic occurring in the late spring of each year lasting for about five weeks. In these annual epidemics, the susceptible fraction falls from c. 35% to 12%. Next (Figure 8(b); 𝑇air,mean=12C), chaos is again produced (cf. Figure 7(d)) with a strong tendency towards annual epidemics: 204 epidemics occur in 200 years. In the last two years of the simulation shown in Figure 8(b), there is a small epidemic on 22 July, followed by larger epidemics on 5 October and 16 May. Then (Figure 8(c); 𝑇air,mean=13C), there is chaos but now with a tendency towards biennial epidemics (169 epidemics occur during 200 years). Last (Figure 8(d); 𝑇air,mean=15C), the system immediately settles down into regular biennial epidemics occurring in early spring of every other year but the amplitudes remain slightly chaotic.

Note that, throughout Figures 7 and 8, as mean annual air temperature threshold 𝑇air,mean is increased, forcing is increased (i.e., the magnitude of the seasonal changes in basic reproductive ratio 𝑅0 ((20), (18), (19), Figure 6) increases), but the mean annual value of 𝑅0, 𝑅0, decreases. This causes between-epidemic recovery time to increase (see discussion of Figure 3(a) above). The frequency of epidemics then decreases.

Increases in mean annual air temperature 𝑇air,mean can be continued, and although the situations simulated are now less realistic, they do contribute to understanding the system. With 𝑇air,mean=19C, the response is chaotic with 37 epidemics in the first 200 years. Mostly five or six years elapse between epidemics; occasionally two smaller epidemics occur in the same year (e.g., in the 179th year). The mean reproductive ratio 𝑅0 is 2.945. When 𝑇air,mean=21C, the response eventually settles down to a regular pattern with 27 epidemics in the first 200 years and eight years between epidemics and an 𝑅0 of 2.499. When 𝑇air,mean=26C the response is chaotic with 15 epidemics in the first 200 years and 𝑅0=1.367. Last, with 𝑇air,mean=27.5C, after departing from the initial steady state, it is c. 110 years before influenza reemerges, and then it is at a low level, eventually settling down to a repeating annual late spring epidemic lasting about two months with a maximum of 0.6 hospital admissions per day and an annual sum of 20 hospital admissions per year. 𝑅0 has a modest mean annual value (1.131) and is below unity for much of the year.

4. Discussion

Many models of influenza are more empirical than mechanistic, and therefore, although they can be and sometimes have been used to fit historic data [54], they are of little value in further understanding or for indicating how future epidemics/pandemics might be handled before they occur ([32, 42, 55, 56]).

Seasonality is an important feature in influenza incidence. There are many ways in which seasonality can be incorporated into an influenza model. In [42] the contact rate is reduced by a factor of ten for the 6-month nonepidemic season. Here a simple representation of UK daily weather allows the impact of mean daily air temperature to be explored. A similar approach could be applied to relative or absolute humidity, radiation, and wind speed. Seasonal forcing gives rise to wide range of dynamics, from regular at various intervals, to chaos, as illustrated in Figures 7 and 8.

Some of simulations of Figures 7 and 8 are similar to influenza incidents which have occurred. For example, in Figure 7(d), the three maxima near year 184 in the spring, autumn and following spring have similarities with the waves of the 1918–1920 influenza pandemic [57]. Note that this is achieved within a single simulation without changing parameterization. In comparison, in [25] the authors fitted these data with a model, applying the model separately to each wave, with different parameters and initial values (loc. cit. Table 1, Figure 4). It is legitimate to ask just what this procedure means. The two principal peaks occurring within a single influenza season in the autumn and spring around year 199 of Figure 8(b) resemble the peaks shown in the 1957-58 pandemic [20, Figure 3A]. The first two peaks illustrated in Figure 8(c) in year 182 occur in late spring (the lesser peak) and the following autumn (the main peak), resembling the 1968–70 pandemic [50].

Apart from Section 3.4 and Figure 5 which applies to an unusually sharply defined context, we did not attempt to fit our model to a wider selection of historical data. The reasons for this include that the model is not sufficiently detailed to make this meaningful in a general context, historical data usually have many lacunae, current understanding of the mechanisms of seasonal impact on influenza is very limited, given the number and nature of chaotic solutions, fitting could be technically difficult, and last, Popper's cogent discussion of historical data-fitting [54], in which he concludes that such exercises are usually not scientifically productive, seems to be particularly pertinent to this investigation. Regrettably, in spite of all the evidence, “parameter twiddling” and fitting historic data are still highly regarded by some investigators, although it is hard to find examples where such work has led to significant progress. Nevertheless, with simplified “proof-of-concept” models, it is important that predictions should be acceptable as has been shown to be the case with the current model.

Where we part company with many influenza modellers is in our use of multiple pools to represent given categories: that is, four pools are used to represent infected latent persons and seven pools for asymptomatically infected persons (Figure 1) (but see [26]). Mathematically, this is a trivial addition requiring some extra programming, but it gives three significant benefits. First, progress through the stages of influenza is clearly sequential, suggesting that to use successive pools is biologically reasonable. Second, the overall transit times of sequential pools are gamma-distributed which is arguably more realistic than given by a single-pool representation. Third, the method opens a path to a more mechanistic picture of observations where quantities such as infectivity (14) and death rates can vary from pool to pool within a category. While simple models are best for elucidating many general principles, there seems to be no alternative to more detailed mechanistic (reductionist) models for serious application. With such models, parameters will be more determined by experiment at the assumption level, rather than making parameter adjustments on the basis of comparison of predictive-level data.

5. Conclusions

A mechanistic influenza model has been constructed in which sequential pools are used for some disease categories, allowing gamma-function-type dynamics with delays, consistent with biological observations. A simplified representation of seasonality is given and is, we believe, the first attempt to include weather explicitly in an influenza model. The model has been “validated” by application to an outbreak of influenza in a school. It has been demonstrated that seasonal forcing gives rise to a rich variety of dynamic disease patterns, from regular with outbreaks at annual, biennial, and longer intervals of time, to chaotic. Some of these predicted patterns seem highly pertinent to mankind's experiences with influenza. It is suggested that seasonality and its effects could usefully be an integral part of influenza epidemiology including the areas of prediction and amelioration. Recognizing that seasonality is important in influenza dynamics, we were surprised by our inability to find more controlled-environment studies of the effects of environmental factors on virus viability or other significant processes in the disease cycle.

Appendices

A. Sequential Pools and the Gamma Function

The aim in this appendix is to explain how the use of sequential pools to represent a given clinical category affects the transit dynamics for that category. We emphasize first that a realistic mechanistic treatment of infectious disease dynamics absolutely requires the use of several, arguably at least three, sequential pools per clinical category, and second that the traditional approach of using a single pool per clinical category cannot be expected to predict credible dynamics or robust predictions. It is to be noted that the empirical use of a gamma function for a single-pool category (e.g., the entire infected box in Figure 1) is based mechanistically on several stages (in this case four) each with first-order (exponential) kinetics. Although these points will be familiar to many workers, this appendix has been written because they are not always appreciated. A general discussion of the topic can be found in, for example, [44], pp. 818–822.

A.1. Single Pool

Assume that a given disease category, for example, “infected”, is represented by a single pool, 𝑎, so the number of pools 𝑚a=1. This method is used in many SIR models. A single pool 𝑎 emptying at rate 𝑘a from initial value of 𝑎0 at time 𝑡=0 obeys the equationsd𝑎d𝑡=𝑘a𝑎;𝑎=𝑎0e𝑘a𝑡,𝑘a=0.25day1,𝑡=0𝑎=𝑎0=1.(A.1) This is drawn in Figure 9: the line labelled “1 pool, 𝑎”. It can be seen that the maximum rate at which pool 𝑎 is depleted occurs at time 𝑡=0 (shown by ●). The mean time for departure is at 𝑡=4 day = (𝑚a=1)/𝑘a (shown by ■). For a single pool, these are far apart. The biological data do not support such dynamics, which would imply, for say the infected category, that it is most probable that an infected person exits the infected category immediately following the infection event. Since our paper is particularly concerned with the detailed dynamics of influenza, including the possible occurrence of chaotic events (here the existence of lags can make a crucial difference to model behaviour), we regard it as essential to depart from the assumption of assigning a single pool to each disease category.

A.2. Two Pools

Next assume the category is represented by two sequential pools, 𝑏1 and 𝑏2, where 𝑏1 empties into 𝑏2, and 𝑏2 is the final pool in the category. The number of pools, 𝑚b=2. In this case, the relevant equations becomed𝑏1d𝑡=𝑘b𝑏1,d𝑏2d𝑡=𝑘b𝑏1𝑘b𝑏2,𝑏1=𝑏0e𝑘b𝑡,𝑏2=𝑘b𝑡e𝑘b𝑡,𝑘b=0.5day1,𝑡=0𝑏1=𝑏0=1,𝑏2=0.(A.2) This is shown in Figure 9: the line labelled “2 pools, 𝑏2”. With two pools, 𝑘b has been doubled relative to 𝑘a in Eqns (A.1), so that the mean time for departure from the final pool in the category 𝑏2 is the same at time 𝑡=4 day = (𝑚b=2)/𝑘b (shown by ■). We see that the maximum rate for departure from the infected category now occurs at time 𝑡=2 day (= (𝑚b1)/𝑘b) (shown by ●), in contrast to 0 day for the single pool case above. In moving from one pool to two pools, the maximum value at 2 day has shifted towards the mean value at 4 day. There is a large qualitative difference between the 1-pool curve for 𝑎 and 2-pools curve for 𝑏2 drawn in Figure 9.

A.3. Three Pools

This case is given in detail as there is another qualitative difference between 2-pool and 3-pool dynamics. Now there are three sequential pools, 𝑐1, 𝑐2, and 𝑐3, with 𝑐1 emptying into 𝑐2, 𝑐2 into 𝑐3, and 𝑐3 is the last pool in the category. The number of pools is 𝑚c=3. The equations for the system ared𝑐1d𝑡=𝑘c𝑐1,d𝑐2d𝑡=𝑘c𝑐1𝑘c𝑐2,d𝑐3d𝑡=𝑘c𝑐2𝑘c𝑐3,𝑐1=𝑐0e𝑘c𝑡,𝑐2=𝑘c𝑡e𝑘c𝑡,𝑐3=𝑘c𝑡2e2!𝑘c𝑡,𝑘c=0.75day1,𝑡=0𝑐1=𝑐0𝑐=1,2=𝑐3=0.(A.3) Now (Figure 9) the curve for the final pool in the category (𝑐3) is sigmoidal at low values of time 𝑡. The time of the maximum (●) (= ((𝑚b=3)1)/𝑘c=8/3) has moved closer to the mean time (■). The mean time is (𝑚c=3)/(𝑘c=0.75) = 4 day as in the other curves. With three pools, there is now a sigmoid departure from time 𝑡=0, giving a more sharply defined biological delay (cf. the two-pool 𝑏2 line).

A.4. Four and More Pools

Adding more sequential pools to the 3-pool situation has the effect on the final pool of increasing the sigmoidicity, and moving the time of maximum (●) closer to the mean time (■). At the same time as adding more pools, the exit rate constant of each pool is increased so that the overall mean transit time is unchanged. For four pools, 𝑑1,,𝑑4, we have𝑚d=4,𝑘d=1day1,𝑑4=𝑘d𝑡𝑚d1md!e1𝑘d𝑡=𝑡36e𝑡,𝑡mean=𝑚d𝑘d=4day,𝑡max=𝑚d1𝑘d=3day.(A.4) This is drawn in Figure 9. There is only a minor quantitative difference between the 3- and 4-pool sequences when they have the same mean transit time 𝑡mean.

For a general number of pools, 𝑚, with state variables, 𝑦1,,𝑦m, and each with outgoing rate constant, 𝑘, the value of the state variable for the final pool in the sequence, 𝑦m, is𝑦m𝑘=𝑓𝑚1𝑡𝑚1(e𝑚1)!𝑘𝑡,𝑡=0𝑦1𝑦=𝑓=1,i(𝑖>1)=0,(A.5) where 𝑓 is a constant. The total outflow from the final pool, 𝑦m, is 𝑓, so that0𝑘𝑦mdt=𝑓0𝑘𝑚𝑡𝑚1e(𝑚1)!𝑘𝑡d𝑡=𝑓=1.(A.6) The statistics on the final pool, 𝑦m, are (using to denote expectation values)𝑚𝑡=𝑘,𝑡2=𝑚(𝑚+1)𝑘2,𝑡variance(𝑡)=2𝑡2=𝑚𝑘2,standarddeviation(𝑡)=𝑚variance(𝑡)=1/2𝑘,𝑡max𝑦m=ismaximum𝑚1𝑘,coecientofvariation(𝑡)=standarddeviation(𝑡)=1𝑡𝑚.(A.7) The 8-pool curve of Figure 9, illustrated by the time course of the final pool in the sequence, 8, is shown because in Figure 1 every path between susceptible and recovered traverses 8 pools each with rate constant 2 day−1, giving an overall average transit time of 4 day (see main text for further discussion of this point).

Finally we note that although there are several equivalent definitions for the gamma function used by mathematicians and others (e.g., [44], pp. 819–821), possibly the most intuitive definition for the biologist is that given in (A.6), namely, where the normalized gamma function, 𝛾(𝑡,𝑚), is given by𝑘𝛾(𝑡,𝑚)=𝑚𝑡𝑚1e(𝑚1)!𝑘𝑡,with0𝛾(𝑡,𝑚)d𝑡=0𝑘𝑚𝑡𝑚1e(𝑚1)!𝑘𝑡d𝑡=1.(A.8) If the number of sequential pools, 𝑚, and also 𝑘, with rate constant 𝑘=𝑚/𝑇, where 𝑇 is the overall mean transit time which is constant, then γ(t, m) approaches a Dirac delta function located at time 𝑡=𝑇.

B. Basic Reproductive Ratio 𝑅0 and Mean Generation Time 𝜏gen

An alternative and useful way of writing (20) for 𝑅0 is as a sum of the individual contributions from the 4 + 4 + 3 + 7 = 18 diseased pools of Figure 1. Using an obvious notation, these 18 terms are written as, first for the four infected pools,𝑅0𝛽(inf1)=𝑘vir,nm+𝑘vir,env𝑘vir𝑘inf𝑘inf𝑤+𝜇inf1𝑘inf=𝑐inf𝑤inf1𝑘inf,𝑐inf=𝛽𝑘vir,nm+𝑘vir,env𝑘vir𝑘inf𝑘inf,𝑅+𝜇0(inf2)=𝑐inf1𝑓asy𝑘inf𝑘inf𝑤+𝜇inf2𝑘inf,𝑅0(inf3)=𝑐inf1𝑓asy𝑘inf𝑘inf+𝜇2𝑤inf3𝑘inf,𝑅0(inf4)=𝑐inf1𝑓asy𝑘inf𝑘inf+𝜇3𝑤inf4𝑘inf,(B.1) next for the four sick pools𝑅0𝛽(sic1)=𝑘vir,nm+𝑘vir,env𝑘vir𝑘inf𝑘inf+𝜇1𝑓asy×𝑘inf𝑘inf+𝜇3𝑘sic𝑘sic+𝜇+𝜇sic𝑤sic1𝑘sic=𝑐sic𝑤sic1𝑘sic,𝑐sic=𝑐inf1𝑓asy𝑘inf𝑘inf+𝜇3𝑘sic𝑘sic+𝜇+𝜇sic,𝑅0(sic2)=𝑐sic1𝑓hos𝑘sic𝑘sic+𝜇+𝜇sic𝑤sic2𝑘sic,𝑅0(sic3)=𝑐sic1𝑓hos𝑘sic𝑘sic+𝜇+𝜇sic2𝑤sic3𝑘sic,𝑅0(sic4)=𝑐sic1𝑓hos𝑘sic𝑘sic+𝜇+𝜇sic3𝑤sic4𝑘sic,(B.2) then, the three hospitalized pools𝑅0𝛽(hos2)=𝑘vir,nm+𝑘vir,env𝑘vir𝑘inf𝑘inf×+𝜇1𝑓asy𝑘inf𝑘inf+𝜇3𝑘sic𝑘sic+𝜇+𝜇sic×𝑓hos𝑘hos𝑘hos+𝜇+𝜇hos𝑤hos2𝑘hos=𝑐hos𝑤hos2𝑘hos,𝑐hos=𝑐inf1𝑓asy𝑘inf𝑘inf+𝜇3×𝑘sic𝑘sic+𝜇+𝜇sic𝑓hos𝑘hos𝑘hos+𝜇+𝜇hos,𝑅0(hos3)=𝑐hos𝑘hos𝑘hos+𝜇+𝜇hos𝑤hos3𝑘hos,𝑅0(hos4)=𝑐hos𝑘hos𝑘hos+𝜇+𝜇hos2𝑤hos4𝑘hos,(B.3) and last the seven asymptomatic pools𝑅0𝛽(asy2)=𝑘vir,nm+𝑘vir,env𝑘vir𝑘inf𝑘inf+𝜇×fasy𝑘asy𝑘asy𝑤+𝜇asy2𝑘asy=𝑐asy𝑤asy2𝑘asy,𝑐asy=𝑐inf𝑓asy𝑘asy𝑘asy,𝑅+𝜇0(asy3)=𝑐asy𝑘asy𝑘asy𝑤+𝜇asy3𝑘asy,𝑅0(asy8)=𝑐asy𝑘asy𝑘asy+𝜇6𝑤asy8𝑘asy.(B.4) Basic reproductive ratio is the sum of these contributions:𝑅0=𝑅0(inf1)++𝑅0(asy8).(B.5) Generation time 𝜏gen (day) can similarly be written in terms of the contributions from the different infectious pools (using the contributions to 𝑅0 defined above). First write down the contributions to the total generation time:Δ𝜏(inf1)=𝑅01(inf1)𝑘inf+1+𝜇𝑘vir,dnm+𝑘vir,env,,Δ𝜏(inf4)=𝑅04(inf4)𝑘inf+1+𝜇𝑘vir,dnm+𝑘vir,env,Δ𝜏(sic1)=𝑅0(4sic1)𝑘inf+1+𝜇𝑘sic+𝜇+𝜇sic+1𝑘vir,dnm+𝑘vir,env,,Δ𝜏(sic4)=𝑅04(sic4)𝑘inf+4+𝜇𝑘sic+𝜇+𝜇sic+1𝑘vir,dnm+𝑘vir,env,Δ𝜏(hos2)=𝑅04(hos2)𝑘inf+1+𝜇𝑘sic+𝜇+𝜇sic+1𝑘hos+𝜇+𝜇hos+1𝑘vir,dnm+𝑘vir,env,,Δ𝜏(hos4)=𝑅04(hos4)𝑘inf+1+𝜇𝑘sic+𝜇+𝜇sic+3𝑘hos+𝜇+𝜇hos+1𝑘vir,dnm+𝑘vir,env,Δ𝜏(asy2)=𝑅01(asy2)𝑘inf+1+𝜇𝑘asy+1+𝜇𝑘vir,dnm+𝑘vir,env,,Δ𝜏(asy8)=𝑅01(asy8)𝑘inf+7+𝜇𝑘asy+1+𝜇𝑘vir,dnm+𝑘vir,env.(B.6)Finally, summing the above, the mean generation time is𝜏gen=Δ𝜏(inf1)++Δ𝜏(asy8)𝑅0.(B.7)

Conflict of Interests

The authors declare that they have no competing interests.

Acknowledgment

The work was funded, in part, through the Canada Research Chairs Program.

References

  1. E. Hak, M. J. Meijboom, and E. Buskens, “Modelling the health-economic impact of the next influenza pandemic in The Netherlands,” Vaccine, vol. 24, no. 44–46, pp. 6756–6760, 2006. View at Publisher · View at Google Scholar · View at Scopus
  2. N. A. M. Molinari, I. R. Ortega-Sanchez, M. L. Messonnier et al., “The annual impact of seasonal influenza in the US: measuring disease burden and costs,” Vaccine, vol. 25, no. 27, pp. 5086–5096, 2007. View at Publisher · View at Google Scholar · View at Scopus
  3. J. S. Nguyen-van-Tam and C. Sellwood, “Avian influenza and the threat of the next human pandemic,” Journal of Hospital Infection, vol. 65, no. 2, pp. 10–13, 2007. View at Publisher · View at Google Scholar · View at Scopus
  4. I. Capua and S. Marangon, “Control of avian influenza in poultry,” Emerging Infectious Diseases, vol. 12, no. 9, pp. 1319–1324, 2006. View at Scopus
  5. R. M. Anderson and R. M. May, Infectious Diseases of Humans: Dynamics and Control, Oxford University Press, Oxford, UK, 1991.
  6. O. Diekmann and J. Hesterbeek, Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation, Wiley, New York, NY, USA, 2000.
  7. M. Nuño, G. Chowell, and A. B. Gumel, “Assessing the role of basic control measures, antivirals and vaccine in curtailing pandemic influenza: scenarios for the US, UK and the Netherlands,” Journal of the Royal Society Interface, vol. 4, no. 14, pp. 505–521, 2007. View at Publisher · View at Google Scholar · View at Scopus
  8. M. Gauthier-Clerc, C. Lebarbenchon, and F. Thomas, “Recent expansion of highly pathogenic avian influenza H5N1: a critical review,” Ibis, vol. 149, no. 2, pp. 202–214, 2007. View at Publisher · View at Google Scholar · View at Scopus
  9. K. Van Reeth, “Avian and swine influenza viruses: our current understanding of the zoonotic risk,” Veterinary Research, vol. 38, no. 2, pp. 243–260, 2007. View at Publisher · View at Google Scholar · View at Scopus
  10. D. J. D. Earn, J. Dushoff, and S. A. Levin, “Ecology and evolution of the flu,” Trends in Ecology and Evolution, vol. 17, no. 7, pp. 334–340, 2002. View at Publisher · View at Google Scholar · View at Scopus
  11. C. M. Pease, “An evolutionary epidemiological mechanism, with applications to type a influenza,” Theoretical Population Biology, vol. 31, no. 3, pp. 422–452, 1987. View at Scopus
  12. O. A. Christophersen and A. Haug, “Why is the world so poorly prepared for a pandemic of hypervirulent avian influenza?” Microbial Ecology in Health and Disease, vol. 18, no. 3-4, pp. 113–132, 2006. View at Publisher · View at Google Scholar · View at Scopus
  13. T. Day, J. B. André, and A. Park, “The evolutionary emergence of pandemic influenza,” Proceedings of the Royal Society B: Biological Sciences, vol. 273, no. 1604, pp. 2945–2953, 2006. View at Publisher · View at Google Scholar · View at Scopus
  14. S. Iwami, Y. Takeuchi, and X. Liu, “Avian-human influenza epidemic model,” Mathematical Biosciences, vol. 207, no. 1, pp. 1–25, 2007. View at Publisher · View at Google Scholar · View at Scopus
  15. M. Nuño, Z. Feng, M. Martcheva, and C. Castillo-Chavez, “Dynamics of two-strain influenza with isolation and partial cross-immunity,” SIAM Journal on Applied Mathematics, vol. 65, no. 3, pp. 964–982, 2005. View at Publisher · View at Google Scholar · View at Scopus
  16. B. Hancioglu, D. Swigon, and G. Clermont, “A dynamical model of human immune response to influenza A virus infection,” Journal of Theoretical Biology, vol. 246, no. 1, pp. 70–86, 2007. View at Publisher · View at Google Scholar · View at Scopus
  17. G. A. Bocharov and A. A. Romanyukha, “Mathematical model of antiviral immune response III. Influenza a virus infection,” Journal of Theoretical Biology, vol. 167, no. 4, pp. 323–360, 1994. View at Publisher · View at Google Scholar · View at Scopus
  18. R. D. Balicer, M. Huerta, N. Davidovitch, and I. Grotto, “Cost-benefit of stockpiling drugs for influenza pandemic,” Emerging Infectious Diseases, vol. 11, no. 8, pp. 1280–1282, 2005. View at Scopus
  19. N. M. Ferguson, D. A. T. Cummings, S. Cauchemez et al., “Strategies for containing an emerging influenza pandemic in Southeast Asia,” Nature, vol. 437, no. 7056, pp. 209–214, 2005. View at Publisher · View at Google Scholar · View at Scopus
  20. R. Gani, H. Hughes, D. Fleming, T. Griffin, J. Medlock, and S. Leach, “Potential impact of antiviral drug use during influenza pandemic,” Emerging Infectious Diseases, vol. 11, no. 9, pp. 1355–1362, 2005. View at Scopus
  21. I. M. Longini Jr., A. Nizam, S. Xu et al., “Containing pandemic influenza at the source,” Science, vol. 309, no. 5737, pp. 1083–1087, 2005. View at Publisher · View at Google Scholar · View at Scopus
  22. B. Pourbohloul, L. A. Meyers, D. M. Skowronski, M. Krajden, D. M. Patrick, and R. C. Brunham, “Modeling control strategies of respiratory pathogens,” Emerging Infectious Diseases, vol. 11, no. 8, pp. 1249–1256, 2005. View at Scopus
  23. F. Carrat, J. Luong, H. Lao, A. V. Sallé, C. Lajaunie, and H. Wackernagel, “A “small-world-like” model for comparing interventions aimed at preventing and controlling influenza pandemics,” BMC Medicine, vol. 4, article 26, 2006. View at Publisher · View at Google Scholar · View at Scopus
  24. Y. H. Hsieh, C. C. King, C. W. S. Chen, M. S. Ho, S. B. Hsu, and Y. C. Wu, “Impact of quarantine on the 2003 SARS outbreak: a retrospective modeling study,” Journal of Theoretical Biology, vol. 244, no. 4, pp. 729–736, 2007. View at Publisher · View at Google Scholar · View at Scopus
  25. G. Chowell, C. E. Ammon, N. W. Hengartner, and J. M. Hyman, “Transmission dynamics of the great influenza pandemic of 1918 in Geneva, Switzerland: assessing the effects of hypothetical interventions,” Journal of Theoretical Biology, vol. 241, no. 2, pp. 193–204, 2006. View at Publisher · View at Google Scholar · View at Scopus
  26. M. Eichner, M. Schwehm, H. P. Duerr, and S. O. Brockmann, “The influenza pandemic preparedness planning tool InfluSim,” BMC Infectious Diseases, vol. 7, article 17, 2007. View at Publisher · View at Google Scholar · View at Scopus
  27. H. P. Duerr, S. O. Brockmann, I. Piechotowski, M. Schwehm, and M. Eichner, “Influenza pandemic intervention planning using InfluSim: pharmaceutical and non-pharmaceutical interventions,” BMC Infectious Diseases, vol. 7, article 76, 2007. View at Publisher · View at Google Scholar · View at Scopus
  28. A. L. Lloyd, “Destabilization of epidemic models with the inclusion of realistic distributions of infectious periods,” Proceedings of the Royal Society B: Biological Sciences, vol. 268, no. 1470, pp. 985–993, 2001. View at Publisher · View at Google Scholar · View at Scopus
  29. J. C. Hopkins and R. J. Leipold, “On the dangers of adjusting the parameter values of mechanism-based mathematical models,” Journal of Theoretical Biology, vol. 183, no. 4, pp. 417–427, 1996. View at Publisher · View at Google Scholar · View at Scopus
  30. S. A. Levin and R. Durrett, “From individuals to epidemics,” Philosophical Transactions of the Royal Society B: Biological Sciences, vol. 351, no. 1347, pp. 1615–1621, 1996. View at Scopus
  31. P. Sebastiani, K. D. Mandl, P. Szolovits, I. S. Kohane, and M. F. Ramoni, “A Bayesian dynamic model for influenza surveillance,” Statistics in Medicine, vol. 25, no. 11, pp. 1803–1816, 2006. View at Publisher · View at Google Scholar · View at Scopus
  32. J. Arino, F. Brauer, P. van den Driessche, J. Watmough, and J. Wu, “Simple models for containment of a pandemic,” Journal of the Royal Society Interface, vol. 3, no. 8, pp. 453–457, 2006. View at Publisher · View at Google Scholar · View at Scopus
  33. J. Shaman and M. Kohn, “Absolute humidity modulates influenza survival, transmission, and seasonality,” Proceedings of the National Academy of Sciences of the United States of America, vol. 106, no. 9, pp. 3243–3248, 2009. View at Publisher · View at Google Scholar · View at Scopus
  34. J. Dushoff, J. B. Plotkin, C. Viboud, D. J. D. Earn, and L. Simonsen, “Mortality due to influenza in the United States—an annualized regression approach using multiple-cause mortality data,” American Journal of Epidemiology, vol. 163, no. 2, pp. 181–187, 2006. View at Publisher · View at Google Scholar · View at Scopus
  35. T. A. Reichert, L. Simonsen, A. Sharma, S. A. Pardo, D. S. Fedson, and M. A. Miller, “Influenza and the winter increase in mortality in the United States, 1959–1999,” American Journal of Epidemiology, vol. 160, no. 5, pp. 492–502, 2004. View at Publisher · View at Google Scholar · View at Scopus
  36. G. A. Poland, “If you could halve the mortality rate, would you do it?” Clinical Infectious Diseases, vol. 35, no. 4, pp. 378–380, 2002. View at Publisher · View at Google Scholar · View at Scopus
  37. K. L. Nichol, J. Nordin, J. Mullooly, R. Lask, K. Fillbrandt, and M. Iwane, “Influenza vaccination and reduction in hospitalizations for cardiac disease and stroke among the elderly,” New England Journal of Medicine, vol. 348, no. 14, pp. 1322–1332, 2003. View at Publisher · View at Google Scholar · View at Scopus
  38. D. L. Crombie, D. M. Fleming, K. W. Cross, and R. J. Lancashire, “Concurrence of monthly variations of mortality related to underlying cause in Europe,” Journal of Epidemiology and Community Health, vol. 49, no. 4, pp. 373–378, 1995. View at Scopus
  39. B. M. Bolker and B. T. Grenfell, “Chaos and biological complexity in measles dynamics,” Proceedings of the Royal Society B: Biological Sciences, vol. 251, no. 1330, pp. 75–81, 1993. View at Scopus
  40. R. Casagrandi, L. Bolzoni, S. A. Levin, and V. Andreasen, “The SIRC model and influenza A,” Mathematical Biosciences, vol. 200, no. 2, pp. 152–169, 2006. View at Publisher · View at Google Scholar · View at Scopus
  41. J. Dushoff, J. B. Plotkin, S. A. Levin, and D. J. D. Earn, “Dynamical resonance can account for seasonality of influenza epidemics,” Proceedings of the National Academy of Sciences of the United States of America, vol. 101, no. 48, pp. 16915–16916, 2004. View at Publisher · View at Google Scholar · View at Scopus
  42. L. A. Rvachev and I. M. Longini, “A mathematical model for the global spread of influenza,” Mathematical Biosciences, vol. 75, no. 1, pp. 3–22, 1985. View at Scopus
  43. H. J. Wearing, P. Rohani, and M. J. Keeling, “Appropriate models for the management of infectious diseases,” PLoS Medicine, vol. 2, no. 7, Article ID e174, pp. 621–627, 2005. View at Publisher · View at Google Scholar · View at Scopus
  44. J. H. M. Thornley and J. France, Mathematical Models in Agriculture, CAB International, Wallingford, UK, 2007.
  45. C. E. Mills, J. M. Robins, and M. Lipsitch, “Transmissibility of 1918 pandemic influenza,” Nature, vol. 432, no. 7019, pp. 904–906, 2004. View at Publisher · View at Google Scholar · View at Scopus
  46. P. Wright and R. G. Webster, “Orthomyxoviruses,” in Fields Virology, D. M. Knipe, P. M. Howley, and D. E. Griffin, Eds., pp. 1533–1579, Lippincott, Williams and Wilkins, Philadelphia, Pa, USA, 2001.
  47. P. Baccam, C. Beauchemin, C. A. Macken, F. G. Hayden, and A. S. Perelson, “Kinetics of influenza A virus infection in humans,” Journal of Virology, vol. 80, no. 15, pp. 7590–7599, 2006. View at Publisher · View at Google Scholar · View at Scopus
  48. A. C. Lowen, S. Mubareka, J. Steel, and P. Palese, “Influenza virus transmission is dependent on relative humidity and temperature,” PLoS Pathogens, vol. 3, no. 10, Article ID e151, pp. 1470–1476, 2007. View at Publisher · View at Google Scholar · View at Scopus
  49. J. H. M. Thornley, Grassland Dynamics: An Ecosystem Simulation Model, CAB International, Wallingford, UK, 1998.
  50. M. P. Taylor, “Influenza 1969-1970. Incidence in general practice based on a population survey,” Journal of the Royal College of General Practitioners, vol. 21, no. 102, pp. 17–22, 1971. View at Scopus
  51. J. Fry, “Epidemic influenza. Patterns over 20 years (1949–1968),” Journal of the Royal College of General Practitioners, vol. 17, no. 79, pp. 100–103, 1969. View at Scopus
  52. B. J. Coburn, B. G. Wagner, and S. Blower, “Modeling influenza epidemics and pandemics: insights into the future of swine flu (H1N1),” BMC Medicine, vol. 7, article 30, 2009. View at Publisher · View at Google Scholar · View at Scopus
  53. Anon, “Influenza in a boarding school,” British Medical Journal, vol. 1, p. 587, 1978.
  54. K. Popper, The Poverty of Historicism, Routledge & Kegan Paul, London, UK, 1957.
  55. A. Flahault, S. Deguen, and A. J. Valleron, “A mathematical model for the European spread of influenza,” European Journal of Epidemiology, vol. 10, no. 4, pp. 471–474, 1994. View at Publisher · View at Google Scholar · View at Scopus
  56. A. Flahault, S. Letrait, P. Blin, S. Hazout, J. Menares, and A. J. Valleron, “Modelling the 1985 influenza epidemic in France,” Statistics in Medicine, vol. 7, no. 11, pp. 1147–1155, 1988. View at Scopus
  57. N. P. Johnson and J. Mueller, “Updating the accounts: global mortality of the 1918–1920 “Spanish” influenza pandemic,” Bulletin of the History of Medicine, vol. 76, no. 1, pp. 105–115, 2002. View at Scopus