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).
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 (); clinically sick, also considered as four stages of sickness (); 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 (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 , to . Infectious persons give rise to virus particles (), 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 (). 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 () would be to assume that weather affects transmission rate Ξ² directly. In the equations to follow, transmission rate Ξ² and virus pool 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 2 dayβ1 to the four rate constants , and (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 (, Table 1).
2.1.1. Some Definitions
First, define some totals in terms of the state variables (Figure 1): The notation of the first five equations is obvious. Total live population is given by the sixth equation and is assumed equal to the fertile population (seventh equation); both of these are equal to the city population, ((2) and following paragraph). Finally, total infected population 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 .
State variables are next treated by categories and pools.
2.1.2. Susceptibles
The differential equation for is The two principal inputs are first from births with the number of fertile persons 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 pool is included, . This is the death rate caused by influenza (Figure 1). is given by (16). This is done so that the live person number, (1), remains constant at value (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 β term in (2)). Also, the instantaneous death rate from influenza, (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 , the virus quantity and is divided by the number of live persons (1). All terms in (2) scale equally with population size.
2.1.3. 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 from the inf1 pool is partly asymptomatic with fraction (Figure 1, Table 1; (7)), the remainder entering the inf2 pool (second equation, first term). The third and fourth equations are straightforward. Rate constant is 2 dayβ1, such that with four pools the mean time from infection to clinical sickness is 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 (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 (Table 1). The time between the onset of sickness and recovery is similarly (to the pools above) gamma-distributed (Appendix A, Figure 9, (A.4)).
The rate at which individuals are admitted to hospital, (Figure 1), is with units of persons . This is often a recorded statistic.
The total flu-related death flux is (input from all four sick pools to the dead-from-influenza box of Figure 1; units: persons )
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 represents hospital admissions (5), a recorded statistic, which may be useful for comparing with data. The flu-induced death rate (Table 1) is assumed to be the same as that for the nonhospitalized clinically sick, (4). This may be justified in that the fraction 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 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 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 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 (day) is the time period during which immunity is lost (three years is assumed; Table 1), and (dayβ1) is the rate constant out of each of the ten compartments, then The standard deviation in the exit times is and the coefficient of variation is 1/β10 () (A.7). The differential equations for the pools are The first term on the right of the 1st equation () is the last term in (10); the second term () is transfer to the next compartment; the last term (ΞΌ) is the natural death rate. The 2nd term on the right of the last equation () is transfer into the susceptible pool (2).
2.1.9. Virus Pool
It is assumed that virus production, , which is a surrogate for βinfectious strengthβ, is calculated with 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 () 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.β of (13) is the input to the virus pool (15).
The outputs from the virus pool (15) are death by natural mortality (, Table 1) and death induced by a suboptimal environment, represented by rate constant (18). The value of 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 . The differential equation for the virus pool 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 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.
2.1.11. Environmentally induced Virus Mortality and Seasonal Forcing
The environmentally dependent function in (15) (Figure 1) is assumed here to depend only on daily mean values of air temperature, . 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, (; Figure 1) is written as is used to switch environmental effects off (0) or on (1). Air temperature above a threshold increases virus mortality according to power . Parameter is assumed equal to two giving a quadratic dependence of temperature above the threshold (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). 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 varies from c. 3 (January) to 17Β°C (July) ([44], p. 270; [49] p. 142). It is assumed that varies sinusoidally, with Annual mean and seasonal variation are and . is the Julian day number (1 on 1 January, leap years are ignored). (d) is the phase of the sinusoid, which is maximum on 25 July. Combining equations (19) with (18) modifies environmentally induced virus death rate, , and hence reproductive ratio (20). The same seasonal pattern is applied every year. These equations for temperature are a good approximation to long-term weather means in the UK (Note: the study reported in [34] suggests that, on an annual timescale, cold weather does not predict winter deaths, but, on a shorter term timescale, cold weather could be a significant trigger).
2.2. Basic Reproductive Ratio, R0, and the Disease Generation Time, (Day)
Two important epidemiological parameters are basic reproductive ratio, , and the disease generation time, (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. is the average time it takes the direct infections which contribute to to arise. 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 ) 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 to and back to (Figure 1). The first term is Start with a single virus particle in the 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 . It is assumed that there is no additional virus death process due to exposure. is a dimensionless weighting factor (14). The third factor () is the number of virus particles produced per person in the inf1 state during his life: is the average lifetime of a person in the inf1 state. The last factor, with in the denominator, is the probability that this lifetime (of ) 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 box and diverting the (primary) infected individuals into an accumulator, rather than allowing them to enter the box, where they can lead to secondary infections. The algebraic and numerical methods agree to six decimal digits, giving a basic reproductive ratio for the default parameters without forcing.
A βdynamicβ reproductive ratio, , can be calculated when the total population is not entirely susceptible by means of The live population, , 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 (see previous section).
The disease generation time, (day), can be computed analytically by (B.7), giving day. A numerical calculation using Runge-Kutta integration gives also day agreeing with the analytical result to nine decimal digits, although Euler integration gives day (in each case the integration interval is 1/32 day; see first paragraph of Section 3). Due to overlapping generations, is not equal to the time constant at which the total infecteds increase. Assume that total infecteds (1) increase initially at an exponential rate according to where (day) is the growth time constant (the proportional growth rate of infecteds is dayβ1). This depends on both and 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 (day); (b) the latent period, namely the time which elapses between the infection event and the appearance of clinical sickness, denoted by (day); (c) the infectious period (day), which is the time during which the subject is infectious (Figure 1βproducing virus); (d) the period of immunity (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, 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, (day). With the infectivity weighting factors in (14), the first infected pool with a mean lifetime of 0.5 day is not infectious but the second infected pool is, and therefore Next the latent period (day) is given by (there are four sequential pools in the infected category of Figure 1) The infectious period (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 , and are all equal (Table 1), say to ; (c) the first infected pool () is assumed not infectious and all sick and hospitalized pools are infectious (14). We therefore write the infectious period as Loss of immunity occurs during a time period of , assumed to be 3 years. This delay is represented by (Figure 1) 10 sequential pools each having an outgoing rate constant of (dayβ1) (ignoring natural mortality (12)). This gives a gamma-distributed delay (Appendix A, (A.5) to (A.8) with ). and are related by Next compare the values in (24)β(27) with the literature. Our value of 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 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 = 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β 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 , 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, (persons dayβ1); Figure 1), as this is often a recorded statistic. An βepidemicβ is deemed to have occurred if shows a maximum (in steps of Ξt) and if at the maximum persons .
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)
(a)
(b)
(c)
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
3.2. Dynamics without Intrinsic Loss of Immunity and without Forcing
Here the model is exercised with parameter ((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 () peaks a few days before hospitalized numbers (), both then falling to very low values. Most (99%) of the population joins the recovered (immune) box (), with 1% escaping infection altogether. 90% of those infected travel via the asymptomatic route (Figure 1; , 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 (22) (not plotted) falling to near zero, as the fraction of susceptibles () becomes small. In the first epidemic, at time day, the initial peak in hospital admissions (Figure 1, (5)), decreasing to less than 100 at the second epidemic. After the first epidemic, (22) slowly recovers as immune numbers () decrease and susceptibles () 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 (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 (with ; Figure 1; Table 1; (3)), and the initial (time zero) immune fraction, . 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 ; 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 (e.g., Figure 2(b)) becomes closer together.
Figure 3(b), where the rates of transit of infected persons through the system are increased (, Figure 1), is less straightforward. Here a low value for the (e.g., 0.25 dayβ1) 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 spikes as in Figure 2(b) are closer together). As the are increased basic reproductive ratio decreases (20) first the epidemic becomes narrower and faster (e.g., dayβ1), but further increases in (e.g., and 8 dayβ1) cause the epidemic to become slower to take hold and less peaked, with fewer hospital admissions (HA). Increasing always causes the 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 (,(1)).
Finally, Figure 3(c) illustrates how increasing the initial immune fraction 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β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 (; 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 dayβ1 (instead of 0) if the variable of Figure 2(a) is replaced by the 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 , Table 1, (1) and following paragraph), with an inverse switching of numbers in the delay compartment, (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 = 1 (22), and the fractions in the three principal categories (Figure 1) of susceptible; all infected categories lumped together (infected, asymptomatic, sick, hospitalized, (1)) and the delayed categories are 0.209, 0.0030, and 0.787. Comparing these values with the situation in which there is no intrinsic loss of immunity (Figure 2(b)), it is seen that hospitalizations have increased by a factor of 20 as has the fraction in all infected categories, although the fractions of susceptibles and recovered or immunes have barely changed.
3.4. Validation
βValidationβ is often a misused and misunderstood concept and is perhaps better described as an evaluation of applicability. Validity is not a property of a model alone, neither is it a βzero or oneβ concept. It describes the relationship between model predictions and a set of data obtained under prescribed conditions. In this section, the model of Figure 1 is βvalidatedβ by fitting the predictions of the model, with minimal parameter adjustment (a perfectly formulated mechanistic model would permit no adjustment of parameters or initial values) to data on an influenza epidemic [53]. Success in this endeavour gives the model some credibility, although it does not make the model valid for general use.
The data in [53] relate to an influenza outbreak at an English boarding school, which provides a simple situation that seems comparable to the single-city homogeneous-mixing model of Figure 1, remembering that the model does not have the stratification which might be needed if a larger region was considered.
Minimal tuning is applied. One infected person is introduced to the school at noon on 18 January, otherwise all are susceptible. Table 1 parameters are altered for the Figure 5 simulation as follows: infection parameter Ξ² is changed to 0.07 (virus units)β1 dayβ1; all birth and death rates are set to zero () for such young persons; total population ; and the asymptomatic fraction, , 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, , is 1.015 dayβ1.
For comparison with the data in [53], we define the number of persons confined to bed, , as Since this definition is somewhat arbitrary, in comparing the predicted from the model with the data from Anon (1978), 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 pools for the E category and pools for the I category. They find a best fit with and . 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 (), relative humidity (), radiation (possibly multicomponent), and wind speed. The effects can be highly complex: for example, [48] which examines the effects of and 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 alone is operative.
Figure 6(a) illustrates the seasonal variation of mean daily air temperature in the southern UK. Using (19) with (18), and , this affects environmentally induced virus death rate (Figure 1, (15)), and thereby basic reproductive ratio (20). Mean daily air temperature varies between 3Β°C (Jan 24) and 17Β°C (Jul 25). When crosses the temperature threshold of = 13Β°C in May, this increases environmentally induced virus mortality, (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, (19), (or equivalently, decreasing threshold temperature , (18)) varies the duration and intensity of seasonal forcing of . Changing mean annual air temperature, (19), changes the average annual value of , , as well as its maximum, minimum and amplitude (maximum-minimum). The dependence of these quantities on is also shown, together with the number of forcing days per year (, a forcing day is one on which mean air temperature is above threshold temperature and virus mortality is reduced). With mean annual air temperature , there is no forcing at all (): air temperature (19) never exceeds threshold temperature = 13Β°C (18); is invariant at its maximum value. Environmentally induced virus death rate ((15), (18)) stays at zero; the system remains in the steady state shown in Figure 4(b) and (this situation is equivalent to and = 17Β°C). As mean annual air temperature increases above 6Β°C, there are more days in the summer months when mean daily temperature > threshold temperature , number of forcing days per year increases, and . The amplitude increases to a maximum when, before decreasing. With , 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 are taken between β4.5 and 16Β°C, assuming always an annual variation 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 (; 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, (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 ) 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 (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); ), 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); ), there is chaos but now with a tendency towards biennial epidemics (169 epidemics occur during 200 years). Last (Figure 8(d); ), 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 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 can be continued, and although the situations simulated are now less realistic, they do contribute to understanding the system. With , 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 , 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 the response is chaotic with 15 epidemics in the first 200 years and . Last, with , 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 . This method is used in many SIR models. A single pool emptying at rate from initial value of at time obeys the equations 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 (shown by β). The mean time for departure is at day = ( (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, . In this case, the relevant equations become This is shown in Figure 9: the line labelled β2 pools, β. With two pools, has been doubled relative to in Eqns (A.1), so that the mean time for departure from the final pool in the category is the same at time day = (shown by β ). We see that the maximum rate for departure from the infected category now occurs at time day (= ()/) (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 . The equations for the system are Now (Figure 9) the curve for the final pool in the category () is sigmoidal at low values of time . The time of the maximum (β) (= (()/) has moved closer to the mean time (β ). The mean time is ()/() = 4 day as in the other curves. With three pools, there is now a sigmoid departure from time , 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 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 .
For a general number of pools, , with state variables, , and each with outgoing rate constant, , the value of the state variable for the final pool in the sequence, , is where is a constant. The total outflow from the final pool, , is , so that The statistics on the final pool, , are (using to denote expectation values) 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β1, giving an overall average transit time of 4 day (see main text for further discussion of this point).
Finally we note that although there are several equivalent definitions for the gamma function used by mathematicians and others (e.g., [44], pp. 819β821), possibly the most intuitive definition for the biologist is that given in (A.6), namely, where the normalized gamma function, , is given by 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
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, next for the four sick pools then, the three hospitalized pools and last the seven asymptomatic pools Basic reproductive ratio is the sum of these contributions: Generation time (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:Finally, summing the above, the mean generation time is
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.