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

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 [5–7]. 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).

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 [12–14].

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, 18–22]. 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⁻¹ day⁻¹ (rather than the customary day⁻¹) (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⁻¹, and the modulation of virus mortality ((15), (18)) may be considered as modulating the traditional transmission rate.

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⁻¹ 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):𝑛𝑛_infsic=𝑛_inf1+𝑛_inf2+𝑛_inf3+𝑛_inf4,𝑛=_hos𝑛_sic1+𝑛_sic2+𝑛_sic3+𝑛_sic4,𝑛_asy=𝑛_hos2+𝑛_hos3+𝑛_hos4,𝑛_del=𝑛_asy2+𝑛_asy3+⋯+𝑛_asy7+𝑛_asy8,𝑛_liv=𝑛_del1+𝑛_del2+⋯+𝑛_del9+𝑛_del10,𝑛_fer=𝑛_sus+𝑛_inf+𝑛_sic+𝑛_hos+𝑛_asy+𝑛_rec+𝑛_del,=𝑛_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𝑛_inf1d𝑡=𝛽𝑛^sus𝑛_vir𝑛_livd𝑛_inf2−𝑘_inf𝑛_inf1−𝜇𝑛_inf1,d𝑡=1−𝑓^asy𝑘_infd𝑛_inf3𝑛_inf1−𝑘_inf𝑛_inf2−𝜇𝑛_inf2,d𝑡=𝑘^infd𝑛_inf4𝑛_inf2−𝑘_inf𝑛_inf3−𝜇𝑛_inf3,d𝑡=𝑘^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⁻¹, 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𝑡=𝑘^infd𝑛_sic2𝑛_inf4−𝑘_sic𝑛_sic1−𝜇+𝜇_sic𝑛_sic1,d𝑡=1−𝑓^hos𝑘_sicd𝑛_sic3𝑛_sic1−𝑘_sic𝑛_sic2−𝜇+𝜇_sic𝑛_sic2,d𝑡=𝑘^sicd𝑛_sic4𝑛_sic2−𝑘_sic𝑛_sic3−𝜇+𝜇_sic𝑛_sic3,d𝑡=𝑘^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 day⁻¹. 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 day⁻¹)𝐼_sic→ded=𝜇_sic𝑛_sic1+𝑛_sic2+𝑛_sic3+𝑛_sic4.(6)

2.1.5. Asymptomatic Infecteds

