Modelling the outbreak of infectious disease following mutation from a non-transmissible strain.

In-host mutation of a cross-species infectious disease to a form that is transmissible between humans has resulted with devastating global pandemics in the past. We use simple mathematical models to describe this process with the aim to better understand the emergence of an epidemic resulting from such a mutation and the extent of measures that are needed to control it. The feared outbreak of a human-human transmissible form of avian influenza leading to a global epidemic is the paradigm for this study. We extend the SIR approach to derive a deterministic and a stochastic formulation to describe the evolution of two classes of susceptible and infected states and a removed state, leading to a system of ordinary differential equations and a stochastic equivalent based on a Markov process. For the deterministic model, the contrasting timescale of the mutation process and disease infectiousness is exploited in two limits using asymptotic analysis in order to determine, in terms of the model parameters, necessary conditions for an epidemic to take place and timescales for the onset of the epidemic, the size and duration of the epidemic and the maximum level of the infected individuals at one time. Furthermore, the basic reproduction number R0 is determined from asymptotic analysis of a distinguished limit. Comparisons between the deterministic and stochastic model demonstrate that stochasticity has little effect on most aspects of an epidemic, but does have significant impact on its onset particularly for smaller populations and lower mutation rates for representatively large populations. The deterministic model is extended to investigate a range of quarantine and vaccination programmes, whereby in the two asymptotic limits analysed, quantitative estimates on the outcomes and effectiveness of these control measures are established.

Introduction

Cross-species transmission of viruses from animals to human and their subsequent evolution have in the past led to epidemics, some with disastrous results on a global scale. HIV/AIDS, severe acute respiratory syndrome (SARS), West Nile virus and the Spanish flu of 1918 are just a few examples Greger, 2007, Tapper, 2006. Viruses continue to present new threats, a recent example being the ebola outbreaks in West Africa (2014–2016, with 28 616 known cases and 11 310 deaths, as of June 2016) and in the Democratic Republic of the Congo (56 cases and 25 deaths, as of early June 2018) (World Health, 2015-2018). Another example is the avian influenza, in particular the H5N1 virus which has been circulating since 1997 in parts of Asia and some regions of Africa with 860 confirmed cases of human infection (454 deaths) in the years between 2003 and 2018 (World Health Organization, 2018). This disease in its current state does not readily infect humans, but can spread rapidly amongst domesticated birds with a mortality rate as high as  (Alexander, 2000). Human cases have resulted from direct or close contact with infected poultry (domesticated chickens, ducks and turkeys). The spread of avian influenza from an infected person to another has rarely been reported and the transmission has not been observed to continue beyond one person. All influenza viruses, however, have the ability to change and adapt. Past examples, such as the Spanish flu, have shown that this could result with a transmissible virus which can spread easily and directly to humans. Such an outbreak could develop into an epidemic on a global scale, potentially taking millions of lives. A better understanding of the dynamics of how events leading to an epidemic can help to control and eradicate the disease and to prevent an infectious disease from developing into an epidemic.

Viruses are highly prone to genetic mutation, influenza viruses being well-studied examples. Minor mutations may involve small changes to surface glycoproteins, resulting with a strain of the same viral subtype that can deceive the adaptive immune response, and reinfect a previously infected host. This kind of mutation is known as antigenic drift and leads to recurrent, often annual, outbreaks of an influenza subtype. More dramatic mutations leading to new viral subtypes are known as antigenic shifts, which can arise in at least two ways. The first way is known as reassortment in which influenza viruses of two different subtypes infect the same host simultaneously, whereby genetic material could be exchanged to form a new, viable virus. This is a feared possibility regarding individuals infected by both the human and avian forms of influenza, leading to a transmissible form of the disease. The second involves direct or indirect introduction of an avian influenza virus into humans which subsequently adapt to the new host Scholtissek, 1994, Ferguson et al., 2003. The new virus subtype, with no antibodies against it present in the human body, will transmit easily from human to human, as a result, a global pandemic may occur.

Mutation of the highly pathogenic avian influenza virus from a non-transmissible form to one that is rapidly transmissible among humans is treated here as an antigenic shift and not a result of a sequence of antigenic drift events, and hence is modelled as a one-step process taking place within its human host. The indirect transmission of the mutated virus into humans via an intermediate host (e.g. pigs) is not considered here, the model developed, however, can be adopted to include an intermediate host without too much difficulty. The relative simplicity of the system we study enables a detailed analysis to be carried out, from which a number of key characteristics of an epidemic following disease mutation can be quantified in terms of the model’s parameters.

Published models of avian flu have focused on the spread of disease between bird populations and farms Boender et al., 2007, Le Menach et al., 2006, Sharke et al., 2008, Wang et al., 2012 and cross-species transmission between bird and human populations Chong et al., 2014, Gumel, 2009, Iwami et al., 2007, Martcheva, 2014. In each of these studies their aim is to understand the dynamics of the disease and the effectiveness of measures to control in both bird and human populations. Each of the bird–human transmission models is extension of the “standard” SIR compartmental formulation (Susceptibles–Infectives–Removed) for an endemic disease (i.e. natural birth and death processes are considered), whereby mutation to the human–human transmissible form emerges via antigenic shift in human hosts. The mutation events are assumed to have occurred either at a fixed rate Chong et al., 2014, Gumel, 2009, Iwami et al., 2007 or as a spontaneous event at a given time point (Martcheva, 2014). The core of the analysis investigates the stability of equilibria and establishing the basic reproduction numbers of the disease in the populations. In contrast to these studies, we analyse in this paper a model that describes the potential development of an epidemic of human–human transmission form of flu by focusing on the human population only, whereby bird–human transmission only occurs to a sub-population of humans (bird handlers) at a fixed rate (assuming an approximately fixed population of diseased birds). In common with the above mentioned avian flu models, we extend the SIR formulation and assume mutation is via antigenic shift in humans, however, we analyse in much more detail the transient behaviour (e.g. timescales for epidemic outbreak, duration of epidemic) as well as the long-time dynamics. Moreover, a stochastic version of the model is investigated to assess the reliability of the predictions from a continuum model in a range of circumstances. We note that modelling epidemics using SIR based models has had a long history since the seminal work of Kermack and McKendrick in 1927 (Kermack and McKendrick, 1927), including applications to influenza, sexually transmitted disease, the “Bombay plague” (Murray, 1993) and West Nile virus (Wonham et al., 2004). Amongst the mathematical extensions include more general infection rate kinetics (Korobeinikov and Maini, 2004), terms with time delays Jin and Ma, 2006, Tchuenche et al., 2007, Cooke, 1979, stochasticity Tornatore et al., 2005, Tuckwell and Williams, 2007 and spatial spread (Källén et al., 1985). We also note there are a number of epidemic models examining mutations and genetic drift (e.g. Boni et al., 2006, Girvan et al., 2002, Lin et al., 2003, Pease, 1987), many of which are more sophisticated than the proposed model, but at a cost at being more difficult to analyse.

