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

Authors: John H. M. Thornley; James France
Journal: ISRN Biomathematics
 (2012)

---

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


---

## Body

## 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 (Figure1) so that the underlying biological dynamics is realistic, to validate the model by applying it to data describing an influenza epidemic (Figure 5), to suggest a possible mechanism for the direct effects of environment (season) on influenza dynamics (Figure 6), and last to apply various levels of environmental forcing to demonstrate how the model, without reparameterization, can give rise to a wide range of dynamics, ranging from regular (e.g., twice yearly, annual, biennial, etc.) epidemics to chaos, where sometimes predictions resemble some of the great influenza pandemics (Figures 7 and 8).Figure 1
Model scheme. The 32 state variables are denoted byn + subscript. All pools other than the virus pool (nvir) and the dead-from-influenza pool (nded) 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 kinf. Similarly, ksic applies to the four nonhospitalized clinically sick pools, khos to the three hospitalized clinically sick pools, kasy to the seven asymptomatic pools, and kdel to the ten pools in the loss-of-immunity delay pipe. fasy and fhos are the fractions of the fluxes between ninf1 and ninf2, and nsic1 and nsic2, which are asymptomatic and hospitalized. μsic and μhos are the death rates (caused by flu) of the clinically sick and hospitalized pools.Before proceeding with this agenda, a brief review of some of the issues and modelling approaches relevant to influenza is given, and our contribution is discussed in this wider context.


### 1.2. Brief Review of Influenza Issues and Modelling Approaches

#### 1.2.1. Avian Influenza

One of the reasons for the ongoing concern with influenza is the continuing threat from avian flu. Avian influenza is endemic in many parts of the world and, from time to time, there are outbreaks of highly pathogenic strains, such as the recent outbreak of the H5N1 strain of the virus [8]. The virus can be transferred to other animals, including horses, pigs, and humans, frequently with fatal consequences [9]. Because of antigenic drift (a high rate of immunologically significant mutations) and shift (reassortment between different strains of influenza within a single host), coupled to other factors such as loss of immunity, there is a continuing threat to the human population [10, 11]. This threat depends on possible changes in the to-date limited capacity of the highly pathogenic avian strains for human-to-human transmission [8, 9]. This risk makes it important that key factors determining virus evolution, epidemic occurrence, and cross-species transmission are well understood, so that effective strategies for containment and control might be designed [12–14].There are also models dealing with the avian-human influenza nexus. For example, [14] presents an ordinary differential equation model which combines an SI (susceptible, infected) avian model with an SIR (susceptible, infected, resistant) human model. They suggest that measures such as both extermination (of avians) and quarantine (for humans) could be needed to avoid a pandemic of influenza.


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

An interesting paper [15] extends what might be called the complex models of simple influenza by presenting a simple deterministic model of a more complex situation: they treat the dynamics of two-strain influenza, focussing on competition and cross-immunity. Isolation period and crossimmunity are critical parameters. Some of their results are similar to those reported here (with variable interepidemic periods from 2 to 10–13 years), although the model and mechanisms are quite different.At a more detailed level, the cell level in a single person, a model of the immune response to the influenza virus which treats innate and adaptive immunity has been proposed [16]. The model has 10 ordinary differential equations, representing interferon, T-cells, killer cells, antibodies, and other states. They explore the impact of initial viral load on disease progression. When this is small, the disease progresses asymptomatically. The model builds on [17].


#### 1.2.3. Network Models and Stochastic Simulations

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


#### 1.2.4. Deterministic Models

Deterministic models of influenza are also numerous, and the model presented here is of this type. While biological variability is a reality, there is a feeling amongst some modellers, especially those with a background in engineering or physics or applied mathematics, that a deterministic approach gets closer to the real science, to understand what is going on, than a stochastic approach. Reflecting this are many quotations, such as “God does not play dice,” and “if you need statistics, then do a better experiment.” Deterministic models are easier to build, easier to understand and use, easier to falsify, and easier for tracing cause and effect. With a deterministic model, a clearly labelled box diagram tells the reader most of what he needs to know. In the next two paragraphs, we refer to two examples of the genre. Although many models have the potential to represent directly seasonal effects of environment, we have not found one that actually does this. Some models represent places where people assemble and interact such as schools and clinics which affect contact rates: indirectly, these can represent a seasonality. But here, effects of weather variables on model parameters are represented directly.The model of [25] is a typical example of a simple deterministic influenza model. The model has seven significant state variables and ordinary differential equations. From our point of view, their use of a single stage to represent each category (e.g., latent, infectious, asymptomatic, hospitalized, etc.) gives a biologically unreasonable representation of the transit time distribution through that category. Since here we represent day-to-day effects of weather on influenza dynamics, it is important that the model should represent the dynamics of influenza realistically. This means using several serial stages for progress through each disease category.A typicallarge 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.


## 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 (Figure1) so that the underlying biological dynamics is realistic, to validate the model by applying it to data describing an influenza epidemic (Figure 5), to suggest a possible mechanism for the direct effects of environment (season) on influenza dynamics (Figure 6), and last to apply various levels of environmental forcing to demonstrate how the model, without reparameterization, can give rise to a wide range of dynamics, ranging from regular (e.g., twice yearly, annual, biennial, etc.) epidemics to chaos, where sometimes predictions resemble some of the great influenza pandemics (Figures 7 and 8).Figure 1
Model scheme. The 32 state variables are denoted byn + subscript. All pools other than the virus pool (nvir) and the dead-from-influenza pool (nded) 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 kinf. Similarly, ksic applies to the four nonhospitalized clinically sick pools, khos to the three hospitalized clinically sick pools, kasy to the seven asymptomatic pools, and kdel to the ten pools in the loss-of-immunity delay pipe. fasy and fhos are the fractions of the fluxes between ninf1 and ninf2, and nsic1 and nsic2, which are asymptomatic and hospitalized. μsic and μhos are the death rates (caused by flu) of the clinically sick and hospitalized pools.Before proceeding with this agenda, a brief review of some of the issues and modelling approaches relevant to influenza is given, and our contribution is discussed in this wider context.


## 1.2. Brief Review of Influenza Issues and Modelling Approaches

### 1.2.1. Avian Influenza

One of the reasons for the ongoing concern with influenza is the continuing threat from avian flu. Avian influenza is endemic in many parts of the world and, from time to time, there are outbreaks of highly pathogenic strains, such as the recent outbreak of the H5N1 strain of the virus [8]. The virus can be transferred to other animals, including horses, pigs, and humans, frequently with fatal consequences [9]. Because of antigenic drift (a high rate of immunologically significant mutations) and shift (reassortment between different strains of influenza within a single host), coupled to other factors such as loss of immunity, there is a continuing threat to the human population [10, 11]. This threat depends on possible changes in the to-date limited capacity of the highly pathogenic avian strains for human-to-human transmission [8, 9]. This risk makes it important that key factors determining virus evolution, epidemic occurrence, and cross-species transmission are well understood, so that effective strategies for containment and control might be designed [12–14].There are also models dealing with the avian-human influenza nexus. For example, [14] presents an ordinary differential equation model which combines an SI (susceptible, infected) avian model with an SIR (susceptible, infected, resistant) human model. They suggest that measures such as both extermination (of avians) and quarantine (for humans) could be needed to avoid a pandemic of influenza.


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

An interesting paper [15] extends what might be called the complex models of simple influenza by presenting a simple deterministic model of a more complex situation: they treat the dynamics of two-strain influenza, focussing on competition and cross-immunity. Isolation period and crossimmunity are critical parameters. Some of their results are similar to those reported here (with variable interepidemic periods from 2 to 10–13 years), although the model and mechanisms are quite different.At a more detailed level, the cell level in a single person, a model of the immune response to the influenza virus which treats innate and adaptive immunity has been proposed [16]. The model has 10 ordinary differential equations, representing interferon, T-cells, killer cells, antibodies, and other states. They explore the impact of initial viral load on disease progression. When this is small, the disease progresses asymptomatically. The model builds on [17].


### 1.2.3. Network Models and Stochastic Simulations

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


### 1.2.4. Deterministic Models

