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

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.


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][6][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).
Before proceeding with this agenda, a brief review of some of the issues and modelling approaches relevant to All pools other than the virus pool (n vir ) and the dead-frominfluenza pool (n 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 k inf .Similarly, k sic applies to the four nonhospitalized clinically sick pools, k hos to the three hospitalized clinically sick pools, k asy to the seven asymptomatic pools, and k del to the ten pools in the loss-of-immunity delay pipe.f asy and f hos are the fractions of the fluxes between n inf1 and n inf2 , and n sic1 and n sic2 , which are asymptomatic and hospitalized.μ sic and μ hos are the death rates (caused by flu) of the clinically sick and hospitalized pools.
influenza is given, and our contribution is discussed in this wider context.

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 todate 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][13][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.

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

Network Models and Stochastic
Simulations.Some recent influenza models are based on networks and stochastic simulations [11,[18][19][20][21][22].The model of [23] is of considerable Infecteds Persons n rec = 0 usually, but see Figure 3  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.

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.

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

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 j, j = 1, . . ., 4); clinically sick, also considered as four stages of sickness (sicj, j = 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 n del1 , n del2 , . .., to n del10 .Infectious persons give rise to virus particles (n 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 (n 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 (n vir ) would be to assume that weather affects transmission rate β directly.In the equations to follow, transmission rate β and virus pool n 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.
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 ISRN Biomathematics 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 (n city , Table 1).
2.1.1.Some Definitions.First, define some totals in terms of the state variables (Figure 1): n liv = n sus + n inf + n sic + n hos + n asy + n rec + n del , n fer = n liv = n city , n iash = n inf + n asy + n sic + n hos . ( The notation of the first five equations is obvious.Total live population n liv is given by the sixth equation and is assumed equal to the fertile population n fer (seventh equation); both of these are equal to the city population, n city ((2) and following paragraph).Finally, total infected population n iash is obtained by adding together the inf, asy, sic, and hos categories.All the variables in ( The two principal inputs are first from births with the number of fertile persons n 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 n sus pool is included, I ded .This is the death rate caused by influenza (Figure 1).I ded is given by ( 16).This is done so that the live person number, n liv (1), remains constant at value n 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 −μn sus term in ( 2)).Also, the instantaneous death rate from influenza, I 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 n sus , the virus quantity n vir and is divided by the number of live persons n liv (1).All terms in (2) scale equally with population size.

Infecteds.
The differential equations for the four sequential infected pools are ( 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 k inf n inf1 from the inf1 pool is partly asymptomatic with fraction f 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 k inf is 2 day −1 , such that with four pools the mean time from infection to clinical sickness is 4/k 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.)

Nonhospitalized Clinically Sick.
Persons who become clinically sick enter the sic1 pool (Figure 1).A fraction f 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 are 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 n inf pools above) gamma-distributed (Appendix A, Figure 9 , (A.4)).
The rate at which individuals are admitted to hospital, I hos2 (Figure 1), is with units of persons day −1 .This is often a recorded statistic.The total flu-related death flux is (input from all four sick pools to the dead-from-influenza n ded box of Figure 1; units: persons day −1 ) 2.1.5.Asymptomatic Infecteds.The differential equations for these seven pools (Figure 1) are . . .
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 are The input term f hos k sic n 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 f 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 2.1.7.Recovered and Immune.The state variable n rec is governed by the equation The positive (input) terms are from the last equations of ( 7), (4), and (8).Death occurs at its natural rate (μ).The rate constant k 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 k rec causes rapid entry into the delay sequence of pools (Figure 1) which leads to loss of immunity.

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 k del (day −1 ) is the rate constant out of each of the ten compartments, then The standard deviation in the exit times is ( √ 10)/k del and the coefficient of variation is 1/ √ 10 (= standard deviation/τ imm ) (A.7).The differential equations for the pools are The first term on the right of the 1st equation (k rec ) is the last term in (10); the second term (k 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 (k del ) is transfer into the susceptible pool (2).
2.1.9.Virus Pool.It is assumed that virus production, s vir , which is a surrogate for "infectious strength", is calculated with k vir = 200 virus units weighted infected person k vir is a production rate constant which is applied to weighted infected persons.The value above gives a basic reproductive ratio R 0 (20) of 4.78 (see Section 2.3).The weighting factors w of ( 13) are (0 ≤ w ≤ 1) 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." s vir of ( 13) is the input to the virus pool (15).The outputs from the virus pool (15) are death by natural mortality (k vir,dnm , Table 1) and death induced by a suboptimal environment, represented by rate constant k vir,env (18).The value of k 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 k vir,dnm .The differential equation for the virus pool n vir is 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 n ded pool (Figure 1) is There are no outputs.Therefore, 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.

Environmentally induced Virus Mortality k vir,env and
Seasonal Forcing.The environmentally dependent function k vir,env in (15) (Figure 1) is assumed here to depend only on daily mean values of air temperature, T 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, k vir,eny (d −1 ; Figure 1) is written as s env is used to switch environmental effects off (0) or on (1).Air temperature T air above a threshold T air0 increases virus mortality according to power q T .Parameter q T is assumed equal to two giving a quadratic dependence of temperature above the threshold T 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 Q 10 factor for the consequences of a 10 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).

Basic
Reproductive Ratio, R 0 , and the Disease Generation Time, τ gen (Day).Two important epidemiological parameters are basic reproductive ratio, R 0 , and the disease generation time, τ gen (day).These are both derived from the basic parameters of the model.R 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 R 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 R 0 can only be calculated numerically), R 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 R 0 ) 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 n vir to n inf1 and back to n vir (Figure 1).The first term is Start with a single virus particle in the n vir box (we could equally well start with a single person in the inf1 firstinfected 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/(k vir,dnm + k vir,env ).It is assumed that there is no additional virus death process due to exposure.w inf1 is a dimensionless weighting factor (14).The third factor (k vir /k inf ) is the number of virus particles produced per person in the inf1 state during his life: 1/k inf is the average lifetime of a person in the inf1 state.The last factor, with k inf + μ in the denominator, is the probability that this lifetime (of 1/k 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 R 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 n inf1 box and diverting the (primary) infected individuals into an accumulator, rather than allowing them to enter the n inf1 box, where they can lead to secondary infections.The algebraic and numerical methods agree to six decimal digits, giving a basic reproductive ratio R 0 = 4.78 for the default parameters without forcing.
A "dynamic" reproductive ratio, R dyn , can be calculated when the total population is not entirely susceptible by means of The live population, n liv , is given by (1).This allows approximately for a slowly changing fraction of susceptible ISRN Biomathematics individuals.However, R 0 may itself be dynamic on a shorter time scale due to seasonal effects on virus death rate k 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 n iash (1) increase initially at an exponential rate according to where τ (day) is the growth time constant (the proportional growth rate of infecteds is 1/τ day −1 ).This depends on both R 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.

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 seasonalforcing 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 R 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 n inf1 with a mean lifetime of 0.5 day is not infectious but the second infected pool n inf2 is, and therefore Next the latent period τ lat (day) is given by (there are four sequential pools in the infected category of Figure 1) 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 k inf , k asy , k sic , and k hos are all equal (Table 1), say to k dis ; (c) the first infected pool (n inf1 ) is assumed not infectious and all sick and hospitalized pools are infectious (14).We therefore write the infectious period τ inf as 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 k del (day −1 ) (ignoring natural mortality (12)).This gives a gamma-distributed delay (Appendix A, (A.5) to (A.8) with m = 10).τ imm and k del are related by 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 w inf2 = w 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 R 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 R 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 R 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 R 0 and therefore will underestimate the level of global control measures . . .that will be needed . ..." [52] reviews R 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 R 0 = 4.78 (20).This can be regarded as an upper baseline applying to an optimum environment (one maximizing R 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.

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

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, I hos2 (persons day −1 ; Figure 1), as this is often a recorded statistic.An "epidemic" is deemed to have occurred if I hos2 shows a maximum (in steps of Δt) and if at the maximum I hos2 ≥ 0.1 persons day −1 .
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.

Dynamics without Intrinsic Loss of Immunity and without
Forcing.Here the model is exercised with parameter k 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.

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 (n inf ) peaks a few days before hospitalized numbers (n hos ), both then falling to very low values.Most (99%) of the population joins the recovered (immune) box (n rec ), with 1% escaping infection altogether.90% of those infected travel via the asymptomatic route (Figure 1; f asy , Table 1, ( 7)).The epidemic in Figure 2(a) is short compared to UK experience,   1.There is no loss of individual immunity (k rec = 0; (10), Figure 1) and no seasonal forcing ((18); s env = 0).In each case, the value of the altered parameter, basic reproductive ratio R 0 (20), and total hospital admissions (HA, the integral of I hos2 , (5)) are given; the curve for the default parameters is thickened and indicated.(a) Changes to infection parameter β.(b) Changes to k inf (= k sic = k hos = k asy ; Figure 1).(c) Changes to the initial (time t = 0) immune fraction in the recovered category, f rec0 .This is achieved by altering the time t = 0 value of n rec state variable (Figure 1) to  Figure 4: Dynamics with loss of individual immunity.Individuals now lose immunity (k rec = 1; Figure 1).Parameters and initial conditions have the values in Table 1.(a) Short-term dynamics: susceptible n sus and delayed n del numbers are calculated with ( 2) and ( 1), the dynamic reproductive ratio R dyn with (22), and hospital admission rate I hos2 with ( 5).(b) Longer-term dynamics.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 R dyn (22) (not plotted) falling to near zero, as the fraction of susceptibles (n sus /n liv ) becomes small.In the first epidemic, at time t ≈ 20 day, the initial peak in hospital admissions I hos2 = 949 (Figure 1, ( 5)), decreasing to less than 100 at the second epidemic.
After the first epidemic, R dyn (22) slowly recovers as immune numbers (n rec ) decrease and susceptibles (n 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 R 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 k inf (with k inf = k sic = k hos = k asy ; Figure 1; Table 1; (3)), and the initial (time zero) immune fraction, f 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 k vir ; both are linear factors of R 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 longterm steady state to be more quickly attained-the spikes of I 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 (k inf = k sic = k hos = k asy = k xxx , Figure 1), is less straightforward.Here a low value for the k xxx (e.g., 0.25 day −1 ) gives a high basic reproductive ratio R 0 (20) (Appendix B) an epidemic which is slow to take hold but (surprisingly) moves more rapidly towards a longterm steady state (i.e., the I hos2 spikes as in Figure 2 15), ( 18)), giving virus inactivation by temperature alone.T air is mean daily air temperature (appropriate to southern Britain) (19); T air0 is an assumed temperature threshold for virus inactivation, here equal to 13 • C (( 18), Table 1)); R 0 is the seasonally modulated basic reproductive ratio (20), calculated assuming temperature is constant.(b) For different values of mean annual air temperature, T air,mean (19), basic reproductive ratio R 0 (20) is seasonally modulated as in (a) via the seasonal modulation of k vir,env , for a constant threshold temperature for virus inactivation of T air0 = 13 • C ( 18), giving for R 0 a maximum, minimum, annual mean, and amplitude (maximum-minimum).N 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).Increases in mean annual air temperature T air,mean increase the level of forcing (( 18), (19), Figure 6(a)).R 0 is the mean annual value of R 0 .
hospital admissions (HA).Increasing k xxx always causes the I 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 (n iash ,(1)).Finally, Figure 3(c) illustrates how increasing the initial immune fraction f 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).

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 k 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 × 365 days; 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 k rec = 1 day −1 (instead of 0) if the n rec variable of Figure 2(a) is replaced by the n 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 n city , Table 1, (1) and following paragraph), with an inverse switching of numbers in the delay compartment, n 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 R dyn = 1 (22), and the fractions in the three principal categories (Figure 1) of susceptible; all infected categories lumped together (infected, asymptomatic, sick, hospitalized, n 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.

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 ISRN Biomathematics 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 homogeneousmixing 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 n city = 763; and the asymptomatic fraction, f 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 R 0 (20) is 6.44, the mean generation time is 2.65 day (A.7), and the initial proportional growth rate of total infecteds, n iash(1) , is 1.015 day −1 .
For comparison with the data in [53], we define the number of persons confined to bed, n ctb , as Since this definition is somewhat arbitrary, in comparing the predicted n ctb from the model with the data from Anon (1978), n ctb (t) is scaled with an adjustable factor so that the comparison line drawn in Figure 5 is 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 m pools for the E category and n pools for the I category.They find a best fit with m = 2 and n = 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 R 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.

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 (T air ), relative humidity (R H ), radiation (possibly multicomponent), and wind speed.The effects can be highly complex: for example, [48] which examines the effects of T air and R 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 T air alone is operative.
Figure 6(a) illustrates the seasonal variation of mean daily air temperature T air in the southern UK.Using (19) with (18), and s env = 1, this affects environmentally induced virus death rate k vir,env (Figure 1, ( 15)), and thereby basic reproductive ratio R 0 (20).Mean daily air temperature T air varies between 3 • C (Jan 24) and 17 • C (Jul 25).When T air crosses the temperature threshold of T air0 = 13 • C in May, this increases environmentally induced virus mortality, k vir,env (18) and decreases basic reproductive ratio R 0 (20).With this formulation, influenza is most likely during the months of October to April when R 0 is highest.
Figure 6(b) shows how varying (e.g., increasing) mean annual air temperature, T air,mean (19), (or equivalently, decreasing threshold temperature T air0 , ( 18)) varies the duration and intensity of seasonal forcing of R 0 .Changing mean annual air temperature, T air,mean (19), changes the average annual value of R 0 , R 0 , as well as its maximum, minimum and amplitude (maximum-minimum).The dependence of these quantities on T ait,mean is also shown, together with the number of forcing days per year (N 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 T air,mean ≤ 6 • C, there is no forcing at all (N fdy = 0): air temperature T air (19) never exceeds threshold temperature T air0 = 13 • C (18); R 0 is invariant at its maximum value.Environmentally induced virus death rate k vir,env (( 15), ( 18)) stays at zero; the system remains in the steady state shown in Figure 4(b) and R 0 = R 0 (this situation is equivalent to T air,mean = 10 • C and T air0 = 17 • C).As mean annual air temperature T air,mean increases above 6 • C, there are more days in the summer months when mean daily temperature T air > threshold temperature T air0 , number of forcing days per year N fdy increases, and R 0 < R 0 .The amplitude increases to a maximum whenT air,mean = 21 • C, before decreasing.With T air,mean > 28.5 • C, influenza epidemics do not occur, because the annual average of R 0 is less than unity.
In the simulations presented in Figures 7 and 8, various values of mean annual air temperature T air,mean are taken between −4.5 and 16 • C, assuming always an annual variation T air,var of ±7 • C (19) as in the UK.Otherwise, parameters have the values in Table 1.In these simulations, the longterm steady state reached in Figure 4, with loss of immunity (k 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 (maximumminimum) 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 Figure 7 illustrates the effects of low levels of seasonal forcing on influenza hospital admissions, I hos2 (Figure 1,( 5)).Forcing is seen as a downward modulation of the basic reproductive ratio R 0 , which decreases the annual average, R 0 .
The lower graph in Figure 7(a) shows that R 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 day −1 ) in a 12month 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 twelvemonth 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 5week 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 T 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); T 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); T 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); T 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 T air,mean is increased, forcing is increased (i.e., the magnitude of the seasonal changes in basic reproductive ratio R 0 ((20), ( 18), (19), Figure 6) increases), but the mean annual value of R 0 , R 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 T air,mean can be continued, and although the situations simulated are now less realistic, they do contribute to understanding the system.With T 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 R 0 is 2.945.When T 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 R 0 of 2.499.When T air,mean = 26 • C the response is chaotic with 15 epidemics in the first 200 years and R 0 = 1.367.Last, with T 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.R 0 has a modest mean annual value (1.131) and is below unity for much of the year.

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

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.

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, a, so the number of pools m a = 1.This method is used in many SIR models.A single pool a emptying at rate k a from initial value of a 0 at time t = 0 obeys the equations This is drawn in Figure 9 : the line labelled "1 pool, a".It can be seen that the maximum rate at which pool a is depleted occurs at time t = 0 (shown by •).The mean time for departure is at t = 4 day = (m a = 1)/k 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.We see that the maximum rate for departure from the infected category now occurs at time t = 2 day (= (m b − 1)/k 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 a and 2-pools curve for b 2 drawn in Figure 9.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 t mean .
For a general number of pools, m, with state variables, y 1 , . . ., y m , and each with outgoing rate constant, k, the value of the state variable for the final pool in the sequence, y m , is y m = f k m−1 t m−1 (m − 1)! e −kt , t = 0 : y 1 = f = 1, where f is a constant.The total outflow from the final pool, y m , is f, so that The statistics on the final pool, y m , are (using to denote expectation values) (A.7) The 8-pool curve of Figure 9 , illustrated by the time course of the final pool in the sequence, h 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, γ(t, m), is given by  Generation time τ gen (day) can similarly be written in terms of the contributions from the different infectious pools (using the contributions to R 0 defined above).First write down the contributions to the total generation time:

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

Figure 2 :
Figure 2: Dynamics without loss of individual immunity or seasonal forcing.There is no loss of individual immunity (k rec = 0; (10), Figure 1) nor seasonal forcing ((18); s env = 0).(a) Short-term dynamics.(b) Long-term dynamics.Parameters and initial conditions have the values in Table1.Successive epidemics are caused by loss of immunity at the population level from the natural birth and death processes (μ, (2), Table1).

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

Figure 6 :
Figure 6: Environmental effects on basic reproductive ratio, R 0 .(a) Seasonal change in virus inactivation rate, k vir,env (Figure 1; (15), (18)), giving virus inactivation by temperature alone.T air is mean daily air temperature (appropriate to southern Britain)(19); T air0 is an assumed temperature threshold for virus inactivation, here equal to 13 • C ((18), Table1)); R 0 is the seasonally modulated basic reproductive ratio(20), calculated assuming temperature is constant.(b) For different values of mean annual air temperature, T air,mean(19), basic reproductive ratio R 0(20) is seasonally modulated as in (a) via the seasonal modulation of k vir,env , for a constant threshold temperature for virus inactivation of T air0 = 13 • C (18), giving for R 0 a maximum, minimum, annual mean, and amplitude (maximum-minimum).N 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).