In this paper we propose a simple deterministic model and a stochastic equivalent to describe in-host viral mutation of a cross-species disease that enables human–human transmission. Whilst the former model is expected to be mainly applicable to large populations, the slow mutation process suggests that stochasticity is important, particularly early in an epidemic. We analyse these models to determine the effects of parameters on a number of features of an epidemic (onset time, epidemic duration etc.), establishing the strengths and weaknesses of the two approaches. Furthermore, the deterministic model is extended to investigate potential measures in disease control. Though the model provides a general framework for in-host mutation more virulent forms of disease, we use the much feared avian flu scenario as a paradigm.

The paper is arranged as follows. In the next section, the mathematical model is derived to describe the spread of a bird–human transmissible virus, which can mutate into a human–human transmissible form. Only a fraction of the population is assumed to be at risk of the bird–human transmissible form, but the whole population is vulnerable to the mutated version. We carry out large-time (Section 3) and asymptotic (Section 4) analysis to obtain predictions, in terms of reduced model solutions, on the possible development of an epidemic. In Section 5 the effects of intervention by means of quarantine and vaccination are investigated. The stochastic version of the model is presented and analysed in Section 6 and we conclude with an overview of the result in Section 7.

Deterministic model

The model proposed below aims to describe in-host mutation of a non-human–human transmissible to a human–human transmissible virus strain and its spread throughout a population. The population is assumed well mixed and of size ranging from , i.e. anywhere from large city to continental scale. Though the development of mutant viral forms may take years to develop, an epidemic may last just a few weeks; we assume over the timescale of interest the demographic turnover is negligible so that the total population remains constant. The entire population is broadly compartmentalised into susceptible, infective and recovered (or removed) classes and the populations vary continuously in time. In the case of avian flu, only a fraction of a population are exposed to birds who are themselves exposed to the bird–bird transmissible form of avian flu. We thus split the susceptible class into two, one class that are bird-handlers (mainly poultry farmers) at risk of catching flu directly from their birds and the other class who are not exposed to potentially diseased birds. We note, that this assumption has general application, for example, where sub-populations have different exposure to a non-human–human transmissible disease based on their activity or their geographical location. We will not be modelling bird populations explicitly and assume that their numbers are constant. Those infected by the bird–human transmissible form of avian flu are not able to infect other susceptibles, unless the virus has mutated to a human–human transmissible form, thus two infected classes are assumed. To summarise, the total population is divided into five subclasses:

: population of the bird-handlers at risk of catching avian flu from birds.

: population of individuals not at risk of catching avian flu from birds.

: population of infectives with only the non-mutated bird–human transmissible avian flu.

: populations of infectives with the mutated, human–human transmissible avian flu.

: population of those recovered/removed from either form of avian flu.

with and conventionally denote the susceptible, infective and removed/recovered. The variable is the population who has survived or died from the disease; taking the mortality fraction to be , then is the number of deaths from the disease (currently  World Health Organization, 2015-2018). We note that we have assumed that recovery from one form of the flu leads to immunity against both forms.

The pathway summarising the interaction between the population classes is illustrated in Fig. 1, and a summary of the modelling assumptions is as follows,

Fig. 1

Schematic of the infection and mutation pathways used in the derivation of the model.

The total population is constant.

Only the bird-handler class can (spontaneously) catch avian flu from birds.

Non-mutated avian flu cannot be transmitted from person to person.

becomes through infection by or through mutation of the virus.

and can be infected by .

The two forms of flu are sufficiently similar for the removal rate to be the same.

Once the infected and are recovered, they become immune to both forms.

The above assumptions do not distinguish between a person infected with a mutated form of the disease and one infected with both non-mutated and mutated forms since the latter affects the disease progress in the same way as the former does. Using these assumptions, the following system of ODEs is obtained: where the sum of the populations is constant The rate constants in the system of the equations are

: disease infection rate constant from birds to humans,

: infection rate constant of the human–human transmissible avian flu,

: mutation rate constant from to ,

: recovery/removal rate constant for both and .

We note that although the disease in the bird population is not explicitly modelled, the proposed system (2.1)–(2.5) can be viewed as a limiting case of such a model. If the disease dynamics of the bird population occurs on a much faster timescale than the processes of bird–human transmission and disease mutation (e.g. as in Iwami et al., 2007), then after a transient period the infected bird population will be close to steady-state (assumed stable) whilst human form of the disease continues to evolve. Since the infected bird population is effectively fixed in this scenario, the rate of bird–human transmission will be proportional to the susceptible bird-handler population, with being the constant of proportionality.

To close the system we suppose that initially, at , the first cases of the bird–human transmissible form of the disease arise; formally, the infection rate is defined to be for and for . Consequently, we assume that initially there is nobody in the infected or recovered classes, hence where is the initial population of those not at risk of catching the bird–human transmissible form of the disease.

In contrast to most ordinary differential equation models in epidemiology, the disease is initiated by the bird–human transmission term in Eqs. (2.1), (2.3), rather than an assumed initial infected population. As a consequence, whilst , there is no disease-free state to determine the basic reproduction number using, for example, the next generation operator method. In this paper we will use the following definition,