Deterministic models of influenza are also numerous, and the model presented here is of this type. While biological variability is a reality, there is a feeling amongst some modellers, especially those with a background in engineering or physics or applied mathematics, that a deterministic approach gets closer to the real science, to understand what is going on, than a stochastic approach. Reflecting this are many quotations, such as “God does not play dice,” and “if you need statistics, then do a better experiment.” Deterministic models are easier to build, easier to understand and use, easier to falsify, and easier for tracing cause and effect. With a deterministic model, a clearly labelled box diagram tells the reader most of what he needs to know. In the next two paragraphs, we refer to two examples of the genre. Although many models have the potential to represent directly seasonal effects of environment, we have not found one that actually does this. Some models represent places where people assemble and interact such as schools and clinics which affect contact rates: indirectly, these can represent a seasonality. But here, effects of weather variables on model parameters are represented directly.The model of [25] is a typical example of a simple deterministic influenza model. The model has seven significant state variables and ordinary differential equations. From our point of view, their use of a single stage to represent each category (e.g., latent, infectious, asymptomatic, hospitalized, etc.) gives a biologically unreasonable representation of the transit time distribution through that category. Since here we represent day-to-day effects of weather on influenza dynamics, it is important that the model should represent the dynamics of influenza realistically. This means using several serial stages for progress through each disease category.A typicallarge 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.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 typicallarge 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 Figure1. 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 (infj,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 ndel1,ndel2,…, to ndel10. Infectious persons give rise to virus particles (nvir), 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 (nvir). 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 (nvir) would be to assume that weather affects transmission rate β directly. In the equations to follow, transmission rate β and virus pool nvir 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.Table 1
State variables, parameters and definitions. Initial values of state variables, principal parameter values and definitions of principal variables. See Figure1 for state variables. Relevant equation numbers are given.
State variable initial values (t = 0)Unitsnsus=ncity-ninf1-nrecSusceptiblesPersonsninf1=1InfectedsPersonsnrec=0 usually, but see Figure 3(c)Recovered (immune)PersonsAll othernxxx,nxxxx=0 (Figure 1)Persons, virusesParameters + valueDescriptionUnitsfasy=0.9Asymptomatic fraction (3)fhos=0.05Fraction of sick hospitalized (4)kasy=2Transfer rate (7)Day−1kdel=10/1095 (see τimm below)Transfer rate (11)Day−1khos=2Transfer rate (8)Day−1kinf=2Transfer rate (3)Day−1krec=0(or1)Rate constant (10)Day−1ksic=2Transfer rate (4)Day−1kvir=200Virus production parameter (13)Virus units person−1 day−1kvir,dnm=5Virus natural mortality (15)Day−1ncity=106Population of city (see 1st line of table, paragraph before (1), and paragraph after (2))PersonsR0= 4.78 by (20)Basic reproductive ratioβ = 0.1Infection transmission rate (2)Persons (virus units)−1 day−1μ=1/(60×365)=4.566×10-5Natural birth/death rate (2)Day−1μhos=0.001Flu-induced death rate (hospital) (8)Day−1μsic=0.001Flu-induced death rate (home) (4)Day−1τgen=2.51 ((B.7); see 2nd paragraph after (22))Generation timeDay−1τimm=3×365=1095Loss of immunity delay (11)DayEnvironmental parametersDescriptionUnitsqT=2Temperature exponent (18)senv=0 (1)Environment switch (18)Tair0=13Temperature threshold (18)°CkTair=0.1Scaling factor (18)(°C)−2d-1VariablesDescriptionUnitskhosHospital admission rate (5)Persons day−1kvir,envEnvironmentally induced death rate ((18), (15))Day−1R0Basic reproductive ratio (20)sinfInfectious strength (virus production rate) (13)Virus person−1 day−1Initial values and parameters are listed in Table1, although some parameters which are only used once are defined in or after the equation where they appear.All routes from susceptible to recovered (Figure1) pass through eight pools. We have assigned the same value of 2 day−1 to the four rate constants kinf,kasy,khos, and ksic (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 (Figure1) 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 (ncity, Table 1).


#### 2.1.1. Some Definitions

First, define some totals in terms of the state variables (Figure1):(1)ninf⁡=ninf1+ninf2+ninf3+ninf4,nsic=nsic1+nsic2+nsic3+nsic4,nhos=nhos2+nhos3+nhos4,nasy=nasy2+nasy3+⋯+nasy7+nasy8,ndel=ndel1+ndel2+⋯+ndel9+ndel10,nliv=nsus+ninf⁡+nsic+nhos+nasy+nrec+ndel,nfer=nliv=ncity,niash=ninf⁡+nasy+nsic+nhos.
The notation of the first five equations is obvious. Total live population nliv is given by the sixth equation and is assumed equal to the fertile population nfer (seventh equation); both of these are equal to the city population, ncity ((2) and following paragraph). Finally, total infected population niash 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 ncity.State variables are next treated by categories and pools.


#### 2.1.2. Susceptibles

The differential equation fornsus is(2)dnsus(t)dt=μnfer+kdelndel10+Ided-μnsus-βnsusnvirnliv.
The two principal inputs are first from births with the number of fertile persons nfer 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 nsus pool is included, Ided. This is the death rate caused by influenza (Figure 1). Ided is given by (16). This is done so that the live person number, nliv (1), remains constant at value ncity (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 −μnsus term in (2)). Also, the instantaneous death rate from influenza, Ided (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 nsus, the virus quantity nvir and is divided by the number of live persons nliv (1). All terms in (2) scale equally with population size.


#### 2.1.3. Infecteds

The differential equations for the four sequential infected pools are(3)dninf1dt=βnsusnvirnliv-kinf⁡ninf1-μninf1,dninf2dt=(1-fasy)kinf⁡ninf1-kinf⁡ninf2-μninf2,dninf3dt=kinf⁡ninf2-kinf⁡ninf3-μninf3,dninf4dt=kinf⁡ninf3-kinf⁡ninf4-μninf4.
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 kinfninf1 from the inf1 pool is partly asymptomatic with fraction fasy (Figure 1, Table 1; (7)), the remainder entering the inf2 pool (second equation, first term). The third and fourth equations are straightforward. Rate constant kinf is 2 day−1, such that with four pools the mean time from infection to clinical sickness is 4/kinf=2 days [45]. This gives a gamma-distributed lag for the overall exit time from the fourth pool (Appendix A, (A.4), Figure 9.)


#### 2.1.4. Nonhospitalized Clinically Sick

Persons who become clinically sick enter the sic1 pool (Figure1). A fraction fhos (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(4)dnsic1dt=kinf⁡ninf4-ksicnsic1-(μ+μsic)nsic1,dnsic2dt=(1-fhos)ksicnsic1-ksicnsic2-(μ+μsic)nsic2,dnsic3dt=ksicnsic2-ksicnsic3-(μ+μsic)nsic3,dnsic4dt=ksicnsic3-ksicnsic4-(μ+μsic)nsic4.
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 ninf pools above) gamma-distributed (Appendix A, Figure 9, (A.4)).The rate at which individuals are admitted to hospital,Ihos2 (Figure 1), is(5)Ihos2=fhosksicnsic1,
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-influenzanded box of Figure 1; units: persons day-1)(6)Isic→ded=μsic(nsic1+nsic2+nsic3+nsic4).


#### 2.1.5. Asymptomatic Infecteds

The differential equations for these seven pools (Figure1) are(7)dnasy2dt=fasykinf⁡ninf1-kasynasy2-μnasy2,dnasy3dt=kasynasy2-kasynasy3-μnasy3,…dnasy8dt=kasynasy7-kasynasy8-μnasy8.
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(8)dnhos2dt=fhosksicnsic1-khosnhos2-(μ+μhos)nhos2,dnhos3dt=khosnhos2-khosnhos3-(μ+μhos)nhos3,dnhos4dt=khosnhos3-khosnhos4-(μ+μhos)nhos4.
The input term fhosksicnsic1 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 fhos 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(9)Ihos→ded=μhos(nhos2+nhos3+nhos4).


#### 2.1.7. Recovered and Immune

The state variablenrec is governed by the equation(10)dnrecdt=kasynasy8+ksicnsic4+khosnhos4-μnrec-krecnrec.
The positive (input) terms are from the last equations of (7), (4), and (8). Death occurs at its natural rate (μ). The rate constant krec 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 krec causes rapid entry into the delay sequence of pools (Figure 1) which leads to loss of immunity.


#### 2.1.8. Delay Pipe Representing Loss of Immunity

Loss of immunity might arise from loss of immunological memory, or from drift and shift in the antigenic character of the virus [11]. This process is represented by ten sequential compartments giving a gamma function delay (Appendix A; also e.g., [44], pp. 818–822). If τimm (day) is the time period during which immunity is lost (three years is assumed; Table 1), and kdel (day−1) is the rate constant out of each of the ten compartments, then(11)kdel=10τimm.
The standard deviation in the exit times is (√10)/kdel and the coefficient of variation is 1/√10 (=standard deviation/τimm) (A.7). The differential equations for the pools are(12)dndel1dt=krecnrec-kdelndel1-μndel1,dndel2dt=kdelndel1-kdelndel2-μndel2,…dndel9dt=kdelndel8-kdelndel9-μndel9,dndel10dt=kdelndel9-kdelndel10-μndel10;
The first term on the right of the 1st equation (krec) is the last term in (10); the second term (kdel) is transfer to the next compartment; the last term (μ) is the natural death rate. The 2nd term on the right of the last equation (kdel) is transfer into the susceptible pool (2).


#### 2.1.9. Virus Pool

It is assumed that virus production,svir, which is a surrogate for “infectious strength”, is calculated with(13)svir(t)=kvir(winf1ninf1+winf2ninf2+winf3ninf3+winf4ninf⁡4+wsic1nsic1+wsic2nsic2+wsic3nsic3+wsic4nsic4+wasy2nasy2+wasy3nasy3+…+wasy8nasy8+whos2nhos2+whos3nhos3+whos4nhos4).kvir=200virusunits(weightedinfectedperson)-1day-1.kvir is a production rate constant which is applied to weighted infected persons. The value above gives a basic reproductive ratio R0 (20) of 4.78 (see Section 2.3). The weighting factors w of (13) are (0≤w≤1)(14)winf⁡1=0,winf2=0.5,winf3=1,winf4=1;wsic1=1,wsic2=1,wsic3=1,wsic4=0.5;wasy2=0.2,wasy3=0.4,wasy4=0.4,wasy5=0.4,wasy6=0.4,wasy7=0.2,wasy8=0;whos2=0.2,whos3=0.2,whos4=0.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.” svir of (13) is the input to the virus pool (15).The outputs from the virus pool (15) are death by natural mortality (kvir,dnm, Table 1) and death induced by a suboptimal environment, represented by rate constant kvir,env (18). The value of kvir,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 kvir,dnm. The differential equation for the virus pool nvir is(15)dnvirdt=svir-(kvir,dnm+kvir,env)nvir.


#### 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 nded pool (Figure 1) is(16)Ided=Isic→ded+Ihos→ded.
There are no outputs. Therefore,(17)dndeddt=Ided,t=0,nded=0.
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 Mortalitykvir,env and Seasonal Forcing

The environmentally dependent functionkvir,env in (15) (Figure 1) is assumed here to depend only on daily mean values of air temperature, Tair. 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, kvir,eny (d-1; Figure 1) is written as(18)kvir,env=senvkTair[(Tair-Tair0)+ABS(Tair-Tair0)2]qT,senv=0(or1),kTair=0.1(C°)-2d-1,Tair0=13C°,qT=2.senv is used to switch environmental effects off (0) or on (1). Air temperature Tair above a threshold Tair0 increases virus mortality according to power qT. Parameter qT is assumed equal to two giving a quadratic dependence of temperature above the threshold Tair0 (Table 1). It is usual for the biological effects of temperature to be nonlinear, sometimes approximating to exponential, as in the use of a Q10 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). kTair 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 temperatureTair varies from c. 3 (January) to 17°C (July) ([44], p. 270; [49] p. 142). It is assumed that Tair varies sinusoidally, with(19)Tair=Tair,mean+Tair,var{sin⁡[2π365(iJulian-τTair)]};Tair,mean=10°C(varied),Tair,var=7°C,τTair=115day.
Annual mean and seasonal variation are Tair,mean and Tair,vr. iJulian is the Julian day number (1 on 1 January, leap years are ignored). τTair (d) is the phase of the sinusoid, which is maximum on 25 July. Combining equations (19) with (18) modifies environmentally induced virus death rate, kvir,env, and hence reproductive ratio (20). The same seasonal pattern is applied every year. These equations for temperature are a good approximation to long-term weather means in the UK (Note: the study reported in [34] suggests that, on an annual timescale, cold weather does not predict winter deaths, but, on a shorter term timescale, cold weather could be a significant trigger).


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

Two important epidemiological parameters are basic reproductive ratio,R0, and the disease generation time, τgen (day). These are both derived from the basic parameters of the model. R0 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 R0 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 caseR0 can only be calculated numerically), R0 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 R0) (20)R0=βkvir,nm+kvir,envkvirkinf⁡kinf⁡+μ×{(kasykasy+μ)6(khoskhos+μ+μhos)3winf1kinf⁡+(1-fasy)×〈(khoskhos+μ+μhos)3kinf⁡kinf⁡+μwinf2kinf⁡+(kinf⁡kinf⁡+μ)2winf3kinf⁡+(kinf⁡kinf⁡+μ)3×⟦(khoskhos+μ+μhos)3winf4kinf⁡+ksicksic+μ+μsic×{(khoskhos+μ+μhos)3(kasykasy+μ)6wsic1ksic+(1-fhos)[ksicksic+μ+μsicwsic2ksic+(ksicksic+μ+μsic)2wsic3ksic+(ksicksic+μ+μsic)3wsic4ksic]+fhos[khoskhos+μ+μhoswhos2khos+(khoskhos+μ+μhos)2whos3khos+(khoskhos+μ+μhos)3whos4khos]}⟧〉+fasykasykasy+μ[wasy2kasy+(kasykasy+μ)wasy3kasy+(kasykasy+μ)2wasy4+⋯+(kasykasy+μ)6wasy8kasy](khoskhos+μ+μhos)3}.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 nvir to ninf1 and back to nvir (Figure 1). The first term is(21)βkvir,nm+kvir,envwinf⁡1kvirkinf⁡kinf⁡kinf⁡+μ.
Start with a single virus particle in the nvir box (we could equally well start with a single person in the inf1 first-infected compartment). The first factor of (21), with units of persons per virus particle, is the probable number of persons infected by a virus particle during its life; the mean lifetime of a virus particle is 1/(kvir,dnm+kvir,env). It is assumed that there is no additional virus death process due to exposure. winf1 is a dimensionless weighting factor (14). The third factor (kvir/kinf) is the number of virus particles produced per person in the inf1 state during his life: 1/kinf is the average lifetime of a person in the inf1 state. The last factor, with kinf+μ in the denominator, is the probability that this lifetime (of 1/kinf) 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 R0 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 ninf1 box and diverting the (primary) infected individuals into an accumulator, rather than allowing them to enter the ninf1 box, where they can lead to secondary infections. The algebraic and numerical methods agree to six decimal digits, giving a basic reproductive ratio R0=4.78 for the default parameters without forcing.A “dynamic” reproductive ratio,Rdyn, can be calculated when the total population is not entirely susceptible by means of(22)Rdyn=R0nsusnliv.
The live population, nliv, is given by (1). This allows approximately for a slowly changing fraction of susceptible individuals. However, R0 may itself be dynamic on a shorter time scale due to seasonal effects on virus death rate kvir,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 niash (1) increase initially at an exponential rate according to(23)niash~et/τ,
where τ (day) is the growth time constant (the proportional growth rate of infecteds is 1/τ day−1). This depends on both R0 and τgen and the structure of the model [6, Section  1.2.3]. This can be extracted numerically (there is a period of constant exponential growth that lasts for some 10 days) to give τ = 1.29 day. This is considerably less than the generation time of 2.51 day.