The differential equations for these seven pools (Figure 1) ared𝑛_asy2d𝑡=𝑓^asyd𝑛_asy3𝑘_inf𝑛_inf1−𝑘_asy𝑛_asy2−𝜇𝑛_asy2,d𝑡=𝑘^asyd𝑛…_asy8𝑛_asy2−𝑘_asy𝑛_asy3−𝜇𝑛_asy3,d𝑡=𝑘^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𝑡=𝑓^hosd𝑛_hos3𝑘_sic𝑛_sic1−𝑘_hos𝑛_hos2−𝜇+𝜇_hos𝑛_hos2,d𝑡=𝑘^hosd𝑛_hos4𝑛_hos2−𝑘_hos𝑛_hos3−𝜇+𝜇_hos𝑛_hos3,d𝑡=𝑘^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𝐼_hos→ded=𝜇_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⁻¹, Table 1, the precise value is unimportant as long as it is of order of 1 day⁻¹ 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⁻¹) 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𝑡=𝑘^recd𝑛_del2𝑛_rec−𝑘_del𝑛_del1−𝜇𝑛_del1,d𝑡=𝑘^deld𝑛…_del9𝑛_del1−𝑘_del𝑛_del2−𝜇𝑛_del2,d𝑡=𝑘^deld𝑛_del10𝑛_del8−𝑘_del𝑛_del9−𝜇𝑛_del9,d𝑡=𝑘^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𝑛_asy8k_vir+𝑤_hos2𝑛_hos2+𝑤_hos3𝑛_hos3+𝑤_hos4𝑛_hos4.=200virusunits(weightedinfectedperson)⁻¹day⁻¹.(13)𝑘_vir is a production rate constant which is applied to weighted infected persons. The value above gives a basic reproductive ratio 𝑅₀ (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=𝐼_sic→ded+𝐼_hos→ded.(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 (d⁻¹; Figure 1) is written as𝑘_vir,env=𝑠_env𝑘_Tair𝑇^air−𝑇_air0+ABS𝑇_air−𝑇_air02^𝑞𝑇𝑠_env,=0(𝑇or_air01),𝑘_Tair=0.1^∘C⁻²d⁻¹,=13^∘C,𝑞_𝑇=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 𝑄₁₀ 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,varsin2𝜋365𝑖^Julian𝑇−𝜏_Tair_air,mean;=10𝜏^∘C(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, R₀, and the Disease Generation Time, 𝜏_gen (Day)

Two important epidemiological parameters are basic reproductive ratio, 𝑅₀, and the disease generation time, 𝜏_gen (day). These are both derived from the basic parameters of the model. 𝑅₀ 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 𝑅₀ 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 𝑅₀ can only be calculated numerically), 𝑅₀ 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 𝑅₀) 𝑅₀=𝛽𝑘_vir,nm+𝑘_vir,env𝑘_vir𝑘_inf×𝑤𝑘_inf+𝜇_inf1𝑘_inf×𝑘+_inf1−𝑓_asy𝑘_inf+𝜇𝑤𝑘^inf2_inf+𝑘^inf𝑘_inf+𝜇²𝑤_inf3𝑘_inf+𝑘^inf𝑘_inf+𝜇³×𝑤^inf4𝑘_inf+𝑘^sic𝑘_sic+𝜇+𝜇_sic×𝑤^sic1𝑘_sic+1−𝑓_hos𝑘^sic𝑘_sic+𝜇+𝜇_sic𝑤_sic2𝑘_sic+𝑘^sic𝑘_sic+𝜇+𝜇_sic²𝑤_sic3𝑘_sic+𝑘^sic𝑘_sic+𝜇+𝜇_sic³𝑤_sic4+𝑓_hos𝑘_sic𝑘_hos+𝑘𝜇^hos+𝜇_hos𝑤_hos2𝑘_hos+𝑘^hos𝑘_hos+𝜇+𝜇_hos²𝑤_hos3𝑘_hos+𝑘^hos𝑘_hos+𝜇+𝜇_hos³𝑤_hos4+𝑓_asy𝑘_hos𝑘_asy𝑘_asy+𝜇𝑤𝑘^asy2_asy+𝑘^asy𝑘_asy+𝜇𝑤𝑘^asy3_asy+𝑘^asy𝑘_asy+𝜇²𝑤_asy4+⋯+𝑘^asy𝑘_asy+𝜇⁶𝑤_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 𝑅₀ 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 𝑅₀=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=𝑅₀𝑛_sus𝑛_liv.(22) The live population, 𝑛_liv, is given by (1). This allows approximately for a slowly changing fraction of susceptible individuals. However, 𝑅₀ 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𝑛_iash∼e^𝑡/𝜏,(23) where 𝜏 (day) is the growth time constant (the proportional growth rate of infecteds is 1/𝜏 day−1). This depends on both 𝑅₀ 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 𝑅₀ (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=2day⁻¹=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=2day⁻¹=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=2day⁻¹=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⁻¹) (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 𝑅₀ (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 𝑅₀ 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 𝑅₀ = 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 𝑅₀ and therefore will underestimate the level of global control measures … that will be needed ….” [52] reviews 𝑅₀ 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 𝑅₀ = 4.78 (20). This can be regarded as an upper baseline applying to an optimum environment (one maximizing 𝑅₀), 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⁻¹); 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 𝐼_hos2≥0.1 persons day⁻¹.

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.

(a)

(b)

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).

(a)

(b)

(c)

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)).

(a)

(b)

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.

Figure 5 

Model validation. The model is validated by comparing predictions (solid line) with data (●) from [53]. See text for details.

(a)

(b)

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).

(a)

(b)

(c)

(d)

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.

(a)

(b)

(c)