The basic reproduction number, , is defined as a formulation of the model parameters such that if then , whilst if then .

In the analysis of Section 4 and Appendix A the term is in the context of a small parameter , whereby in the limit a formula for is derived.

Non-dimensionalisation

In the analysis to follow it is convenient to non-dimensionalise the model using the rescalings These rescalings imply that all the non-dimensional dependent variables represent population fractions of each of the classes and that removal/recovery of an infected individual occurs in time. The non-dimensional forms of the rate constants are we expect that , reflects the fact that most individuals with the non-mutated form will be recovered or died before the virus strain has a chance to mutate to the transmissible form, and for the case of avian flu , as so few individuals ( 1000) have caught the disease in the last decade. It turns out in the analysis to follow that is not necessarily equal to the basic reproduction number , as is the case in the standard SIR model (discussed in Section 4.2). The hats will henceforth be dropped for brevity and the dimensionless system of equations are with and the initial conditions are where .

Large-time analysis

It is clear in the way the model is constructed that all susceptibles in the bird-handler class will eventually be removed as and that the only steady-state of the system belongs to the family of solutions where . The steady-state is the population of non-bird-handler susceptibles, , as . The value of , as with the standard SIR model, is dependent on the initial conditions and is discussed further below. Linear stability analysis on this steady-state results with a Jacobian that has a zero eigenvalue, with all others being negative provided that thereby (3.1) is neutrally stable. By analogy to the standard SIR model, this result suggests that the basic reproduction number for this system is given by . This is true when the parameters and are constants, whereby, the mutation happens sufficiently fast for the whole population to be vulnerable to the human–human transmissible form of the disease. However, as discussed in Section 4.2, will change as . Nevertheless, the inequality (3.2) suggests that in large time will tend to a population satisfying the inequality . As expected, the fraction of susceptibles that remain following an epidemic will be reduced on increased virulence .

One of the key aims of the modelling is to try to predict, in terms of the model parameters, the total mortality following a human–human transmissible disease epidemic, i.e. , recalling that is the disease’s mortality fraction. For the standard SIR model this is straightforward to do, resulting with being the solution of an implicit equation Andreasen, 2011, Miller, 2012. A similar approach can be used in the present model, though the details are more complicated. In the analysis below, we will establish a readily computable bound on , from which we can deduce bounds on .

Eliminating using Eqs. (2.9), (2.10) leads to, on integration, Eliminating using (2.9), (2.11) gives which can be solved and, following the substitution of (3.3), to get The case of is not singular in the qualitative sense and leads to the solution ; this case does not produce any new behaviour of interest and will not be discussed further by assuming in what follows that . As a side note, the system (2.1)–(2.4) can be reduced to a nonlinear 2nd-order differential equation by using from (2.10) and Eqs. (3.3), (3.5), namely However, to determine the long time behaviour it is simpler to add Eqs. (2.9), (2.10), (2.11), (2.12) and substituting from (2.10) to yield on integration, Substituting (3.1) and using (3.5) lead to the following equation for which is a solution, where is the integral operator which has been defined to have the property for constant . In the standard SIR model, the steady-state value of , say , is the solution of , hence the integral term results from the additional assumptions in the current model.

Finding the solution satisfying is not possible without full knowledge of the solution , however, upper and lower bounds can be determined. Viewing as a variable and as a function of , we can deduce from (3.7) that since from (3.2). So the solution of increases on increasing . Moreover, for any non-negative since the bracketed term in the integral (3.8) has the same sign as (neglecting the case when this is zero). From (2.10), is a monotonic decreasing functions, so , which means that is bounded by using the property of stated above. Using (3.9), (3.10) the following bound on can be deduced where constants and satisfy and, respectively, both being straightforward to calculate numerically. Of more interest is the bound on , namely where and . Fig. 2 shows plots of the bounds, given by (3.12), and numerically determined values of against and . The upper bound appears to provide a very good approximation of for all values of ; this was found to be true for small , but as approaches , tends to the lower bound. The reason for this is discussed in Section 4.2. The functions that determine the upper and lower bounds are independent of , hence these values are constant in the right-hand figure. As increases, tends smoothly from the upper to the lower bound, with the lower bound seeming to be a very good approximation from .

Fig. 2

Plots of (solid) and the upper and lower bounds determined from (3.12) (dashed lines). Plot (a) shows against (using ) and (b) against (using ). Note the line is super-imposed over the upper bound line in the left-hand plot.

It is worth noting from (3.9) that , which suggests that applying the standard SIR model to describe transmissible avian-flu spread may significantly over-predict the extent of an epidemic, particularly if .

Numerical results and asymptotic analysis

Fig. 3 shows the evolution of the population fractions and using parameters , and , noting that the mutation rate is small in the simulations with two distinct cases of bird–human infection rate and . We observe that there is an initial time period (), in which when the majority of the population is unaffected by the disease, whereas for the susceptible bird-handlers, , are steadily acquiring the bird–human transmissible form of the disease. At around the time indicated by the vertical line (), the epidemic of the mutated transmissible form begins affecting the whole population. One of the interesting features in these simulations is this apparent time lag before the onset of the epidemic, moreover the occurrence of the bird–human transmissible form appears to be negligible during the epidemic. We note increasing moves the line to the left, whereby for there is no delay in the emergence of the mutated flu epidemic.

Fig. 3

The evolution of and where (a) and (b) with , , and . Note that for case (b) is negligibly small throughout the entire time.

We investigated the model in two biologically relevant cases using singular perturbation theory, namely (1) and and (2) and . Since case (1) is most relevant for avian flu, whereby the rates of bird–human disease transmission and mutation are very low, we present the analysis in some detail in Section 4.1, whilst that for case (2) is briefly summarised in Section 4.2. The method used provides a systematic process by which the solutions of the model are separated as a sequence of timescales, within which the dominant processes are identified and represented by reduced versions of the equations. The aim is to determine, in terms of the model’s parameters, key quantities that characterise an epidemic, namely