5 T 5 TTFigure 7 : 5 TT 5 T
Figure7: Weak environmental forcing.Seasonal forcing(18) combined with loss of individual immunity (Figure4; Figure1, k rec = 1 day −1 ).The system is initially in a steady state without forcing (Figure4(b)); forcing is applied after 1 year.The effects of seasonal forcing on the basic reproductive rate R 0(20) are shown in parts (a)-(c) (right-hand y-axis).Increases in mean annual air temperature T air,mean increase the level of forcing ((18),(19), Figure6(a)).R 0 is the mean annual value of R 0 .Parts (c) and (d) are from a single simulation and have the same level of forcing.

Figure 8 :
Figure 8: Moderate environmental forcing.Seasonal forcing(18) combined with loss of individual immunity (Figure4; Figure1, k rec = 1 day −1 ).The system is initially in a steady state without forcing (Figure4(b)); forcing is applied after 1 year.The effects of seasonal forcing on the basic reproductive rate R 0(20) are shown (right-hand y-axis).Increases in mean annual air temperature T air,mean increase the level of forcing ((18),(19), Figure6(a)).R 0 is the mean annual value of R 0 .

d 4 8 pools, h 8 2 pools, b 2 Figure 9 :
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 (b 2 ), three pools (c 3 ), four pools (d 4 ), and eight pools (h 8 ) are given, showing in each case the value of the state variable for the final pool in the sequence.For m pools, the rate constant k m has been set to give the same mean overall transit time of m/k m = 4 day (indicated by ).For each curve, the maximum value (•) denotes the most probable overall transit time.
Figure6(a) illustrates the seasonal variation of mean daily air temperature T air in the southern UK.Using(19) with(18), and s env = 1, this affects environmentally induced virus death rate k vir,env (Figure1,(15)), and thereby basic reproductive ratio R 0(20).Mean daily air temperature T air varies between 3 • C (Jan 24) and 17 • C (Jul 25).When T air crosses the temperature threshold of T air0 = 13 • C in May, this increases environmentally induced virus mortality, k vir,env(18) and decreases basic reproductive ratio R 0(20).With this formulation, influenza is most likely during the months of October to April when R 0 is highest.Figure6(b) shows how varying (e.g., increasing) mean annual air temperature, T air,mean(19), (or equivalently, decreasing threshold temperature T air0 , (18)) varies the duration and intensity of seasonal forcing of R 0 .Changing mean annual air temperature, T air,mean(19), changes the average annual value of R 0 , R 0 , as well as its maximum, minimum and amplitude (maximum-minimum).The dependence of these quantities on T ait,mean is also shown, together with the number of forcing days per year (N 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 T air,mean ≤ 6 • C, there is no forcing at all (N fdy = 0): air temperature T air(19) never exceeds threshold temperature T air0 = 13 • C (18); R 0 is invariant at its maximum value.Environmentally induced virus death rate k vir,env ((15), (18)) stays at zero; the system remains in the steady state shown in Figure4(b) and R 0 = R 0 (this situation is equivalent to T air,mean = 10 • C and T air0 = 17 • C).As mean annual air temperature T air,mean increases above 6 • C, there are more days in the summer months when mean daily temperature T air > threshold temperature T air0 , number of forcing days per year N fdy increases, and R 0 < R 0 .The amplitude increases to a maximum whenT air,mean = 21 • C, before decreasing.With T air,mean > 28.5 • C, influenza epidemics do not occur, because the annual average of R 0 is less than unity.In the simulations presented in Figures7 and 8, various values of mean annual air temperature T air,mean are taken between −4.5 and 16 • C, assuming always an annual variation T air,var of ±7 • C (19) as in the UK.Otherwise, parameters have the values in Table1.In these simulations, the longterm steady state reached in Figure4, with loss of immunity (k rec = 1; Figure1; (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 (maximumminimum) of basic reproductive ratio can be read off Figure6(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) (Figure6).Therefore, of the simulations described in Figures7 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.Figure7illustrates the effects of low levels of seasonal forcing on influenza hospital admissions, I hos2 (Figure1,(5)).Forcing is seen as a downward modulation of the basic

A. 2 .
Two Pools.Next assume the category is represented by two sequential pools, b 1 and b 2 , where b 1 empties into b 2 , and b 2 is the final pool in the category.The number of pools, m b = 2.In this case, the relevant equations become db 1 dt = −k b b 1 , db 2 dt = k b b 1 − k b b 2 , b 1 = b 0 e −kbt , b 2 = k b te −kbt , k b = 0.5 day −1 , t = 0 : b 1 = b 0 = 1, b 2 = 0. (A.2) This is shown in Figure 9 : the line labelled "2 pools, b 2 ".With two pools, k b has been doubled relative to k a in Eqns (A.1), so that the mean time for departure from the final pool in the category b 2 is the same at time t = 4 day = (m b = 2)/k b (shown by ).

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, c 1 , c 2 , and c 3 , with c 1 emptying into c 2 , c 2 into c 3 , and c 3 is the last pool in the category.The number of pools is m c = 3.The equations for the system aredc 1 dt = −k c c 1 , dc 2 dt = k c c 1 − k c c 2 , dc 3 dt = k c c 2 − k c c 3 , c 1 = c 0 e −kct , c 2 = k c te −kct , c 3 = (k c t) 2 2! e −kct , k c = 0.75 day −1 , t = 0 : c 1 = c 0 = 1,c 2 = c 3 = 0. (A.3) Now (Figure 9) the curve for the final pool in the category (c 3 ) is sigmoidal at low values of time t.The time of the maximum (•) (= ((m b = 3) − 1)/k c = 8/3) has moved closer to the mean time ( ).The mean time is (m c = 3)/(k c = 0.75) = 4 day as in the other curves.With three pools, there is now a sigmoid departure from time t = 0, giving a more sharply defined biological delay (cf. the two-pool b 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, d 1 , . . ., d 4 , we have m d = 4, k d = 1 day −1 , d 4 = (k d t) md−1 (m d − 1)! e −kdt = t 3 6 e −t , t mean = m d k d = 4 day, t max = m d − 1 k d = 3 day.(A.4)

8 )
If the number of sequential pools, m → ∞, and also k → ∞, with rate constant k = m/T, where T is the overall mean transit time which is constant, then γ(t, m) approaches a Dirac delta function located at time t = T.B. Basic Reproductive Ratio R 0 and Mean Generation Time τ genAn alternative and useful way of writing(20) for R 0 is as a sum of the individual contributions from the 4 + 4 + 3 + 7 = 18 diseased pools of Figure1.Using an obvious notation, these 18 terms are written as, first for the four infected pools,R 0 (inf1) = β k vir,nm + k vir,env k vir k inf k inf + μ w inf1 k inf = c inf w inf1 k inf ,

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.State variable initial values (t = 0) Units n sus = n city − n inf1 − n rec Susceptibles Persons n inf1 = 1 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.
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 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 Tair 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 T air varies from c. 3 (January) to 17 • C (July) ([44], p. 270; [49] p. 142).It is assumed that T air varies sinusoidally, with (18)temperature rise on chemical reaction rate, or the application of the Arrhenius equation for chemical reactions (e.g.,[44], pp.103-105).k•C(varied),Tair,var=7 • C, τ Tair = 115 day.(19)Annualmean and seasonal variation are T air,mean and T air,vr .iJulian is the Julian day number (1 on 1 January, leap years are ignored).τTair(d)isthe phase of the sinusoid, which is maximum on 25 July.Combining equations(19)with(18)modifies environmentally induced virus death rate, k vir,env , in Table1.Successive epidemics are caused by loss of immunity at the population level from the natural birth and death processes (μ, (2), Table1).
Figure 3: Effect of varying key parameters.Parameters and initial conditions have the values in Table