### 2.3. Parameterization

Parameters have been introduced while developing the model. Their values are listed for reference in Table1. Here we summarize the evidential basis for the parameter values used.As a preliminary, statements from two physicians are quoted. In [50] a general practitioner in the Doncaster (UK) area describes his study of the 1969-1970 pandemic as it affected his urban practice. He said “the true incidence of influenza during an epidemic is probably impossible to assess.” This is perhaps equally true today and sets the scene for the significance of processes such as “validation” and data fitting (see Figure 5 and Section 3.4). Another general practitioner, this time in Kent (UK) [51], states that “an epidemic of influenza tends to last in this area for between two and three months. Beginning slowly the epidemic reaches its peak in four to five weeks and then subsides slowly. The extent and severity of any attack will depend on factors such as the strain of influenza virus, on the state of the host-immunity of the population and on the timing of the epidemic; the fatality and complication rates are always higher during the cold and foggy winter months.” His statement agrees with many of our seasonal-forcing simulations (Figures 6 and 7).Five key epidemiological quantities for influenza are (a) the time which elapses after the infection event until the subject becomes infectious to others, denoted by theτinf0 (day); (b) the latent period, namely the time which elapses between the infection event and the appearance of clinical sickness, denoted by τlat (day); (c) the infectious period τinf (day), which is the time during which the subject is infectious (Figure 1—producing virus); (d) the period of immunity τimm (day)—the time period after recovery from influenza during which the subject has immunity, before gradually losing it and returning to the susceptible pool (Figure 1, τimm is represented by the box at the bottom of the diagram with 10 sequential pools); (e) the basic reproductive ratio R0 (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 ninf1 with a mean lifetime of 0.5 day is not infectious but the second infected pool ninf2 is, and therefore(24)τinf⁡0=1(kinf⁡=2day-1)=0.5day.
Next the latent period τlat (day) is given by (there are four sequential pools in the infected category of Figure 1)(25)τlat=4(kinf⁡=2day-1)=2day.
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 kinf,kasy,ksic, and khos are all equal (Table 1), say to kdis; (c) the first infected pool (ninf1) is assumed not infectious and all sick and hospitalized pools are infectious (14). We therefore write the infectious period τinf as(26)τinf⁡=7(kdis=2day-1)=3.5day.
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 kdel (day−1) (ignoring natural mortality (12)). This gives a gamma-distributed delay (Appendix A, (A.5) to (A.8) with m=10). τimm and kdel are related by(27)kdel=10τimmwithτimm=1095day=3year.
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 winf2=wasy2=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 ratioR0 (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 R0 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 R0 = 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 R0 and therefore will underestimate the level of global control measures … that will be needed ….” [52] reviews R0 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 R0 = 4.78 (20). This can be regarded as an upper baseline applying to an optimum environment (one maximizing R0), 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 Figure1, 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.


## 2.1. Model Scheme

The scheme is drawn in Figure1. 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 (infj,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 ndel1,ndel2,…, to ndel10. Infectious persons give rise to virus particles (nvir), 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 (nvir). 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 (nvir) would be to assume that weather affects transmission rate β directly. In the equations to follow, transmission rate β and virus pool nvir 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.Table 1
State variables, parameters and definitions. Initial values of state variables, principal parameter values and definitions of principal variables. See Figure1 for state variables. Relevant equation numbers are given.
State variable initial values (t = 0)Unitsnsus=ncity-ninf1-nrecSusceptiblesPersonsninf1=1InfectedsPersonsnrec=0 usually, but see Figure 3(c)Recovered (immune)PersonsAll othernxxx,nxxxx=0 (Figure 1)Persons, virusesParameters + valueDescriptionUnitsfasy=0.9Asymptomatic fraction (3)fhos=0.05Fraction of sick hospitalized (4)kasy=2Transfer rate (7)Day−1kdel=10/1095 (see τimm below)Transfer rate (11)Day−1khos=2Transfer rate (8)Day−1kinf=2Transfer rate (3)Day−1krec=0(or1)Rate constant (10)Day−1ksic=2Transfer rate (4)Day−1kvir=200Virus production parameter (13)Virus units person−1 day−1kvir,dnm=5Virus natural mortality (15)Day−1ncity=106Population of city (see 1st line of table, paragraph before (1), and paragraph after (2))PersonsR0= 4.78 by (20)Basic reproductive ratioβ = 0.1Infection transmission rate (2)Persons (virus units)−1 day−1μ=1/(60×365)=4.566×10-5Natural birth/death rate (2)Day−1μhos=0.001Flu-induced death rate (hospital) (8)Day−1μsic=0.001Flu-induced death rate (home) (4)Day−1τgen=2.51 ((B.7); see 2nd paragraph after (22))Generation timeDay−1τimm=3×365=1095Loss of immunity delay (11)DayEnvironmental parametersDescriptionUnitsqT=2Temperature exponent (18)senv=0 (1)Environment switch (18)Tair0=13Temperature threshold (18)°CkTair=0.1Scaling factor (18)(°C)−2d-1VariablesDescriptionUnitskhosHospital admission rate (5)Persons day−1kvir,envEnvironmentally induced death rate ((18), (15))Day−1R0Basic reproductive ratio (20)sinfInfectious strength (virus production rate) (13)Virus person−1 day−1Initial values and parameters are listed in Table1, although some parameters which are only used once are defined in or after the equation where they appear.All routes from susceptible to recovered (Figure1) pass through eight pools. We have assigned the same value of 2 day−1 to the four rate constants kinf,kasy,khos, and ksic (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 (Figure1) 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 (ncity, Table 1).


### 2.1.1. Some Definitions

First, define some totals in terms of the state variables (Figure1):(1)ninf⁡=ninf1+ninf2+ninf3+ninf4,nsic=nsic1+nsic2+nsic3+nsic4,nhos=nhos2+nhos3+nhos4,nasy=nasy2+nasy3+⋯+nasy7+nasy8,ndel=ndel1+ndel2+⋯+ndel9+ndel10,nliv=nsus+ninf⁡+nsic+nhos+nasy+nrec+ndel,nfer=nliv=ncity,niash=ninf⁡+nasy+nsic+nhos.
The notation of the first five equations is obvious. Total live population nliv is given by the sixth equation and is assumed equal to the fertile population nfer (seventh equation); both of these are equal to the city population, ncity ((2) and following paragraph). Finally, total infected population niash 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 ncity.State variables are next treated by categories and pools.


### 2.1.2. Susceptibles

The differential equation fornsus is(2)dnsus(t)dt=μnfer+kdelndel10+Ided-μnsus-βnsusnvirnliv.
The two principal inputs are first from births with the number of fertile persons nfer 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 nsus pool is included, Ided. This is the death rate caused by influenza (Figure 1). Ided is given by (16). This is done so that the live person number, nliv (1), remains constant at value ncity (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 −μnsus term in (2)). Also, the instantaneous death rate from influenza, Ided (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 nsus, the virus quantity nvir and is divided by the number of live persons nliv (1). All terms in (2) scale equally with population size.


### 2.1.3. Infecteds

The differential equations for the four sequential infected pools are(3)dninf1dt=βnsusnvirnliv-kinf⁡ninf1-μninf1,dninf2dt=(1-fasy)kinf⁡ninf1-kinf⁡ninf2-μninf2,dninf3dt=kinf⁡ninf2-kinf⁡ninf3-μninf3,dninf4dt=kinf⁡ninf3-kinf⁡ninf4-μninf4.
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 kinfninf1 from the inf1 pool is partly asymptomatic with fraction fasy (Figure 1, Table 1; (7)), the remainder entering the inf2 pool (second equation, first term). The third and fourth equations are straightforward. Rate constant kinf is 2 day−1, such that with four pools the mean time from infection to clinical sickness is 4/kinf=2 days [45]. This gives a gamma-distributed lag for the overall exit time from the fourth pool (Appendix A, (A.4), Figure 9.)


### 2.1.4. Nonhospitalized Clinically Sick

Persons who become clinically sick enter the sic1 pool (Figure1). A fraction fhos (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(4)dnsic1dt=kinf⁡ninf4-ksicnsic1-(μ+μsic)nsic1,dnsic2dt=(1-fhos)ksicnsic1-ksicnsic2-(μ+μsic)nsic2,dnsic3dt=ksicnsic2-ksicnsic3-(μ+μsic)nsic3,dnsic4dt=ksicnsic3-ksicnsic4-(μ+μsic)nsic4.
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 ninf pools above) gamma-distributed (Appendix A, Figure 9, (A.4)).The rate at which individuals are admitted to hospital,Ihos2 (Figure 1), is(5)Ihos2=fhosksicnsic1,
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-influenzanded box of Figure 1; units: persons day-1)(6)Isic→ded=μsic(nsic1+nsic2+nsic3+nsic4).


### 2.1.5. Asymptomatic Infecteds

The differential equations for these seven pools (Figure1) are(7)dnasy2dt=fasykinf⁡ninf1-kasynasy2-μnasy2,dnasy3dt=kasynasy2-kasynasy3-μnasy3,…dnasy8dt=kasynasy7-kasynasy8-μnasy8.
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(8)dnhos2dt=fhosksicnsic1-khosnhos2-(μ+μhos)nhos2,dnhos3dt=khosnhos2-khosnhos3-(μ+μhos)nhos3,dnhos4dt=khosnhos3-khosnhos4-(μ+μhos)nhos4.
The input term fhosksicnsic1 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 fhos 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(9)Ihos→ded=μhos(nhos2+nhos3+nhos4).


### 2.1.7. Recovered and Immune

The state variablenrec is governed by the equation(10)dnrecdt=kasynasy8+ksicnsic4+khosnhos4-μnrec-krecnrec.
The positive (input) terms are from the last equations of (7), (4), and (8). Death occurs at its natural rate (μ). The rate constant krec 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 krec causes rapid entry into the delay sequence of pools (Figure 1) which leads to loss of immunity.


### 2.1.8. Delay Pipe Representing Loss of Immunity

Loss of immunity might arise from loss of immunological memory, or from drift and shift in the antigenic character of the virus [11]. This process is represented by ten sequential compartments giving a gamma function delay (Appendix A; also e.g., [44], pp. 818–822). If τimm (day) is the time period during which immunity is lost (three years is assumed; Table 1), and kdel (day−1) is the rate constant out of each of the ten compartments, then(11)kdel=10τimm.
The standard deviation in the exit times is (√10)/kdel and the coefficient of variation is 1/√10 (=standard deviation/τimm) (A.7). The differential equations for the pools are(12)dndel1dt=krecnrec-kdelndel1-μndel1,dndel2dt=kdelndel1-kdelndel2-μndel2,…dndel9dt=kdelndel8-kdelndel9-μndel9,dndel10dt=kdelndel9-kdelndel10-μndel10;
The first term on the right of the 1st equation (krec) is the last term in (10); the second term (kdel) is transfer to the next compartment; the last term (μ) is the natural death rate. The 2nd term on the right of the last equation (kdel) is transfer into the susceptible pool (2).


### 2.1.9. Virus Pool

It is assumed that virus production,svir, which is a surrogate for “infectious strength”, is calculated with(13)svir(t)=kvir(winf1ninf1+winf2ninf2+winf3ninf3+winf4ninf⁡4+wsic1nsic1+wsic2nsic2+wsic3nsic3+wsic4nsic4+wasy2nasy2+wasy3nasy3+…+wasy8nasy8+whos2nhos2+whos3nhos3+whos4nhos4).kvir=200virusunits(weightedinfectedperson)-1day-1.kvir is a production rate constant which is applied to weighted infected persons. The value above gives a basic reproductive ratio R0 (20) of 4.78 (see Section 2.3). The weighting factors w of (13) are (0≤w≤1)(14)winf⁡1=0,winf2=0.5,winf3=1,winf4=1;wsic1=1,wsic2=1,wsic3=1,wsic4=0.5;wasy2=0.2,wasy3=0.4,wasy4=0.4,wasy5=0.4,wasy6=0.4,wasy7=0.2,wasy8=0;whos2=0.2,whos3=0.2,whos4=0.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.” svir of (13) is the input to the virus pool (15).The outputs from the virus pool (15) are death by natural mortality (kvir,dnm, Table 1) and death induced by a suboptimal environment, represented by rate constant kvir,env (18). The value of kvir,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 kvir,dnm. The differential equation for the virus pool nvir is(15)dnvirdt=svir-(kvir,dnm+kvir,env)nvir.


### 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 nded pool (Figure 1) is(16)Ided=Isic→ded+Ihos→ded.
There are no outputs. Therefore,(17)dndeddt=Ided,t=0,nded=0.
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 Mortalitykvir,env and Seasonal Forcing

The environmentally dependent functionkvir,env in (15) (Figure 1) is assumed here to depend only on daily mean values of air temperature, Tair. 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, kvir,eny (d-1; Figure 1) is written as(18)kvir,env=senvkTair[(Tair-Tair0)+ABS(Tair-Tair0)2]qT,senv=0(or1),kTair=0.1(C°)-2d-1,Tair0=13C°,qT=2.senv is used to switch environmental effects off (0) or on (1). Air temperature Tair above a threshold Tair0 increases virus mortality according to power qT. Parameter qT is assumed equal to two giving a quadratic dependence of temperature above the threshold Tair0 (Table 1). It is usual for the biological effects of temperature to be nonlinear, sometimes approximating to exponential, as in the use of a Q10 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). kTair 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 temperatureTair varies from c. 3 (January) to 17°C (July) ([44], p. 270; [49] p. 142). It is assumed that Tair varies sinusoidally, with(19)Tair=Tair,mean+Tair,var{sin⁡[2π365(iJulian-τTair)]};Tair,mean=10°C(varied),Tair,var=7°C,τTair=115day.
Annual mean and seasonal variation are Tair,mean and Tair,vr. iJulian is the Julian day number (1 on 1 January, leap years are ignored). τTair (d) is the phase of the sinusoid, which is maximum on 25 July. Combining equations (19) with (18) modifies environmentally induced virus death rate, kvir,env, and hence reproductive ratio (20). The same seasonal pattern is applied every year. These equations for temperature are a good approximation to long-term weather means in the UK (Note: the study reported in [34] suggests that, on an annual timescale, cold weather does not predict winter deaths, but, on a shorter term timescale, cold weather could be a significant trigger).


## 2.1.1. Some Definitions

First, define some totals in terms of the state variables (Figure1):(1)ninf⁡=ninf1+ninf2+ninf3+ninf4,nsic=nsic1+nsic2+nsic3+nsic4,nhos=nhos2+nhos3+nhos4,nasy=nasy2+nasy3+⋯+nasy7+nasy8,ndel=ndel1+ndel2+⋯+ndel9+ndel10,nliv=nsus+ninf⁡+nsic+nhos+nasy+nrec+ndel,nfer=nliv=ncity,niash=ninf⁡+nasy+nsic+nhos.
The notation of the first five equations is obvious. Total live population nliv is given by the sixth equation and is assumed equal to the fertile population nfer (seventh equation); both of these are equal to the city population, ncity ((2) and following paragraph). Finally, total infected population niash 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 ncity.State variables are next treated by categories and pools.


## 2.1.2. Susceptibles

The differential equation fornsus is(2)dnsus(t)dt=μnfer+kdelndel10+Ided-μnsus-βnsusnvirnliv.
The two principal inputs are first from births with the number of fertile persons nfer 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 nsus pool is included, Ided. This is the death rate caused by influenza (Figure 1). Ided is given by (16). This is done so that the live person number, nliv (1), remains constant at value ncity (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 −μnsus term in (2)). Also, the instantaneous death rate from influenza, Ided (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 nsus, the virus quantity nvir and is divided by the number of live persons nliv (1). All terms in (2) scale equally with population size.


## 2.1.3. Infecteds

The differential equations for the four sequential infected pools are(3)dninf1dt=βnsusnvirnliv-kinf⁡ninf1-μninf1,dninf2dt=(1-fasy)kinf⁡ninf1-kinf⁡ninf2-μninf2,dninf3dt=kinf⁡ninf2-kinf⁡ninf3-μninf3,dninf4dt=kinf⁡ninf3-kinf⁡ninf4-μninf4.
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 kinfninf1 from the inf1 pool is partly asymptomatic with fraction fasy (Figure 1, Table 1; (7)), the remainder entering the inf2 pool (second equation, first term). The third and fourth equations are straightforward. Rate constant kinf is 2 day−1, such that with four pools the mean time from infection to clinical sickness is 4/kinf=2 days [45]. This gives a gamma-distributed lag for the overall exit time from the fourth pool (Appendix A, (A.4), Figure 9.)


## 2.1.4. Nonhospitalized Clinically Sick

Persons who become clinically sick enter the sic1 pool (Figure1). A fraction fhos (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(4)dnsic1dt=kinf⁡ninf4-ksicnsic1-(μ+μsic)nsic1,dnsic2dt=(1-fhos)ksicnsic1-ksicnsic2-(μ+μsic)nsic2,dnsic3dt=ksicnsic2-ksicnsic3-(μ+μsic)nsic3,dnsic4dt=ksicnsic3-ksicnsic4-(μ+μsic)nsic4.
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 ninf pools above) gamma-distributed (Appendix A, Figure 9, (A.4)).The rate at which individuals are admitted to hospital,Ihos2 (Figure 1), is(5)Ihos2=fhosksicnsic1,
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-influenzanded box of Figure 1; units: persons day-1)(6)Isic→ded=μsic(nsic1+nsic2+nsic3+nsic4).


## 2.1.5. Asymptomatic Infecteds

The differential equations for these seven pools (Figure1) are(7)dnasy2dt=fasykinf⁡ninf1-kasynasy2-μnasy2,dnasy3dt=kasynasy2-kasynasy3-μnasy3,…dnasy8dt=kasynasy7-kasynasy8-μnasy8.
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(8)dnhos2dt=fhosksicnsic1-khosnhos2-(μ+μhos)nhos2,dnhos3dt=khosnhos2-khosnhos3-(μ+μhos)nhos3,dnhos4dt=khosnhos3-khosnhos4-(μ+μhos)nhos4.
The input term fhosksicnsic1 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 fhos 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(9)Ihos→ded=μhos(nhos2+nhos3+nhos4).


## 2.1.7. Recovered and Immune

The state variablenrec is governed by the equation(10)dnrecdt=kasynasy8+ksicnsic4+khosnhos4-μnrec-krecnrec.
The positive (input) terms are from the last equations of (7), (4), and (8). Death occurs at its natural rate (μ). The rate constant krec 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 krec causes rapid entry into the delay sequence of pools (Figure 1) which leads to loss of immunity.


## 2.1.8. Delay Pipe Representing Loss of Immunity

Loss of immunity might arise from loss of immunological memory, or from drift and shift in the antigenic character of the virus [11]. This process is represented by ten sequential compartments giving a gamma function delay (Appendix A; also e.g., [44], pp. 818–822). If τimm (day) is the time period during which immunity is lost (three years is assumed; Table 1), and kdel (day−1) is the rate constant out of each of the ten compartments, then(11)kdel=10τimm.
The standard deviation in the exit times is (√10)/kdel and the coefficient of variation is 1/√10 (=standard deviation/τimm) (A.7). The differential equations for the pools are(12)dndel1dt=krecnrec-kdelndel1-μndel1,dndel2dt=kdelndel1-kdelndel2-μndel2,…dndel9dt=kdelndel8-kdelndel9-μndel9,dndel10dt=kdelndel9-kdelndel10-μndel10;
The first term on the right of the 1st equation (krec) is the last term in (10); the second term (kdel) is transfer to the next compartment; the last term (μ) is the natural death rate. The 2nd term on the right of the last equation (kdel) is transfer into the susceptible pool (2).


## 2.1.9. Virus Pool

It is assumed that virus production,svir, which is a surrogate for “infectious strength”, is calculated with(13)svir(t)=kvir(winf1ninf1+winf2ninf2+winf3ninf3+winf4ninf⁡4+wsic1nsic1+wsic2nsic2+wsic3nsic3+wsic4nsic4+wasy2nasy2+wasy3nasy3+…+wasy8nasy8+whos2nhos2+whos3nhos3+whos4nhos4).kvir=200virusunits(weightedinfectedperson)-1day-1.kvir is a production rate constant which is applied to weighted infected persons. The value above gives a basic reproductive ratio R0 (20) of 4.78 (see Section 2.3). The weighting factors w of (13) are (0≤w≤1)(14)winf⁡1=0,winf2=0.5,winf3=1,winf4=1;wsic1=1,wsic2=1,wsic3=1,wsic4=0.5;wasy2=0.2,wasy3=0.4,wasy4=0.4,wasy5=0.4,wasy6=0.4,wasy7=0.2,wasy8=0;whos2=0.2,whos3=0.2,whos4=0.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.” svir of (13) is the input to the virus pool (15).The outputs from the virus pool (15) are death by natural mortality (kvir,dnm, Table 1) and death induced by a suboptimal environment, represented by rate constant kvir,env (18). The value of kvir,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 kvir,dnm. The differential equation for the virus pool nvir is(15)dnvirdt=svir-(kvir,dnm+kvir,env)nvir.


## 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 nded pool (Figure 1) is(16)Ided=Isic→ded+Ihos→ded.
There are no outputs. Therefore,(17)dndeddt=Ided,t=0,nded=0.
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 Mortalitykvir,env and Seasonal Forcing

The environmentally dependent functionkvir,env in (15) (Figure 1) is assumed here to depend only on daily mean values of air temperature, Tair. 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, kvir,eny (d-1; Figure 1) is written as(18)kvir,env=senvkTair[(Tair-Tair0)+ABS(Tair-Tair0)2]qT,senv=0(or1),kTair=0.1(C°)-2d-1,Tair0=13C°,qT=2.senv is used to switch environmental effects off (0) or on (1). Air temperature Tair above a threshold Tair0 increases virus mortality according to power qT. Parameter qT is assumed equal to two giving a quadratic dependence of temperature above the threshold Tair0 (Table 1). It is usual for the biological effects of temperature to be nonlinear, sometimes approximating to exponential, as in the use of a Q10 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). kTair 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 temperatureTair varies from c. 3 (January) to 17°C (July) ([44], p. 270; [49] p. 142). It is assumed that Tair varies sinusoidally, with(19)Tair=Tair,mean+Tair,var{sin⁡[2π365(iJulian-τTair)]};Tair,mean=10°C(varied),Tair,var=7°C,τTair=115day.
Annual mean and seasonal variation are Tair,mean and Tair,vr. iJulian is the Julian day number (1 on 1 January, leap years are ignored). τTair (d) is the phase of the sinusoid, which is maximum on 25 July. Combining equations (19) with (18) modifies environmentally induced virus death rate, kvir,env, and hence reproductive ratio (20). The same seasonal pattern is applied every year. These equations for temperature are a good approximation to long-term weather means in the UK (Note: the study reported in [34] suggests that, on an annual timescale, cold weather does not predict winter deaths, but, on a shorter term timescale, cold weather could be a significant trigger).


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

Two important epidemiological parameters are basic reproductive ratio,R0, and the disease generation time, τgen (day). These are both derived from the basic parameters of the model. R0 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 R0 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 caseR0 can only be calculated numerically), R0 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 R0) (20)R0=βkvir,nm+kvir,envkvirkinf⁡kinf⁡+μ×{(kasykasy+μ)6(khoskhos+μ+μhos)3winf1kinf⁡+(1-fasy)×〈(khoskhos+μ+μhos)3kinf⁡kinf⁡+μwinf2kinf⁡+(kinf⁡kinf⁡+μ)2winf3kinf⁡+(kinf⁡kinf⁡+μ)3×⟦(khoskhos+μ+μhos)3winf4kinf⁡+ksicksic+μ+μsic×{(khoskhos+μ+μhos)3(kasykasy+μ)6wsic1ksic+(1-fhos)[ksicksic+μ+μsicwsic2ksic+(ksicksic+μ+μsic)2wsic3ksic+(ksicksic+μ+μsic)3wsic4ksic]+fhos[khoskhos+μ+μhoswhos2khos+(khoskhos+μ+μhos)2whos3khos+(khoskhos+μ+μhos)3whos4khos]}⟧〉+fasykasykasy+μ[wasy2kasy+(kasykasy+μ)wasy3kasy+(kasykasy+μ)2wasy4+⋯+(kasykasy+μ)6wasy8kasy](khoskhos+μ+μhos)3}.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 nvir to ninf1 and back to nvir (Figure 1). The first term is(21)βkvir,nm+kvir,envwinf⁡1kvirkinf⁡kinf⁡kinf⁡+μ.
Start with a single virus particle in the nvir box (we could equally well start with a single person in the inf1 first-infected compartment). The first factor of (21), with units of persons per virus particle, is the probable number of persons infected by a virus particle during its life; the mean lifetime of a virus particle is 1/(kvir,dnm+kvir,env). It is assumed that there is no additional virus death process due to exposure. winf1 is a dimensionless weighting factor (14). The third factor (kvir/kinf) is the number of virus particles produced per person in the inf1 state during his life: 1/kinf is the average lifetime of a person in the inf1 state. The last factor, with kinf+μ in the denominator, is the probability that this lifetime (of 1/kinf) 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 R0 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 ninf1 box and diverting the (primary) infected individuals into an accumulator, rather than allowing them to enter the ninf1 box, where they can lead to secondary infections. The algebraic and numerical methods agree to six decimal digits, giving a basic reproductive ratio R0=4.78 for the default parameters without forcing.A “dynamic” reproductive ratio,Rdyn, can be calculated when the total population is not entirely susceptible by means of(22)Rdyn=R0nsusnliv.
The live population, nliv, is given by (1). This allows approximately for a slowly changing fraction of susceptible individuals. However, R0 may itself be dynamic on a shorter time scale due to seasonal effects on virus death rate kvir,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 niash (1) increase initially at an exponential rate according to(23)niash~et/τ,
where τ (day) is the growth time constant (the proportional growth rate of infecteds is 1/τ day−1). This depends on both R0 and τgen and the structure of the model [6, Section  1.2.3]. This can be extracted numerically (there is a period of constant exponential growth that lasts for some 10 days) to give τ = 1.29 day. This is considerably less than the generation time of 2.51 day.