the epidemic’s basic reproduction number , as defined just before Section 2.1.

the timescale for the onset of the epidemic following the first cases of the bird–human transmissible form of the disease.

the maximum level of infected individuals at one time .

the fraction of the population to have contracted the disease given by the final size of the removed class .

the duration of the epidemic , defined to be the time in which the fraction of population who are infected by the mutated disease, , is above a threshold value .

The analysis of Section 3 goes some way to answer the fourth quantity, however, the analysis below will give good approximations to the model predictions for all five quantities in the cases of , as demonstrated in the next section, and and .

The case and

As demonstrated in the numerical solutions shown in Fig. 3, the low rates of infection and mutation will result in an initial phase of low epidemiological activity, before the epidemic erupts and dies away. These phases will be made explicit by the analysis below. The first phase involves the equilibration in population of those infected by the bird–human transmissible form of avian flu over an timescale. We will focus mainly on the leading order solutions of the variables in each timescale, which enables and to be independently small with the formal restriction of . This restriction states that the timescale for a mutation event () is much less than the exponential of the timescale for cross-species transmission and arises from the time at which the solutions of the first timescale breakdown (see Section 4.1.1). The analysis for the case follows that of the and limit summarised in Section 4.2 and the distinguishing features of these two limiting cases merge when (analysed in Appendix A to determine an expression for ). We note, in order to distinguish variables for particular timescales, different superscript symbols are used.

This timescale captures the initial phase of avian flu contraction from birds, though the fraction that acquires the disease is very small (size ). Cases of the mutant form are relatively negligible (size ). Denoting all variables in this timescale with a superscript “”, the appropriate rescalings are which leads to the system with , subject to and .

Formally, we write each dependent variable as a power series expansion, e.g. , substitute them into the system of differential equations, equate the terms with the same powers of and and solve the resulting equations. Doing this we obtain the following leading approximations, up to where . The expansions show that initially the susceptible populations are mostly unaffected by both flu strains, with a fraction bird-handlers and non-bird-handlers infected. We note that the carriers of the human–human transmissible form of the disease will grow exponentially if ; this indicates that in the combined limits of and . For , these expansions are valid whilst the incidence rate of catching the mutated form of the disease is significantly less than that of the bird–human transmissible form. Under the restriction of , these solutions breakdown when in Eq. (4.1), whereby the rates of obtaining the two forms become approximately the same. More precisely, in large time , for constant , implies that breakdown occurs at as , leading to a new timescale of events in which human–human transmission begins to be the dominant means of propagating the mutated form of the disease. We note that if the solutions breakdown at ) as at leading order, so the population noticeably drops before the human–human transmission takes noticeable effect (we will not discuss this case further). The matching conditions for the next timescale’s solutions are as .

For the relatively trivial case of , in large time the level of both forms of the disease reaches a maximum, namely and (i.e.  and constant). In a timescale of the balance in the expansions breaks down leading to the second timescale of importance, in which and . Eventually every member of bird-handler class will be infected by the bird–bird transmissible form, but never sufficiently for the mutated form to ever take hold. Here, and as .

The case will not be discussed further and we will now focus on the .

This timescale represents the phase in which contraction of the human–human transmissible form becomes the dominant means of catching the disease. Using the superscript “” to denote the variables in this timescale, we write noting that is now . The rescaled equations are We note in contrast to the previous timescale, that the leading-order contribution to the human–human transmissible class is now from infection of susceptible individuals rather than mutation from . The key leading terms in the expansion of this system, using the matching conditions (4.5) as , are as . In time the spread of the human–human transmissible form of the disease will become dominant and an epidemic will ensue. This will occur when (i.e. ) in a timescale , whereby the transmission rates will become O(1) and the expansions for this timescale are no longer valid. For matching with the next timescale solutions, we have as .

This timescale represents the epidemic phase of the disease in which, for a time, the fraction of infected individuals becomes . Denoting the variables for this timescale with the superscript “”, we write which gives We note that in this timescale decouples from the system at leading order. Defining , then the leading terms and of the series expansions for and , respectively, satisfy the classic SIR model, subject to the conditions and as , using (4.10). Solving this system we get Since implies , then the leading-order maximum extent of the epidemic is given by . Using the same approach, as described in Appendix B, the following estimate for the duration of the epidemic can be determined as , recalling that the epidemic is “defined” to be when . As , the epidemic declines and , then tends to a constant, say, given by Since the human–human transmissible form of the disease does not distinguish between the two susceptible classes, it follows that and as .

Fig. 4 compares the numerical solutions for and (solid lines) with the leading order asymptotic approximations and (dashed lines). Here, Eq. (4.15) is solved numerically on substitution of (4.17) for , using initial conditions derived by matching with (4.10) as , namely at we impose for suitably large ( in the figures). As to be expected there is excellent agreement for (fig. (a)), since both “small’ parameters are very small, however, agreement continues to be good well beyond the expected range of validity as shown for (fig. (b)). The predicted epidemic onset time (dotted line in fig. (a)) also agrees well with solutions of the full problem.

Fig. 4

Comparison of numerical (solid) and the asymptotic approximations (dashed) presented in Section 4.1.3 for (a) and (b) . The dotted line in plot (a) shows . The other parameters are and .

The analysis is not quite complete as the balance on thistimescale breaks down when the epidemic declines to the point when . Here, the infection rates of the class by both forms of the disease are once again in balance; the appropriate system for this fourth timescale is the same as (4.6)–(4.9). The timescale at which depends on as well as ; for example, in the limit it occurs when . The leading order behaviour has and, the result of note in this new timescale, as . From (4.19), the total fraction of the population that will get either form of the disease satisfies , where is given by and here, is the basic reproduction number, see Section 4.2. It is straightforward to show that , so the fraction that will become infected increases with . Furthermore, (4.20) has the following behaviour in two limits of In the first case, the effects of the human–human transmissible form are second order, with the contribution to the class being from the infection of the by the non-transmissible form. For the highly infectious case, only an exponentially small fraction of the population avoids catching the disease.