(d)

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.

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 𝑅₀ (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⁻¹) gives a high basic reproductive ratio 𝑅₀ (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 𝑅₀ decreases (20) first the epidemic becomes narrower and faster (e.g., 𝑘_xxx=2 day⁻¹), but further increases in 𝑘_xxx (e.g., 𝑘_xxx=5 and 8 day⁻¹) 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⁻¹, 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⁻¹ (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)⁻¹ day⁻¹; 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 𝑅₀ (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(_hos2𝑡)+𝑛_sic2(𝑡)+n_sic3(𝑡)+𝑛_sic4(𝑡)(𝑡)+𝑛_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 𝑅₀ 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 𝑅₀ (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 𝑅₀ (20). With this formulation, influenza is most likely during the months of October to April when 𝑅₀ 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 𝑅₀. Changing mean annual air temperature, 𝑇_air,mean (19), changes the average annual value of 𝑅₀, 𝑅₀, 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,mean≤6^∘C, there is no forcing at all (𝑁_fdy=0): air temperature 𝑇_air (19) never exceeds threshold temperature 𝑇_air0 = 13°C (18); 𝑅₀ 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 𝑅₀=𝑅₀ (this situation is equivalent to 𝑇_air,mean=10^∘C 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 𝑅₀<𝑅₀. The amplitude increases to a maximum when𝑇_air,mean=21^∘C, before decreasing. With 𝑇_air,mean>28.5^∘C, influenza epidemics do not occur, because the annual average of 𝑅₀ 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 𝑅₀, which decreases the annual average, 𝑅₀.

The lower graph in Figure 7(a) shows that 𝑅₀ 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 day⁻¹) 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=11^∘C (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=12^∘C), 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=13^∘C), there is chaos but now with a tendency towards biennial epidemics (169 epidemics occur during 200 years). Last (Figure 8(d); 𝑇_air,mean=15^∘C), 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 𝑅₀ ((20), (18), (19), Figure 6) increases), but the mean annual value of 𝑅₀, 𝑅₀, 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=19^∘C, 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 𝑅₀ is 2.945. When 𝑇_air,mean=21^∘C, the response eventually settles down to a regular pattern with 27 epidemics in the first 200 years and eight years between epidemics and an 𝑅₀ of 2.499. When 𝑇_air,mean=26^∘C the response is chaotic with 15 epidemics in the first 200 years and 𝑅₀=1.367. Last, with 𝑇_air,mean=27.5^∘C, 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. 𝑅₀ 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 𝑎₀ at time 𝑡=0 obeys the equationsd𝑎d𝑡=−𝑘^a𝑘_a𝑎;𝑎=𝑎₀e^−𝑘a𝑡,=0.25day⁻¹,𝑡=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, 𝑏₁ and 𝑏₂, where 𝑏₁ empties into 𝑏₂, and 𝑏₂ is the final pool in the category. The number of pools, 𝑚_b=2. In this case, the relevant equations becomed𝑏₁d𝑡=−𝑘^b𝑏₁,d𝑏²d𝑡=𝑘^b𝑏₁𝑏₁−𝑘_b𝑏₂,𝑘_b=𝑏₀e^−𝑘b𝑡,𝑏₂=𝑘_b𝑡e^−𝑘b𝑡,=0.5day⁻¹,𝑡=0∶𝑏₁=𝑏₀=1,𝑏₂=0.(A.2) This is shown in Figure 9: the line labelled “2 pools, 𝑏₂”. 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 𝑏₂ 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 (= (𝑚_b−1)/𝑘_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 𝑏₂ 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, 𝑐₁, 𝑐₂, and 𝑐₃, with 𝑐₁ emptying into 𝑐₂, 𝑐₂ into 𝑐₃, and 𝑐₃ is the last pool in the category. The number of pools is 𝑚_c=3. The equations for the system ared𝑐₁d𝑡=−𝑘^c𝑐₁,d𝑐²d𝑡=𝑘^cd𝑐₃𝑐₁−𝑘_c𝑐₂,d𝑡=𝑘^c𝑐₁𝑐₂−𝑘_c𝑐₃,=𝑐₀e^−𝑘c𝑡,𝑐₂=𝑘_c𝑡e^−𝑘c𝑡,𝑐₃=𝑘^c𝑡²2!e^−𝑘c𝑡𝑘_c,=0.75day𝑐⁻¹,𝑡=0∶𝑐₁=𝑐₀=1,₂=𝑐₃=0.(A.3) Now (Figure 9) the curve for the final pool in the category (𝑐₃) 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 𝑏₂ 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, 𝑑₁,…,𝑑₄, we have𝑚_d=4,𝑘_d=1day⁻¹,𝑑₄=𝑘^d𝑡^𝑚d−1m_d−1!e^−𝑘d𝑡=𝑡³6e^−𝑡𝑡_mean,=𝑚^d𝑘_d=4day,𝑡_max=𝑚^d−1𝑘_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, 𝑦₁,…,𝑦_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(𝑚−1)!e^−𝑘𝑡𝑦,𝑡=0∶𝑦₁=𝑓=1,_i(𝑖>1)=0,(A.5) where 𝑓 is a constant. The total outflow from the final pool, 𝑦_m, is 𝑓, so that^∞0𝑘𝑦_mdt=𝑓^∞0𝑘^𝑚𝑡^𝑚−1(𝑚−1)!e^−𝑘𝑡d𝑡=𝑓=1.(A.6) The statistics on the final pool, 𝑦_m, are (using ⟨⟩ to denote expectation values)⟨𝑡⟩=𝑚𝑘,𝑡²=𝑚(𝑚𝑘²+1)variance(𝑡)=𝑡₂,−⟨𝑡⟩²=𝑚𝑘²standarddeviation(𝑡)=√,variance(𝑡)=𝑚^1/2𝑡_max𝑘,coefficientofvaria𝑦ti^monis(𝑡m)=axismtaunmdard=d𝑚ev𝑘−iat1io,n(𝑡)⟨𝑡⟩=√1𝑚.(A.7) The 8-pool curve of Figure 9, illustrated by the time course of the final pool in the sequence, ℎ₈, is shown because in Figure 1 every path between susceptible and recovered traverses 8 pools each with rate constant 2 day⁻¹, 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𝛾(𝑡,𝑚)=𝑘^𝑚𝑡^𝑚−1(𝑚−1)!e^−𝑘𝑡with_∞0,𝛾(𝑡,𝑚)d𝑡=^∞0𝑘^𝑚𝑡^𝑚−1(𝑚−1)!e^−𝑘𝑡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 𝑅₀ and Mean Generation Time 𝜏_gen

An alternative and useful way of writing (20) for 𝑅₀ 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,𝑅₀(inf1)=𝛽𝑘_vir,nm+𝑘_vir,env𝑘_vir𝑘_inf𝑘_inf+𝜇𝑤𝑘^inf1_inf=𝑐_inf𝑤_inf1𝑘_inf𝑐_inf,=𝛽𝑘_vir,nm+𝑘_vir,env𝑘_vir𝑘_inf𝑅𝑘_inf+𝜇₀,(inf2)=𝑐_inf1−𝑓_asy𝑘^inf𝑘_inf+𝜇𝑤𝑘^inf2_inf𝑅₀,(inf3)=𝑐_inf1−𝑓_asy𝑘^inf𝑘_inf+𝜇²𝑤_inf3𝑘_inf𝑅₀,(inf4)=𝑐_inf1−𝑓_asy𝑘^inf𝑘_inf+𝜇³𝑤_inf4𝑘_inf,(B.1) next for the four sick pools𝑅₀(sic1)=𝛽𝑘_vir,nm+𝑘_vir,env𝑘_vir𝑘_inf×𝑘_inf𝑘_inf+𝜇1−𝑓^asy𝑘_inf+𝜇³𝑘_sic𝑘_sic+𝜇+𝜇_sic𝑤_sic1𝑘_sic=𝑐_sic𝑤_sic1𝑘_sic𝑐_sic,=𝑐_inf1−𝑓_asy𝑘^inf𝑘_inf+𝜇³𝑘_sic𝑘_sic+𝜇+𝜇_sic𝑅₀,(sic2)=𝑐_sic1−𝑓_hos𝑘^sic𝑘_sic+𝜇+𝜇_sic𝑤_sic2𝑘_sic𝑅₀,(sic3)=𝑐_sic1−𝑓_hos𝑘^sic𝑘_sic+𝜇+𝜇_sic²𝑤_sic3𝑘_sic𝑅₀,(sic4)=𝑐_sic1−𝑓_hos𝑘^sic𝑘_sic+𝜇+𝜇_sic³𝑤_sic4𝑘_sic,(B.2) then, the three hospitalized pools𝑅₀(hos2)=𝛽𝑘_vir,nm+𝑘_vir,env𝑘_vir𝑘_inf×1−𝑓_asy𝑘_inf+𝜇𝑘𝑘_inf+^inf𝜇³𝑘_sic𝑘_sic+𝜇+𝜇_sic×𝑓_hos𝑘_hos𝑘_hos+𝜇+𝜇_hos𝑤_hos2𝑘_hos=𝑐_hos𝑤_hos2𝑘_hos𝑐_hos,=𝑐_inf1−𝑓_asy𝑘^inf𝑘_inf+𝜇³×𝑘^sic𝑘_sic+𝜇+𝜇_sic𝑓_hos𝑘_hos𝑘_hos+𝜇+𝜇_hos𝑅₀,(hos3)=𝑐_hos𝑘_hos𝑘_hos+𝜇+𝜇_hos𝑤_hos3𝑘_hos𝑅₀,(hos4)=𝑐_hos𝑘_hos+𝑘𝜇^hos+𝜇_hos²𝑤_hos4𝑘_hos,(B.3) and last the seven asymptomatic pools𝑅₀(asy2)=𝛽𝑘_vir,nm+𝑘_vir,env𝑘_vir𝑘_inf×f_asy𝑘_inf+𝜇𝑘_asy𝑘_asy+𝜇𝑤𝑘^asy2_asy=𝑐_asy𝑤_asy2𝑘_asy𝑐_asy,=𝑐_inf𝑓_asy𝑘_asy𝑅𝑘_asy+𝜇₀,(asy3)=𝑐_asy𝑘_asy𝑘_asy+𝜇𝑤𝑘^asy3_asy𝑅₀…,(asy8)=𝑐_asy𝑘𝑘_asy+^asy𝜇⁶𝑤_asy8𝑘_asy.(B.4) Basic reproductive ratio is the sum of these contributions:𝑅₀=𝑅₀(inf1)+…+𝑅₀(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 𝑅₀ defined above). First write down the contributions to the total generation time:Δ𝜏(inf1)=𝑅₀(inf1)1Δ𝜏(inf4)=…𝑅₀,𝑘_inf+𝜇+𝑘_vir,dnm1+𝑘_vir,env,(inf4)4Δ𝜏(sic1)=𝑅₀𝑘_inf+𝜇+𝑘_vir,dnm1+𝑘_vir,env,(sic1)4𝑘_inf+𝜇+𝑘_sic+𝜇1+𝜇_sic+1Δ𝜏(sic4)=𝑅₀𝑘_vir,dnm+𝑘_vir,env,…,(sic4)4𝑘_inf+𝜇+𝑘_sic+𝜇4+𝜇_sic+1Δ𝜏(hos2)=𝑅₀𝑘_vir,dnm+𝑘_vir,env,(hos2)4𝑘_inf+𝜇+𝑘_sic+𝜇1+𝜇_sic+1𝑘_hos+𝜇+𝜇_hos+1Δ𝜏(hos4)=…𝑅₀,𝑘_vir,dnm+𝑘_vir,env,(hos4)4𝑘_inf+𝜇+𝑘_sic+𝜇1+𝜇_sic+3𝑘_hos+𝜇+𝜇_hos+1Δ𝜏(asy2)=𝑅₀𝑘_vir,dnm+𝑘_vir,env,(asy2)1Δ𝜏(asy8)=𝑅₀𝑘_inf+𝜇+𝑘_asy1+𝜇+𝑘_vir,dnm1+𝑘_vir,env,…,(asy8)1𝑘_inf+𝜇+𝑘_asy7+𝜇+𝑘_vir,dnm1+𝑘_vir,env.(B.6)Finally, summing the above, the mean generation time is𝜏_gen=∑Δ𝜏(inf1)+⋯+Δ𝜏(asy8)𝑅₀.(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.