## 2.3. Parameterization

Parameters have been introduced while developing the model. Their values are listed for reference in Table1. Here we summarize the evidential basis for the parameter values used.As a preliminary, statements from two physicians are quoted. In [50] a general practitioner in the Doncaster (UK) area describes his study of the 1969-1970 pandemic as it affected his urban practice. He said “the true incidence of influenza during an epidemic is probably impossible to assess.” This is perhaps equally true today and sets the scene for the significance of processes such as “validation” and data fitting (see Figure 5 and Section 3.4). Another general practitioner, this time in Kent (UK) [51], states that “an epidemic of influenza tends to last in this area for between two and three months. Beginning slowly the epidemic reaches its peak in four to five weeks and then subsides slowly. The extent and severity of any attack will depend on factors such as the strain of influenza virus, on the state of the host-immunity of the population and on the timing of the epidemic; the fatality and complication rates are always higher during the cold and foggy winter months.” His statement agrees with many of our seasonal-forcing simulations (Figures 6 and 7).Five key epidemiological quantities for influenza are (a) the time which elapses after the infection event until the subject becomes infectious to others, denoted by theτinf0 (day); (b) the latent period, namely the time which elapses between the infection event and the appearance of clinical sickness, denoted by τlat (day); (c) the infectious period τinf (day), which is the time during which the subject is infectious (Figure 1—producing virus); (d) the period of immunity τimm (day)—the time period after recovery from influenza during which the subject has immunity, before gradually losing it and returning to the susceptible pool (Figure 1, τimm is represented by the box at the bottom of the diagram with 10 sequential pools); (e) the basic reproductive ratio R0 (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 ninf1 with a mean lifetime of 0.5 day is not infectious but the second infected pool ninf2 is, and therefore(24)τinf⁡0=1(kinf⁡=2day-1)=0.5day.
Next the latent period τlat (day) is given by (there are four sequential pools in the infected category of Figure 1)(25)τlat=4(kinf⁡=2day-1)=2day.
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 kinf,kasy,ksic, and khos are all equal (Table 1), say to kdis; (c) the first infected pool (ninf1) is assumed not infectious and all sick and hospitalized pools are infectious (14). We therefore write the infectious period τinf as(26)τinf⁡=7(kdis=2day-1)=3.5day.
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 kdel (day−1) (ignoring natural mortality (12)). This gives a gamma-distributed delay (Appendix A, (A.5) to (A.8) with m=10). τimm and kdel are related by(27)kdel=10τimmwithτimm=1095day=3year.
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 winf2=wasy2=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 ratioR0 (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 R0 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 R0 = 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 R0 and therefore will underestimate the level of global control measures … that will be needed ….” [52] reviews R0 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 R0 = 4.78 (20). This can be regarded as an upper baseline applying to an optimum environment (one maximizing R0), 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 Figure1, 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,Ihos2 (persons day−1); Figure 1), as this is often a recorded statistic. An “epidemic” is deemed to have occurred if Ihos2 shows a maximum (in steps of Δt) and if at the maximum Ihos2≥0.1 persons day-1.State variable initial values and parameter values are as listed in Table1 (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 loss of individual immunity or seasonal forcing. There is no loss of individual immunity (krec=0; (10), Figure 1) nor seasonal forcing ((18); senv=0). (a) Short-term dynamics. (b) Long-term dynamics. Parameters and initial conditions have the values in Table 1. Successive epidemics are caused by loss of immunity at the population level from the natural birth and death processes (μ, (2), Table 1).
(a)(b)Effect of varying key parameters. Parameters and initial conditions have the values in Table1. There is no loss of individual immunity (krec=0; (10), Figure 1) and no seasonal forcing ((18); senv=0). In each case, the value of the altered parameter, basic reproductive ratio R0 (20), and total hospital admissions (HA, the integral of Ihos2, (5)) are given; the curve for the default parameters is thickened and indicated. (a) Changes to infection parameter β. (b) Changes to kinf(=ksic=khos=kasy; Figure 1). (c) Changes to the initial (time t = 0) immune fraction in the recovered category, frec0. This is achieved by altering the time t = 0 value of nrec state variable (Figure 1) to nrec(t=0)=frec0ncity (Table 1, (22)).
(a)(b)(c)Dynamics with loss of individual immunity. Individuals now lose immunity (krec=1; Figure 1). Parameters and initial conditions have the values in Table 1. (a) Short-term dynamics: susceptible nsus and delayed ndel numbers are calculated with (2) and (1), the dynamic reproductive ratio Rdyn with (22), and hospital admission rate Ihos2 with (5). (b) Longer-term dynamics.
(a)(b)Figure 5
Model validation. The model is validated by comparing predictions (solid line) with data (●) from [53]. See text for details.Environmental effects on basic reproductive ratio,R0. (a) Seasonal change in virus inactivation rate, kvir,env (Figure 1; (15), (18)), giving virus inactivation by temperature alone. Tair is mean daily air temperature (appropriate to southern Britain) (19); Tair0 is an assumed temperature threshold for virus inactivation, here equal to 13°C ((18), Table 1)); R0 is the seasonally modulated basic reproductive ratio (20), calculated assuming temperature is constant. (b) For different values of mean annual air temperature, Tair,mean (19), basic reproductive ratio R0 (20) is seasonally modulated as in (a) via the seasonal modulation of kvir,env, for a constant threshold temperature for virus inactivation of Tair0 = 13°C (18), giving for R0 a maximum, minimum, annual mean, and amplitude (maximum-minimum). Nfdy denotes the number of forcing days per year (a “forcing day” is a day on which air temperature is sufficiently high as to increase virus mortality).
(a)(b)Weak environmental forcing. Seasonal forcing (18) combined with loss of individual immunity (Figure 4; Figure 1, krec=1 day−1). The system is initially in a steady state without forcing (Figure 4(b)); forcing is applied after 1 year. The effects of seasonal forcing on the basic reproductive rate R0 (20) are shown in parts (a)–(c) (right-hand y-axis). Increases in mean annual air temperature Tair,mean increase the level of forcing ((18), (19), Figure 6(a)). R̅0 is the mean annual value of R0. Parts (c) and (d) are from a single simulation and have the same level of forcing.
(a)(b)(c)(d)Moderate environmental forcing. Seasonal forcing (18) combined with loss of individual immunity (Figure 4; Figure 1, krec=1 day−1). The system is initially in a steady state without forcing (Figure 4(b)); forcing is applied after 1 year. The effects of seasonal forcing on the basic reproductive rate R0 (20) are shown (right-hand y-axis). Increases in mean annual air temperature Tair,mean increase the level of forcing ((18), (19), Figure 6(a)). R̅0 is the mean annual value of R0.
(a)(b)(c)(d)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 (b2), three pools (c3), four pools (d4), and eight pools (h8) 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 km has been set to give the same mean overall transit time of m/km=4 day (indicated by ■). For each curve, the maximum value (●) denotes the most probable overall transit time.


### 3.2. Dynamics without Intrinsic Loss of Immunity and without Forcing

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


#### 3.2.1. Short-Term Dynamics

These are illustrated in Figure2(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 (ninf) peaks a few days before hospitalized numbers (nhos), both then falling to very low values. Most (99%) of the population joins the recovered (immune) box (nrec), with 1% escaping infection altogether. 90% of those infected travel via the asymptomatic route (Figure 1; fasy, 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 Figure2(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 Rdyn (22) (not plotted) falling to near zero, as the fraction of susceptibles (nsus/nliv) becomes small. In the first epidemic, at time t≈20 day, the initial peak in hospital admissions Ihos2=949 (Figure 1, (5)), decreasing to less than 100 at the second epidemic. After the first epidemic, Rdyn (22) slowly recovers as immune numbers (nrec) decrease and susceptibles (nsus) 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 Rdyn (22) = 1.


#### 3.2.3. Responses to Key Parameters

Figure3 illustrates the effects of changing three key parameters: infectivity parameter β (Figure 1; (2); Table 1), rate parameter kinf (with kinf=ksic=khos=kasy; Figure 1; Table 1; (3)), and the initial (time zero) immune fraction, frec0. 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 kvir; both are linear factors of R0 (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 Ihos2 (e.g., Figure 2(b)) becomes closer together.Figure3(b), where the rates of transit of infected persons through the system are increased (kinf=ksic=khos=kasy=kxxx, Figure 1), is less straightforward. Here a low value for the kxxx (e.g., 0.25 day−1) gives a high basic reproductive ratio R0 (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 Ihos2 spikes as in Figure 2(b) are closer together). As the kxxx are increased basic reproductive ratio R0 decreases (20) first the epidemic becomes narrower and faster (e.g., kxxx=2 day−1), but further increases in kxxx (e.g., kxxx=5 and 8 day−1) cause the epidemic to become slower to take hold and less peaked, with fewer hospital admissions (HA). Increasing kxxx always causes the Ihos2 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 (niash,(1)).Finally, Figure3(c) illustrates how increasing the initial immune fraction freco 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 Figure1 giving loss of immunity is switched on by makingkrec=1 day−1, which causes recovered individuals to be moved quite rapidly into the delay sequence of ten compartments where immunity is lost after three years (τimm=3×365days; Table 1, (11)). Immunity is also being lost due to the natural birth and death processes, as newborns are assumed to be susceptible. The responses of Figure 2(a) are little changed by taking krec=1 day−1 (instead of 0) if the nrec variable of Figure 2(a) is replaced by the ndel 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 ncity, Table 1, (1) and following paragraph), with an inverse switching of numbers in the delay compartment, ndel (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 Rdyn = 1 (22), and the fractions in the three principal categories (Figure 1) of susceptible; all infected categories lumped together (infected, asymptomatic, sick, hospitalized, niash (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 Figure1 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. Table1 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 ncity=763; and the asymptomatic fraction, fasy, 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 R0 (20) is 6.44, the mean generation time is 2.65 day (A.7), and the initial proportional growth rate of total infecteds, niash(1), is 1.015 day−1.For comparison with the data in [53], we define the number of persons confined to bed, nctb, as(28)nctb(t)=nsic1(t)+nsic2(t)+nsic3(t)+nsic4(t)+nhos2(t)+nhos3(t)+nhos4(t).
Since this definition is somewhat arbitrary, in comparing the predicted nctb from the model with the data from Anon (1978), nctb(t) is scaled with an adjustable factor so that the comparison line drawn in Figure 5 is(29)1.35nctb(t).
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 Figure5 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 Figure1 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 R0 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 Figure4 ((18), (19)). Four weather factors which could impinge on virus longevity are air temperature (Tair), relative humidity (RH), radiation (possibly multicomponent), and wind speed. The effects can be highly complex: for example, [48] which examines the effects of Tair and RH 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 Tair alone is operative.Figure6(a) illustrates the seasonal variation of mean daily air temperature Tair in the southern UK. Using (19) with (18), and senv=1, this affects environmentally induced virus death rate kvir,env (Figure 1, (15)), and thereby basic reproductive ratio R0 (20). Mean daily air temperature Tair varies between 3°C (Jan 24) and 17°C (Jul 25). When Tair crosses the temperature threshold of Tair0 = 13°C in May, this increases environmentally induced virus mortality, kvir,env (18) and decreases basic reproductive ratio R0 (20). With this formulation, influenza is most likely during the months of October to April when R0 is highest.Figure6(b) shows how varying (e.g., increasing) mean annual air temperature, Tair,mean (19), (or equivalently, decreasing threshold temperature Tair0, (18)) varies the duration and intensity of seasonal forcing of R0. Changing mean annual air temperature, Tair,mean (19), changes the average annual value of R0, R̅0, as well as its maximum, minimum and amplitude (maximum-minimum). The dependence of these quantities on Tait,mean is also shown, together with the number of forcing days per year (Nfdy, a forcing day is one on which mean air temperature is above threshold temperature and virus mortality is reduced). With mean annual air temperature Tair,mean≤6°C, there is no forcing at all (Nfdy=0): air temperature Tair (19) never exceeds threshold temperature Tair0 = 13°C (18); R0 is invariant at its maximum value. Environmentally induced virus death rate kvir,env ((15), (18)) stays at zero; the system remains in the steady state shown in Figure 4(b) and R̅0=R0 (this situation is equivalent to Tair,mean=10°C and Tair0 = 17°C). As mean annual air temperature Tair,mean increases above 6°C, there are more days in the summer months when mean daily temperature Tair> threshold temperature Tair0, number of forcing days per year Nfdy increases, and R̅0<R0. The amplitude increases to a maximum whenTair,mean=21°C, before decreasing. With Tair,mean>28.5°C, influenza epidemics do not occur, because the annual average of R0 is less than unity.In the simulations presented in Figures7 and 8, various values of mean annual air temperature Tair,mean are taken between −4.5 and 16°C, assuming always an annual variation Tair,var of ±7°C (19) as in the UK. Otherwise, parameters have the values in Table 1. In these simulations, the long-term steady state reached in Figure 4, with loss of immunity (krec=1; Figure 1; (10)) but without seasonal forcing, is used for initial values. Forcing is applied after one calendar year. The aim is to illustrate the wide variety in dynamic behaviour which results from this type of seasonal forcing (which decreases the reproductive ratio). In each case, the mean value, its maximum, minimum and amplitude (maximum-minimum) of basic reproductive ratio can be read off Figure 6(b). In the UK the influenza season is considered to be over by May, when the mean daily temperature ranges from 10.7 (1 May) to 14.1°C (31 May) (Figure 6). Therefore, of the simulations described in Figures 7 and 8, those given in Figures 7(b), 7(c),7(d), 8(a), 8(b), and 8(c) are more relevant to the UK.Figure7 illustrates the effects of low levels of seasonal forcing on influenza hospital admissions, Ihos2 (Figure 1,(5)). Forcing is seen as a downward modulation of the basic reproductive ratio R0, which decreases the annual average, R̅0.The lower graph in Figure7(a) shows that R0 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 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 Figure7(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 (Figures7(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 Figure8. First, in Figure 8(a) with mean annual air temperature Tair,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); Tair,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); Tair,mean=13°C), there is chaos but now with a tendency towards biennial epidemics (169 epidemics occur during 200 years). Last (Figure 8(d); Tair,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 Figures7 and 8, as mean annual air temperature threshold Tair,mean is increased, forcing is increased (i.e., the magnitude of the seasonal changes in basic reproductive ratio R0 ((20), (18), (19), Figure 6) increases), but the mean annual value of R0, 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 temperatureTair,mean can be continued, and although the situations simulated are now less realistic, they do contribute to understanding the system. With Tair,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 Tair,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 Tair,mean=26°C the response is chaotic with 15 epidemics in the first 200 years and R̅0=1.367. Last, with Tair,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.


## 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,Ihos2 (persons day−1); Figure 1), as this is often a recorded statistic. An “epidemic” is deemed to have occurred if Ihos2 shows a maximum (in steps of Δt) and if at the maximum Ihos2≥0.1 persons day-1.State variable initial values and parameter values are as listed in Table1 (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 loss of individual immunity or seasonal forcing. There is no loss of individual immunity (krec=0; (10), Figure 1) nor seasonal forcing ((18); senv=0). (a) Short-term dynamics. (b) Long-term dynamics. Parameters and initial conditions have the values in Table 1. Successive epidemics are caused by loss of immunity at the population level from the natural birth and death processes (μ, (2), Table 1).
(a)(b)Effect of varying key parameters. Parameters and initial conditions have the values in Table1. There is no loss of individual immunity (krec=0; (10), Figure 1) and no seasonal forcing ((18); senv=0). In each case, the value of the altered parameter, basic reproductive ratio R0 (20), and total hospital admissions (HA, the integral of Ihos2, (5)) are given; the curve for the default parameters is thickened and indicated. (a) Changes to infection parameter β. (b) Changes to kinf(=ksic=khos=kasy; Figure 1). (c) Changes to the initial (time t = 0) immune fraction in the recovered category, frec0. This is achieved by altering the time t = 0 value of nrec state variable (Figure 1) to nrec(t=0)=frec0ncity (Table 1, (22)).
(a)(b)(c)Dynamics with loss of individual immunity. Individuals now lose immunity (krec=1; Figure 1). Parameters and initial conditions have the values in Table 1. (a) Short-term dynamics: susceptible nsus and delayed ndel numbers are calculated with (2) and (1), the dynamic reproductive ratio Rdyn with (22), and hospital admission rate Ihos2 with (5). (b) Longer-term dynamics.
(a)(b)Figure 5
Model validation. The model is validated by comparing predictions (solid line) with data (●) from [53]. See text for details.Environmental effects on basic reproductive ratio,R0. (a) Seasonal change in virus inactivation rate, kvir,env (Figure 1; (15), (18)), giving virus inactivation by temperature alone. Tair is mean daily air temperature (appropriate to southern Britain) (19); Tair0 is an assumed temperature threshold for virus inactivation, here equal to 13°C ((18), Table 1)); R0 is the seasonally modulated basic reproductive ratio (20), calculated assuming temperature is constant. (b) For different values of mean annual air temperature, Tair,mean (19), basic reproductive ratio R0 (20) is seasonally modulated as in (a) via the seasonal modulation of kvir,env, for a constant threshold temperature for virus inactivation of Tair0 = 13°C (18), giving for R0 a maximum, minimum, annual mean, and amplitude (maximum-minimum). Nfdy denotes the number of forcing days per year (a “forcing day” is a day on which air temperature is sufficiently high as to increase virus mortality).
(a)(b)Weak environmental forcing. Seasonal forcing (18) combined with loss of individual immunity (Figure 4; Figure 1, krec=1 day−1). The system is initially in a steady state without forcing (Figure 4(b)); forcing is applied after 1 year. The effects of seasonal forcing on the basic reproductive rate R0 (20) are shown in parts (a)–(c) (right-hand y-axis). Increases in mean annual air temperature Tair,mean increase the level of forcing ((18), (19), Figure 6(a)). R̅0 is the mean annual value of R0. Parts (c) and (d) are from a single simulation and have the same level of forcing.
(a)(b)(c)(d)Moderate environmental forcing. Seasonal forcing (18) combined with loss of individual immunity (Figure 4; Figure 1, krec=1 day−1). The system is initially in a steady state without forcing (Figure 4(b)); forcing is applied after 1 year. The effects of seasonal forcing on the basic reproductive rate R0 (20) are shown (right-hand y-axis). Increases in mean annual air temperature Tair,mean increase the level of forcing ((18), (19), Figure 6(a)). R̅0 is the mean annual value of R0.
(a)(b)(c)(d)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 (b2), three pools (c3), four pools (d4), and eight pools (h8) 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 km has been set to give the same mean overall transit time of m/km=4 day (indicated by ■). For each curve, the maximum value (●) denotes the most probable overall transit time.


## 3.2. Dynamics without Intrinsic Loss of Immunity and without Forcing

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


### 3.2.1. Short-Term Dynamics

These are illustrated in Figure2(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 (ninf) peaks a few days before hospitalized numbers (nhos), both then falling to very low values. Most (99%) of the population joins the recovered (immune) box (nrec), with 1% escaping infection altogether. 90% of those infected travel via the asymptomatic route (Figure 1; fasy, 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 Figure2(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 Rdyn (22) (not plotted) falling to near zero, as the fraction of susceptibles (nsus/nliv) becomes small. In the first epidemic, at time t≈20 day, the initial peak in hospital admissions Ihos2=949 (Figure 1, (5)), decreasing to less than 100 at the second epidemic. After the first epidemic, Rdyn (22) slowly recovers as immune numbers (nrec) decrease and susceptibles (nsus) 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 Rdyn (22) = 1.


### 3.2.3. Responses to Key Parameters

Figure3 illustrates the effects of changing three key parameters: infectivity parameter β (Figure 1; (2); Table 1), rate parameter kinf (with kinf=ksic=khos=kasy; Figure 1; Table 1; (3)), and the initial (time zero) immune fraction, frec0. 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 kvir; both are linear factors of R0 (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 Ihos2 (e.g., Figure 2(b)) becomes closer together.Figure3(b), where the rates of transit of infected persons through the system are increased (kinf=ksic=khos=kasy=kxxx, Figure 1), is less straightforward. Here a low value for the kxxx (e.g., 0.25 day−1) gives a high basic reproductive ratio R0 (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 Ihos2 spikes as in Figure 2(b) are closer together). As the kxxx are increased basic reproductive ratio R0 decreases (20) first the epidemic becomes narrower and faster (e.g., kxxx=2 day−1), but further increases in kxxx (e.g., kxxx=5 and 8 day−1) cause the epidemic to become slower to take hold and less peaked, with fewer hospital admissions (HA). Increasing kxxx always causes the Ihos2 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 (niash,(1)).Finally, Figure3(c) illustrates how increasing the initial immune fraction freco 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.2.1. Short-Term Dynamics

These are illustrated in Figure2(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 (ninf) peaks a few days before hospitalized numbers (nhos), both then falling to very low values. Most (99%) of the population joins the recovered (immune) box (nrec), with 1% escaping infection altogether. 90% of those infected travel via the asymptomatic route (Figure 1; fasy, 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 Figure2(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 Rdyn (22) (not plotted) falling to near zero, as the fraction of susceptibles (nsus/nliv) becomes small. In the first epidemic, at time t≈20 day, the initial peak in hospital admissions Ihos2=949 (Figure 1, (5)), decreasing to less than 100 at the second epidemic. After the first epidemic, Rdyn (22) slowly recovers as immune numbers (nrec) decrease and susceptibles (nsus) 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 Rdyn (22) = 1.


## 3.2.3. Responses to Key Parameters

Figure3 illustrates the effects of changing three key parameters: infectivity parameter β (Figure 1; (2); Table 1), rate parameter kinf (with kinf=ksic=khos=kasy; Figure 1; Table 1; (3)), and the initial (time zero) immune fraction, frec0. 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 kvir; both are linear factors of R0 (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 Ihos2 (e.g., Figure 2(b)) becomes closer together.Figure3(b), where the rates of transit of infected persons through the system are increased (kinf=ksic=khos=kasy=kxxx, Figure 1), is less straightforward. Here a low value for the kxxx (e.g., 0.25 day−1) gives a high basic reproductive ratio R0 (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 Ihos2 spikes as in Figure 2(b) are closer together). As the kxxx are increased basic reproductive ratio R0 decreases (20) first the epidemic becomes narrower and faster (e.g., kxxx=2 day−1), but further increases in kxxx (e.g., kxxx=5 and 8 day−1) cause the epidemic to become slower to take hold and less peaked, with fewer hospital admissions (HA). Increasing kxxx always causes the Ihos2 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 (niash,(1)).Finally, Figure3(c) illustrates how increasing the initial immune fraction freco 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 Figure1 giving loss of immunity is switched on by makingkrec=1 day−1, which causes recovered individuals to be moved quite rapidly into the delay sequence of ten compartments where immunity is lost after three years (τimm=3×365days; Table 1, (11)). Immunity is also being lost due to the natural birth and death processes, as newborns are assumed to be susceptible. The responses of Figure 2(a) are little changed by taking krec=1 day−1 (instead of 0) if the nrec variable of Figure 2(a) is replaced by the ndel 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 ncity, Table 1, (1) and following paragraph), with an inverse switching of numbers in the delay compartment, ndel (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 Rdyn = 1 (22), and the fractions in the three principal categories (Figure 1) of susceptible; all infected categories lumped together (infected, asymptomatic, sick, hospitalized, niash (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 Figure1 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. Table1 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 ncity=763; and the asymptomatic fraction, fasy, 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 R0 (20) is 6.44, the mean generation time is 2.65 day (A.7), and the initial proportional growth rate of total infecteds, niash(1), is 1.015 day−1.For comparison with the data in [53], we define the number of persons confined to bed, nctb, as(28)nctb(t)=nsic1(t)+nsic2(t)+nsic3(t)+nsic4(t)+nhos2(t)+nhos3(t)+nhos4(t).
Since this definition is somewhat arbitrary, in comparing the predicted nctb from the model with the data from Anon (1978), nctb(t) is scaled with an adjustable factor so that the comparison line drawn in Figure 5 is(29)1.35nctb(t).
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 Figure5 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 Figure1 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 R0 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 Figure4 ((18), (19)). Four weather factors which could impinge on virus longevity are air temperature (Tair), relative humidity (RH), radiation (possibly multicomponent), and wind speed. The effects can be highly complex: for example, [48] which examines the effects of Tair and RH 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 Tair alone is operative.Figure6(a) illustrates the seasonal variation of mean daily air temperature Tair in the southern UK. Using (19) with (18), and senv=1, this affects environmentally induced virus death rate kvir,env (Figure 1, (15)), and thereby basic reproductive ratio R0 (20). Mean daily air temperature Tair varies between 3°C (Jan 24) and 17°C (Jul 25). When Tair crosses the temperature threshold of Tair0 = 13°C in May, this increases environmentally induced virus mortality, kvir,env (18) and decreases basic reproductive ratio R0 (20). With this formulation, influenza is most likely during the months of October to April when R0 is highest.Figure6(b) shows how varying (e.g., increasing) mean annual air temperature, Tair,mean (19), (or equivalently, decreasing threshold temperature Tair0, (18)) varies the duration and intensity of seasonal forcing of R0. Changing mean annual air temperature, Tair,mean (19), changes the average annual value of R0, R̅0, as well as its maximum, minimum and amplitude (maximum-minimum). The dependence of these quantities on Tait,mean is also shown, together with the number of forcing days per year (Nfdy, a forcing day is one on which mean air temperature is above threshold temperature and virus mortality is reduced). With mean annual air temperature Tair,mean≤6°C, there is no forcing at all (Nfdy=0): air temperature Tair (19) never exceeds threshold temperature Tair0 = 13°C (18); R0 is invariant at its maximum value. Environmentally induced virus death rate kvir,env ((15), (18)) stays at zero; the system remains in the steady state shown in Figure 4(b) and R̅0=R0 (this situation is equivalent to Tair,mean=10°C and Tair0 = 17°C). As mean annual air temperature Tair,mean increases above 6°C, there are more days in the summer months when mean daily temperature Tair> threshold temperature Tair0, number of forcing days per year Nfdy increases, and R̅0<R0. The amplitude increases to a maximum whenTair,mean=21°C, before decreasing. With Tair,mean>28.5°C, influenza epidemics do not occur, because the annual average of R0 is less than unity.In the simulations presented in Figures7 and 8, various values of mean annual air temperature Tair,mean are taken between −4.5 and 16°C, assuming always an annual variation Tair,var of ±7°C (19) as in the UK. Otherwise, parameters have the values in Table 1. In these simulations, the long-term steady state reached in Figure 4, with loss of immunity (krec=1; Figure 1; (10)) but without seasonal forcing, is used for initial values. Forcing is applied after one calendar year. The aim is to illustrate the wide variety in dynamic behaviour which results from this type of seasonal forcing (which decreases the reproductive ratio). In each case, the mean value, its maximum, minimum and amplitude (maximum-minimum) of basic reproductive ratio can be read off Figure 6(b). In the UK the influenza season is considered to be over by May, when the mean daily temperature ranges from 10.7 (1 May) to 14.1°C (31 May) (Figure 6). Therefore, of the simulations described in Figures 7 and 8, those given in Figures 7(b), 7(c),7(d), 8(a), 8(b), and 8(c) are more relevant to the UK.Figure7 illustrates the effects of low levels of seasonal forcing on influenza hospital admissions, Ihos2 (Figure 1,(5)). Forcing is seen as a downward modulation of the basic reproductive ratio R0, which decreases the annual average, R̅0.The lower graph in Figure7(a) shows that R0 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 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 Figure7(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 (Figures7(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 Figure8. First, in Figure 8(a) with mean annual air temperature Tair,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); Tair,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); Tair,mean=13°C), there is chaos but now with a tendency towards biennial epidemics (169 epidemics occur during 200 years). Last (Figure 8(d); Tair,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 Figures7 and 8, as mean annual air temperature threshold Tair,mean is increased, forcing is increased (i.e., the magnitude of the seasonal changes in basic reproductive ratio R0 ((20), (18), (19), Figure 6) increases), but the mean annual value of R0, 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 temperatureTair,mean can be continued, and although the situations simulated are now less realistic, they do contribute to understanding the system. With Tair,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 Tair,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 Tair,mean=26°C the response is chaotic with 15 epidemics in the first 200 years and R̅0=1.367. Last, with Tair,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.


## 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 Figures7 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 Section3.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 (Figure1) (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.



---

*Source: 471653-2012-02-29.xml*