Eq. (4.20) is the same as that for in Section 3, i.e. . Thus, in the limit and , the fraction of the population that acquires the disease is , which is in agreement with the right plot of Fig. 2 as .

Summary of the asymptotic analysis

In addition to the limit, we used the same approach to examine the cases (1) and (details not shown) and (2) the distinguished limit as to determine an expression for  (Appendix A). Much of the singular perturbation analysis for these limits follows exactly the same path as the above and we omit for brevity the details for case (1), and present a partially complete analysis for case (2) in Appendix A. In case (1), the non-mutated form of the disease readily infects the exposed portion of the susceptible population (though not relevant to avian flu, it may be applicable to other diseases). A notable feature for this case is the exposed portion of the population (i.e. the class) is infected rapidly by the non-transmissible form in comparison to disease mutation, the effective population that is vulnerable to the mutated form are those that remain (i.e. the class); consequently, the effective reproduction number for this new limit is . This value of and that of derived in Section 4.1.1 is bridged by case (2), where in Appendix A it is determined to be in the limit . Here on , we obtain for and for as expected from the analysis of Section 4.1 and for case (1). Fig. 5 shows a heat map of the converged fraction of susceptibles as for and . As expected to the left of the diagram and drops down as increases passed the solid line. The solid line is the plot of Eq. (4.21) and tracks the edge where remarkably well considering the estimate is only logarithmically accurate. In what follows we will focus our discussion on the two principle limits (i) and (ii) and , and not discuss the distinguish limit further.

Fig. 5

Heat map of in space and the approximate position of (solid line) as given by Eq. (4.21). The other parameters are and .

The asymptotic analysis of the model Eqs. (2.9), (2.10), (2.11), (2.12), (2.13) enables the estimation of other key quantities that characterise epidemics in terms of the parameters using relatively simple formula; these are presented for the two principle limits in Table 1. Both these limits produce results which are qualitatively similar, but the details are slightly different. The main predictions of the analysis are

Table 1

Table summarising estimates of the key characteristic quantities for an epidemic from the human–humantransmissible form of the disease, in terms of the model parameters, as predicted by the model from the analysis in Sections 4.1 and from the equivalent analysis for and . If as then , otherwise is a good approximation (see Appendix B). The function is defined in Eq. (4.20).

Epidemilogical features	ε→0,μ=O(1)	ε→0,μ→0
Basic reproduction number R_0	βs_0	β
Epidemic lag time T_L	1βs_0−1ln(1ε)	1β−1ln(1εμ)
Maximum infected fraction I_max	s_0−1β+1βln(1βs_0)	1−1β+1βln(1β)
Epidemic Duration T_σ (β→∞)	(1+1βs0∗)ln(s0∗σ)	(1+1β)ln(1σ)
Infected Fraction R_∞	F(R_∞,βs_0)=0	F(R_∞,β)=0

The basic reproduction number given by (4.21) differs in the two limits due to the division in population classes of the susceptibles . The analysis shows that this is due to the timescales and which form of the diseases each of the susceptible classes becomes infected by. For , members of the class are infected more by the bird–human transmissible form, and consequently have all moved to the class, to leading order, by the time the epidemic of the human–human transmissible form takes hold. For , nearly the whole population is susceptible at the onset of the epidemic.

There will be a significant time lag from when the bird–human transmissible has reached its apparent peak, before the onset of the epidemic of the human–human transmissible form. This lag phase is lengthened by decreased mutation and infection rates, becoming infinite as . These results are demonstrated in Fig. 6, showing the effect on the evolution of (left) and (right).

Fig. 6

The evolution of and for and and .

During the epidemic, the maximum infected fraction and epidemic duration increase with (see Table 1). This indicates the perhaps unintuitive result that if the bird–human transmissible is less infectious (smaller ) then the epidemic from the human–human transmissible form will be worse.

The leading order system that results during the epidemic timescale, Eqs. (4.15), (4.16), is the classic SIR(Susceptible–Infectious–Removed) model; this is true for and case as well, though is replaced with . The usual discussion of the SIR model involves a finite initial number of infectives, say, whilst in contrast the current system involves an infinitesimally small initial state, i.e. . This means that the extent of the epidemic corresponds to the least severe case of that predicted by the standard SIR model.

The upper and lower bounds derived for in Section 3 correspond to Eq. (4.20) for and cases, respectively, as . This explains why the is nearly equal to the upper bound solutions in the left of Fig. 2.

The first two points offer some insights into how the disease or epidemic could be managed. Though the mutation rate cannot be altered, the infection rate and the fraction of susceptible population could be manipulated by quarantining and vaccination programmes. This is discussed in the next section.

Eradication and control: quarantine and vaccination

The analysis of the previous section showed that having some form of direct or indirect control over the infection rates and and the population fraction of susceptibles , could significantly reduce the impact of the disease. In this section, we discuss simple extensions to the model to investigate the possible outcomes of a range of quarantine and vaccination measures. We note that similar studies were undertaken on the avian flu models using numerical simulation (Chong et al., 2014) and sensitivity analysis (Martcheva, 2014). Here, we are able to derive analytical expressions to assess the efficacy of treatment strategies.

Quarantine

It is clear that quarantining will help reduce the impact of the disease, so the aim here is to quantify its effects. We model quarantining by simply removing members of the infected class at a rate proportional to the fraction of infected individuals. Eqs. (2.1), (2.2) are unchanged by this assumption and Eqs. (2.3), (2.4), (2.5) become where is the quarantining rate constant. Eq. (5.2) implies that the average timescale an infective individual is infectious to others is reduced from to . Non-dimensionalising using the rescalings (2.7) and (2.8) and writing leads to the obvious modification of (5.1)–(5.3), (the hat being dropped in what follows). Using the analytical approaches described in 4.1 and the equivalent for the case, we deduce the following

The basic reproduction number becomes and, for , the epidemic time lag is approximately .

Here, and, for , the time lag is approximately .

In both cases, the total infected fraction and the maximal extent of the disease are reduced; satisfying in (4.20), but with . It is clear that quarantining will help reduce the impact of the disease.

Vaccination

Vaccination essentially moves a proportion of the susceptible population to the removed or recovered class. The timing or application of vaccination, however, may differ depending on the circumstances and the accessibility of the vaccines. Here we briefly model three rather ideal scenarios assuming an “all-or-nothing” mechanism for vaccine protection, making references to the two limit cases of and both . Since those vaccinated are assumed to be all alive they should be distinguished from those in class (recall a proportion of which will have died); we will introduce as the vaccinated population fraction.

A pre-infection vaccination programme

Here we consider the case when a fraction of the population has been vaccinated against the disease by . Assuming the vaccination was not targeted, this has the effect of modifying the initial conditions to and . In both the cases, this will have the effect of reducing the impact of the disease, particularly by reducing the basic reproduction numbers in the and cases to and , respectively, thereby increasing the chance of preventing an epidemic.

Vaccination proportional to susceptible population fraction

For this case, the vaccination programme is applied throughout at a rate proportional to the susceptible fraction, this is incorporated into the non-dimensionalised model through Eqs. (2.9), (2.10), where is the dimensionless vaccination rate constant; the equations for and remain unchanged.

To avoid redoing the analysis of Section 4 to study this problem, we justify the results below using a simple argument. In both the limits analysed, during the lag phase prior to the epidemic; the current modification will lead to . For the vaccination programme to be effective, we need so that or more, hence we may deduce

The programme will be totally effective at preventing an epidemic if , and will lessen its impact if .

Same again, but the critical quantity for is about .

Since , the critical level is a small quantity, which indicates that the lag phase for the onset of the epidemic buys enough time to allow a “low level” vaccination rate to be effective. Similar investigations based on the standard SIR model would predict that needs to be in order to be effective; this could be a significant over-estimate, which may compromise the application of this vaccination approach in managing the disease. We note that there is no change to in this analysis.

Fig. 7 (a) shows the recovered/removed fraction against time to demonstrate how the vaccination rate effects the total proportion of the population that had caught the disease. Here, the relevant case corresponds to , and the values used are for and 0.6. The figure suggests that even at (or ) the vaccination programme significantly reduces the size of the epidemic. The effect of the vaccination rate on the fraction of the population that catches the disease, , is shown in Fig. 7 (b). The analysis predicts that the vaccination programme will be effective when , and the numerical results support this. Interestingly, as the infection rate increases, the removed fraction decreases slightly beyond . The general pattern of the profiles shifting to the right on increasing ends at about , as increases further (unlikely to be relevant in biological terms) the profiles retreat to the left due to the increased effects of the vaccination rate (as for fixed ) outpacing that of the rate of infection.

Fig. 7

Plot of (a) the recovered/removed fraction against time and (b) the total fraction that acquire the disease, , against vaccination rate (where ). The curves, from left to right, in plot (a) correspond to and 0.6 and in plot (b) they are and . The parameters, unless otherwise stated, are and .

Fixed rate vaccination programme

In this case, we consider vaccines that are applied to the susceptible population at a fixed rate , starting at when the first case of the bird–human transmissible form is identified. This leads to the following modification to the dimensionless equations (2.9), (2.10), (2.13), where is the Heaviside step function defined as for and otherwise; this is to ensure that the solutions of and remain non-negative.

Using this vaccination programme, the amount of removed susceptibles is proportional to time so will only have an effect if or larger. For both the asymptotic limits investigated we can make the same conclusions as those in Section 5.2.2, exchanging in the discussion with .

A stochastic model

In the early stages of an epidemic, where only a few individuals have the mutated form of the disease, progression is likely to be very stochastic. For instance, if the rate of catching bird flu and the mutation rate to the transmissible form are both very small, then the population of will be small and individuals with bird flu will most likely be recovered/removed before a mutation event occurs. In such circumstances, a deterministic model may not capture the early stages well and we hence investigate in this section a stochastic form of the model proposed in Section 2. In particular, we will investigate how stochasticity effects features of an epidemic such as the maximum number of individuals infected with the transmissible form, total number of those recovered or removed and the timescales for the onset of the epidemic and its duration.

The proposed stochastic model will take the form of a continuous time Markov process, whereby future states are only dependent on the present one. Using the same variable names as before, the four state variables are and , the variable being derived from (2.6). We will use the term “event” to describe a single transition from one state variable to another. Let be a state-space vector and let be a population vector, then Table 2 shows the transition probabilities , where , for each of the seven events as .

Table 2

Transition probabilities from for each of the seven possible events leading to population as ; each transitional probability has an error.

Event	Transitional probability, p_k	I^′
S_b infected by non-transmittable birdflu	μi_1δt	(i_1−1,i_2,i_3+1,i_4)
I_b mutated to I_m	εi_3δt	(i_1,i_2,i_3−1,i_4+1)
S_b infected by mutated flu	(β∕N)i_1i_4δt	(i_1−1,i_2,i_3,i_4+1)
S_n infected by mutated flu	(β∕N)i_2i_4δt	(i_1,i_2−1,i_3,i_4+1)
I_b infected by mutated flu	(β∕N)i_3i_4δt	(i_1,i_2,i_3−1,i_4+1)
I_b recovery/removal	γi_3δt	(i_1,i_2,i_3−1,i_4)
I_m recovery/removal	γi_4δt	(i_1,i_2,i_3,i_4−1)

We can, in principle, use these transition probabilities to construct master equations and generating functions to formulate ODEs for the moments (Allen, 2010). However, as is usual, master equations are difficult to analyse and the formulated ODEs suffer from the identification of suitably reliable moment closure conditions; see Chalub and Souza (2011), Isham (1991), Krishnarajah et al. (2005) and references therein for examples in epidemiological modelling. We will instead analyse the Markov process computationally.

Comparison between continuum and stochastic model

The stochastic model is simulated using the “the first reaction method” (Toral and Colet, 2014). Here, at each time point a time step, , is randomly calculated from an exponential distribution with mean for each of the seven events. The event that is chosen to occur corresponds to the time step . We then adjust the population according to event and move from and repeat the process. The randomly chosen is computed in the simulations using pseudo-random numbers from a uniform distribution, , converted to one satisfying an exponential distribution via inverse transform sampling (Toral and Colet, 2014). In all simulations below, we scale the transition probability parameters so that and , hence the values used for and are the non-dimensional version of the continuum model.

Fig. 8, Fig. 9 compare the mean time profiles from 1000 stochastic simulations with those of the deterministic model. Here, the mean time profile for variable , for example, is defined as where is the population of those infected with the mutant flu at time in realisation ; the corresponding term for the other variables is similarly defined. In Fig. 8, the time profiles of and are shown for (top) and (bottom) for different populations chosen to demonstrate the varied agreement with the deterministic solutions; here and (i.e.  as ) are fixed. As expected, increasing improves the agreement, whilst increasing the parameters for infection rate from birds and the mutation rate also improves agreement for fixed . The latter observation is further demonstrated in Fig. 9. The profiles for in Fig. 8 (left) show that the point of divergence between solutions occurs when the mutated form appears to take hold and infect the entire susceptible class, where for smaller and/or the drop in mean density is less steep. Here, the time at which reaches epidemic levels is highly varied (see Fig. 10) resulting with a profile for that is more smeared out.

Fig. 8

Comparison of the results from deterministic (solid curves) and stochastic (mean of 1000 realisations) simulations. The populations have been normalised to fractions for direct comparison across different total populations . The left plots compare and and the right ones and , whilst the top and bottom plots are for and , respectively. The populations in the top plots are (dashed curves), (dotted curves) and (dash-dotted curves) and for the bottom they are (dashed curves), (dotted curves) and (dash-dotted curves). The other parameters are and for all cases.

Fig. 9

Comparison of stochastic (mean of 1000 realisations, left) and deterministic (right) simulations showing the evolution and for (solid curves), (dotted curves) and (dashed curves) with , and in all cases.

Fig. 10

Comparison of the deterministic model solution with 4 realisations of the stochastic model for (left) and (right). Parameters used are , , and .

Fig. 10 shows comparisons between the deterministic solutions to four instances of the stochastic simulations, using a parameter set that leads to a highly varied epidemic onset time. On the right of Fig. 10 the profiles of the mutated fraction are all very similar, albeit at different time-points at which they reach their peak. Where reaches their peak, the corresponding profile drops down to near zero as shown on the left side of the figure. We note that the deterministic model’s prediction of the rising profile of is in good agreement with each of the stochastic runs, where in the latter stochasticity leads to noise about the steady fraction of . The plots reveal that stochasticity makes little difference with regard to the maximum and duration of infection, this being verified in the results shown in Table 3. Furthermore, the table shows that total number of individuals in the recovered/removed class as is in very good agreement. We note that the standard deviations of the stochastic model are less than 1% of the mean in the first three quantities listed. However, where the deterministic model is a poor predictor of stochastic solutions is with regard to the onset of the epidemic, where for smaller values of and it is a severe under-estimate. For this reason, we analyse the stochastic model further to determine analytical estimates for the epidemic onset time in Section 6.2. We investigated the distribution of onset times for a number of values for , and ; two examples of which are presented in Fig. 11, where the probability density functions for are approximated from 10 000 realisations in each case. Here, the onset times is split into equal time intervals (or bins) and the frequency of elements from the raw data lying in each bin is normalised so that there is unit area under the curve; the time points in the plot are the mid values in the bins. For larger and , the shape of the distribution is near bell-shaped, which broadens and skews to the left as and gets smaller. For the case of (fig. (a)), the profile happens to fit very well to that of the Gumbel distribution (here, the distribution’s two parameters are formulated using the mean and standard deviation of the data). However, on comparing all of the probability distributions available in the mathematics package MATLAB with all of the results generated from the stochastic simulations, no single distribution fits the data well in all cases (indeed all these distributions poorly fitted that shown in the right plot); hence this was not pursued further. Fig. 11 (b) shows the case , where, interestingly, the modal value of in the plot is much closer to that predicted by the deterministic model.

Table 3

Table showing comparisons between the deterministic solutions and the means of 10 000 realisations of the stochastic model, with one standard deviation shown. Here, is defined so that where ; other parameters are , and . All realisations had an epidemic phase.

	Deterministic	Stochastic
Maximum I_m	4.0333×10^6	4.0237×10^6  ±  1.0146×10^4
Duration of epidemic, I_m>σN	9.8698	9.8742  ± 0.084654
Total infected, R_∞	9.8413×10^6	9.8404×10^6±1.0129×10^3
Epidemic lag time T_L	5.2031	70.548  ± 66.949

Fig. 11

Approximations of the probability density function for epidemic onset time, , generated from 10000 realisations of the stochastic model. Plot (a) shows the case of (solid, bin width of 0.1) together with Gumbel distribution (dashed) formulated from the mean and standard deviation of the realisations. Plot (b) shows the case of , plotted with a bin width of 4. The vertical dotted line indicates the mean. The remaining parameters are , and .

Expected time for the epidemic onset

Described in Appendix C are two formulations to estimate the expected time of the epidemic onset . Both formulations assume that during the progression of the epidemic when , the infected populations are negligible in comparison to the uninfected classes, hence is assumed. The first formulation, described in Appendix C.1, assumes that the rise in and occur independently and the estimate of is made by summing the time for to reach a pseudo-steady state , giving say and the time for to reach given , say ; hence . This can be calculated easily, but is a flawed estimate since the rising of and are not independent. The second formulation considers changes in both and leading to a two-variable recurrence relation to obtain . This formulation, if solved directly, is very computationally intensive, particularly for large and less small and ; described in Appendix C.2 is a method to obtain a good estimate of the solution more efficiently.

The bullets in Fig. 12 show the mean epidemic onset time from 1000 trials of stochastic model simulations at points between to for . This shows the rapid increase in as parameters decrease. The dotted line shows that the deterministic model is a very poor approximation for small , but becomes very good beyond . Conversely, the first formulation, curve (b), estimates well for very small , but fails to be a good prediction beyond ; this suggests that assuming independence of and is reasonable when are very small. However, the estimate from the second formulation, curve (a), is in very good agreement across the entire range, suggesting that the constant uninfected population assumption is reasonable. We note, that by increasing improves the deterministic model’s prediction of , as expected; for example, using and (currently unrealistic) , the stochastic simulation mean from a 1000 trials gives is close to the deterministic model’s estimate of .

Fig. 12

Comparison of estimates for the onset of an epidemic defined to be when with and , as a function of . The plot shows the mean onset time from 1000 trials of the stochastic model (bullets), that predicted by the deterministic model (dotted line), and the two estimates of described in Appendix C, namely from Appendix C.1 (dot-dashed line, curve (b)) and described in Appendix C.2 (solid line, curve (a)). The other parameters are , and (see Appendix C).

Conclusion

In this paper, a simple mathematical model is constructed to study the outbreak of an infectious disease following mutation of a non-transmissible strain for which only a subset of the population are at risk of catching. This scenario may arise in avian flu, in which currently only people working with poultry are likely to catch the disease, but a mutation could lead to a human–human spreading strain, with potentially devastating results. We derived a mathematical model extending an SIR type approach that divides each of the susceptible and infectious classes into two subclasses. The susceptibles were divided into those who are (e.g. poultry farmers) or are not at risk of catching bird–human strain and those who are infected with the bird–human and human–human transmissible strains. We note, the assumed division of susceptibles into two (or more) classes has very wide application, where sub-populations with distinct disease exposure characteristics to the non-mutated form could be based on activity, geography or demography. The model has relatively few parameters and four dimensionless ones ( and ), of these could be obtainable using suitable census data and using ordinary flu for the infection rates we may expect to be in the range  Chowell et al., 2006, Lin et al., 2003. Estimation of and will be more difficult and will require clinical data and perhaps laboratory data for estimates of mutation rates.

Analysis of the model was aimed at predicting the severity of an epidemic in terms of the model’s parameter values. In particular, we sought to establish the necessary conditions for an epidemic to occur and key characteristics such as the timescale for the onset of the epidemic, the size and duration of the epidemic and the maximum level of the infected individuals at one time. For general parameter values, we derived in Section 3 bounds on the predicted total fraction of the population who will catch the disease. Further progress was made in two biologically feasible asymptotic limits in Section 4 using singular perturbation theory. Here, analytical approximations to the model solutions were made by solving systematically deduced reduced systems for different timescales, making it easier to identify the important processes at various stages of the epidemic and to express key epidemic features in terms of simple formulae (as summarised in Table 1). In particular, it was demonstrated that the naive assertion of the reproduction number (dimensionless) would not be appropriate for small and or larger (Eq. (4.21)). Furthermore, the analysis highlighted the role of parameters and on the logarithmic dependence on the timescale for epidemic onset, , whilst other aspects of an epidemic (duration, maximum infectives etc.) are largely governed by an SIR model and parameters and . However, it is in the prediction of , for very small and and insufficiently large populations, where the deterministic model fails, as demonstrated in Section 6 using a stochastic model; the rarity of mutation events combined with failure to pass on the infection result in extended onset times. This highlights the need for caution when using a deterministic model to describe the progress of a disease in “large populations”, as the presence of small parameters may necessitate the population to be much larger than expected (perhaps beyond that which is applicable) for the model to represent well the mean of a stochastic process. However, in all other key features of an epidemic, the results of the deterministic and stochastic models were in very good agreement.

In Section 5, the model is extended to investigate the effect of a range of quarantine and vaccination programmes. The simple formulations resulting from asymptotic analysis made it a relatively simple task to predict the outcome of such measures and, in particular, indicate how intensively the measures need to be applied in order for them to be effective. They target two characteristics of the disease predicted by the model: (1) by reducing the effective basic reproduction number (specifically in the cases of quarantine and pre-infection vaccination programmes) and (2) by exploiting the time lag prior to the epidemic, whereby a significant proportion of the population will be vaccinated before the outbreak. Of course, the situation is more complicated than that assumed in the model, with the existence of many and newly evolving flu virus strains and the constant need to develop new vaccines. Without suitable estimates of all the parameters, it is difficult to assert which programme will be the most effective. However, the analysis confirms what is surely well known, that all measures will help and should be implemented.

The infection–mutation pathway assumed in the formulation of the models is intentionally simple for analytical tractability, enabling a number of insights to be drawn in the possible dynamics of epidemics resulting from a mutation; though, for example, the asymptotic analysis of Section 4 can in principle be applied to much larger systems. Of course, such modelling simplifications are made at a cost of a number of potentially important details. The model can be refined to consider other phases (e.g. incubation phase), population classes (e.g. bird populations), endemic scenario (population growth and migration), time dependent infectiousness and recovery (e.g. integro- and/or delay-differential equation systems), and spatial effects (either as a continuum or in a patch network). These refinements will doubtless lead to models with better predictive capabilities in more specific scenarios. Nevertheless, the current model serves as a simple approximation to describe the emergence of mutated, spreadable disease, that is simple enough to analyse in some detail and so provide insight into the timing and appropriate use of control measures should the threat of an epidemic becomes imminent; for avian flu, we hope this will never be the case.

Table 4

Possible state changes of with associated transitional probabilities. They correspond to those shown in Table 2, except that .

(I_b,I_m) at time t	(I_b,I_m) at time t+δt	Transitional probability
(i,j)	(i,j+1)	p_1=βjδt
(i,j)	(i+1,j)	p_2=μ(1−s_0)Nδt
(i,j)	(i,j−1)	p_3=γjδt
(i,j)	(i−1,j)	p_4=γiδt
(i,j)	(i−1,j+1)	p_5=((β∕N)ij+εi)δt
(i,j)	(i,j)	1−∑k=15p_k