|Home | About | Journals | Submit | Contact Us | Français|
In this paper, we demonstrate the uses of some simple mathematical models for the study of disease dynamics in a pandemic situation with a focus on influenza. These models are employed to evaluate the effectiveness of various control programs via vaccination and antiviral treatment. We use susceptible-, infectious-, recovered-type epidemic models consisting of ordinary differential equations. These models allow us to derive threshold conditions that can be used to assess the effectiveness of vaccine and drug use and to determine disease outcomes. Simulations are helpful for examining the potential consequences of control options under different scenarios. Particularly, results from models with constant parameters and models with time-dependent parameter functions are compared, demonstrating the significant differences in model outcomes. Results suggest that the effectiveness of vaccination and drug treatment can be very sensitive to factors including the time of introduction of the pathogen into the population, the beginning time of control programs, and the levels of control measures. More importantly, in some cases, the benefits of vaccination and antiviral use might be significantly compromised if these control programs are not designed appropriately. Mathematical models can be very useful for understanding the effects of various factors on the spread and control of infectious diseases. Particularly, the models can help identify potential adverse effects of vaccination and drug treatment in the case of pandemic influenza.
Mathematical models may be helpful for studying the transmission and control of infectious diseases. This area of study has been highlighted recently due to the advent of a new strain of pandemic influenza, worries about the possible reintroduction of smallpox and the emergence of new pathogens such as the SARS coronavirus. The mathematical theory pioneered by Ross, MacDonald, Kermack, McKendrick, and others has played a major role in the prevention and control of infectious diseases (see also (1)). More recently, researchers have used mathematical models to better understand influenza dynamics and to investigate how to effectively control influenza via various strategies including antiviral treatment and vaccination (see, for example, (2–6)).
We have previously used a simple mathematical model to forecast the course of the 2009 influenza A (H1N1) pandemic in the USA (see (7)). In that study, we fit the model to the Centers for Disease Control and Prevention (CDC) influenza data collected during early months of the year, optimizing the model parameters associated with the time-dependent transmission, and the time at which the initial case was introduced. The optimized model was then used to forecast the timing of the autumn wave of infection (see (8)). The most striking feature of the model is that it accurately predicted the peak time of the pandemic. According to the CDC 2009 H1N1 surveillance reports (9), the peak of the fall wave was reached at the end of October (which is between weeks 42 and 43, see the left-hand plot in Fig. 1), which is consistent with our model result (see the right-hand plot in Fig. 1).
In this paper, we present some simple mathematical models that can be used to study the impact of vaccination and antiviral drug treatment on the spread and control of infectious diseases such as influenza. We focus on models consisting of differential equations that are of the susceptible, infectious, recovered (SIR) type with time-dependent functions representing vaccination and treatment policies. A seasonally forced disease transmission rate is also included to reflect the fact that the transmission rate may be higher in some seasons than others.
For the case of a constant transmission rate, the basic reproduction number and the control reproduction number have been shown to determine the prevalence of infection and the severity of disease outbreaks (1). The formula is a function of the levels of vaccination and antiviral use, and reveals that these control measures are always beneficial in mitigating the disease outbreak (see “The Case of Constant Parameters” section).
However, for the case of periodic transmission rate, our simulation results show that whether and how much these control measures are helpful depends on many factors, including the timing and intensity of the programs, as well as the time of introduction of infections. Particularly, our results suggest that in some cases, an increased level of vaccination or antiviral use may actually lead to higher epidemic peaks and/or larger final epidemic sizes (see results presented in “The Case of Time-Dependent Parameters” section).
This paper is organized as follows; in the “Model Descriptions” section, we introduce the SIR models used for our analysis and simulations. The “Results from Model Analysis and Simulations” section describes the main results and their implications for public health. Comparisons of effects of various vaccination and treatment programs on disease control are also presented in the “Results from Model Analysis and Simulations” section. The “Discussion” section includes a summary and further discussion of the results.
One of the simplest epidemiological models is the SIR epidemic model, in which the total population is divided into three epidemiological classes: susceptible (S), infected (and infectious, I), and recovered (R) individuals. Let S=S(t), I=I(t), and R=R(t) denote the numbers of individuals at time t in the corresponding classes, and let N=S+I+R denote the total population size. A deterministic epidemic model using ordinary differential equations consists of the following equations
with initial conditions
where, I0>0 is a constant. In the initial conditions (2.2), t0>0 denotes the time of introduction of the infection into the population, and it will be shown later that t0 is a critical parameter of the model when the transmission rate β(t) is a periodic function due to seasonality. Note that t=0 corresponds to the first day of a calendar year and t0 is the number of days from t=0. The model (2.1) ignores vital dynamics (births and deaths), which is a reasonable assumption when modeling a pandemic.
If vaccines are available before the epidemic starts, a certain level of population immunity can be achieved via vaccination. Let p0 denote the level of population immunity at time t0. Then the disease dynamics can be modeled by the equations in (2.1) with modified initial conditions:
We do not include vaccinated individuals in the R class because the value of R(t) at the end of the disease outbreak will be used to measure the final epidemic size. In system (2.1), a standard incidence form is used for new infections and the function β(t) represents the rate at which a susceptible individual becomes infected when contacting an infectious individual. The duration of infection is assumed to follow an exponential distribution with the mean period 1/γ; and thus, γ is the per capita recovery rate.
In most deterministic epidemic models, β(t)=β0 is assumed to be constant. However, for many diseases, seasonal variation in the transmission rate can be important. Any periodic function can be expressed as an infinite sum of sines and cosines. In this analysis, we express the periodic transmission rate as the first-order harmonic with a 1-year period:
where, β0 is a constant representing a background transmission rate and ε is a constant related to the magnitude of the seasonal fluctuation. The function β(t) given in (2.4) has its maximum at the beginning of a calendar year. When the transmission rate is expressed in this form, the time of introduction of the pathogen to the population t0 is a crucial parameter of the seasonal model because the final size and shape of the epidemic curve (including how many peaks the epidemic exhibits within one calendar year) can depend quite strongly on this parameter (whereas for models with constant β, the final size and shape are independent of t0). Further, the shape of the epidemic curve also strongly depends on the initial fraction of susceptibles in the population at t0 when β(t) is periodic.
All the variables and parameters as well as their definitions are listed in Table I.
The model (2.1) can be extended to include both vaccination and antiviral drug treatment. Let f0 denote the fraction of infected individuals who will receive treatment at time t; Iu and Itr denote the numbers of untreated and treated infected individuals, respectively; and let the infectious period for an untreated individual is 1/γu. Assume that treatment reduces infectiousness by a factor σ and reduces the infectious period to 1/γtr<1/γu. Following the approach of Lipsitch et al. (6) to model drug treatment, we can extend the model (2.1) as
with initial conditions
where Iu0 is a positive constant representing the initial number of infected people. A transition diagram between the epidemiological classes is shown in Fig. 2.
Where there is a delay in the supply of vaccines, we could still use equations in (2.5) with modifications. For example, let tv>t0 denote the time at which the vaccination program starts and assume that a fraction p0 of the remaining susceptibles will be vaccinated at the time tv. We can first use equations in (2.5) and the initial conditions (2.6) with p0=0 for simulations in the interval . For the time after tv, we can continue to use the equations in (2.5), but with new initial conditions (1−p0)S(tv), I(tv), and R(tv).
In the next section, we will present some analysis and simulations of the models discussed above and then use the results to examine the effects of vaccination and treatment on disease mitigation and control.
In this section, we study models (2.1) and (2.5) and present results under two scenarios: case 1 is for constant rates and case 2 is for periodic transmission rate and time-dependent vaccination and treatment.
In epidemiological models, one of the most important quantities is the basic reproduction number denoted by . The definition of for models with constant parameters is the average number of secondary infections produced by one infected individual during his/her entire period of infection in a wholly susceptible population. This number determines whether there will be an outbreak of the disease when an infectious person is introduced into the population and how severe the outbreak will be in the absence of any programs of disease control and prevention.
Figure 3 illustrates the disease outcomes of system (2.1) with the initial condition (2.2) and constant transmission rate β0. It demonstrates that the threshold condition determines whether there will be an outbreak or not . It also shows clearly that the severity of an epidemic depends on the magnitude of the reproduction number .
Besides the kind of deterministic outcomes presented in Fig. 3, the modeling structure shown in (2.1) can also be used to simulate the situations in which the disease events (e.g., infection and recovery) may occur randomly (while being governed by the event rates determined by the equations in (2.1)). By incorporating the randomness of disease transmission process, the models will be capable of providing helpful information about how likely a certain disease outcome will be expected under a given condition. One of the approaches for stochastic simulations is the event/time approach (10). Consider the model (2.1) with constant transmission rate β. There are only two possible events at time ti (i=0, 1, 2, ···); a new infection which occurs at the rate r1(ti)=βS(ti)I(ti)/N and a recovery which occurs at the rate r2(ti)=γI(ti). Let denote the total event rate at time ti. Then the “time to next event”, Ti, can be determined by
where, Y1 is a random number generated from the exponential distribution with parameter rT(ti). To determine which of the two events occurs at time , we use another random number, Y2, generated from a uniform distribution on [0, 1]. For example, if
then the first event (infection) occurs, in which case
Some of these results using this approach are illustrated in Fig. 4. It shows 20 realizations (repeated runs with the same set of parameters), which are represented by different colors in the top two figures. For demonstration purposes, the parameter values used are β=0.45 and γ=1/3, (for which ) and I0=5 with N=20,000. Effects of model parameters (e.g., ) on the disease outcomes can be examined by computing the statistics generated from the stochastic simulations. For example, among the 20 realizations, each of which showing different peak and final sizes and the time at the peak, the means are: epidemic peak=34% of the population, peak time=50 days, and the final size=40% of the population. The histograms show the variations of the 20 epidemic peak sizes and peak times around the respective means.
We now consider the model (2.5). When control programs are implemented, the corresponding reproduction number is commonly called the control reproduction number and denoted by (c for control). One objective of control programs is to reduce to below one, and the best strategy is usually the one that reduces the most and the fastest for given resources. In the case when β(t)=β0 is a constant for all t≥t0, and when vaccines are available at t0, the reproduction numbers for the model (2.5) are given by
It is clear that in the absence of vaccination and treatment, i.e., when p0=f0=0. For p0>0 and/or f0>0, since σ≤1 and γtr≥γu,
In fact, is a decreasing function of p0 and f0 as shown in Fig. 5. An explicit formula for such as the one given in (3.1) can be very helpful for designing control programs. For example, if the objective of a control measure is to reduce to below a presubscribed value , then we can use this formula to determine the combined efforts of vaccination and treatment needed to achieve the goal . The levels of vaccination p0 and treatment f0 which will lead to are identified in Fig. 5. Several curves corresponding to different values are labeled with the thicker curve corresponding to . We observe that all points (p0, f0) above the curve will lead to , in which case the disease cannot take off.
Note that the parameter p0 in the function denotes the proportion of successfully vaccinated people. If the efficacy of vaccine is not 100%, then where η and denote respectively the efficacy of vaccine and the proportion of vaccinated population. In this case, the formula for can also be considered as a function of η, and f0 given by
This formula allows us to evaluate the effect of improved vaccines. For example, in the absence of treatment (i.e., f0=0) the control reproduction number is a function of η and only. This relation is plotted in Fig. 6. The curve corresponding to the threshold condition, , is identified in the contour plot. For a given proportion , we can use this curve to determine the minimum efficacy ηmin such that for all η>ηmin. The contour curves can also help us identify optimal combinations of η and to achieve certain objectives. For example, suppose that the costs of vaccination and vaccine improvement are known and the total cost to successfully vaccinate a proportion of the susceptibles can be determined, denoted by. Then for a presubscribed value , we can determine the optimal levels ηopt and that minimizes the total cost under the constraint
When model parameters vary with time, analytical results such as the threshold conditions described in (3.1) will be very difficult to obtain. Most results are based on numerical simulations. In this section, we present examples in which the transmission rate β(t) has the form as in (2.4), with possible delayed vaccination uses, i.e., tv≥t0. These scenarios are based on the consideration that vaccines and antiviral drugs may not be available at the beginning of an epidemic. We will again focus on three important measures when assessing the effect of a control program: (a) peak size of the epidemic curve (the maximum number of infections during the course of a pandemic), (b) peak time (the time at which the peak occurs), and (c) final size (the total number of infections at the end of a pandemic). Main objectives of effective control strategies should include lowering the peak size to keep demand for facilities below available supply, lowering the final size to reduce morbidity, and delaying the peak to provide more time for response.
Unlike in the case of constant parameters, where vaccination and treatment will always help reduce morbidity (see Fig. 5), the case of periodic transmission rate β(t) may generate non-intuitive results. That is, the model can exhibit outcomes in which an increased use of vaccination or antiviral drugs will lead to a higher morbidity. To demonstrate this, we first consider the simpler case in which the model (2.1) with the initial conditions given in (2.3) is used. This represents the case when vaccination starts immediately when new infections are introduced into the population. One such example is presented in Fig. 7.
In Fig. 7, both the epidemic curve and curve representing the cumulative infections are plotted. For the transmission function β(t), the parameter values used are those of an H1N1-like disease with β0=0.5 and ε=0.35. We assume 1/γ=3 days. These values are also used in other figures except when specified differently. The figure shows two sets of simulations corresponding to two different times (t0) of initial introduction of pathogens. One case is for t0=30 and the other is for t0=40, and for demonstration purposes the vaccination level is chosen to be p0=0.1.
Some interesting observations can be made from Fig. 7. We can first compare the outcomes between the model without vaccination (the solid curves) and the model with vaccination (the dashed curves). It demonstrates the dependence of the model outcomes on the time of introduction (t0) of initial infections. Particularly, we observe that although the peak and final sizes are both decreased for the case of t0=30, the case of t0=40 differs significantly. Although the first wave is reduced, the second wave is also generated. Moreover, the final epidemic size is even higher than that without vaccination, demonstrating a scenario in which the use of vaccine may have a detrimental consequence.
In the next example, we use the model (2.1) with initial conditions given in (2.3) to examine the interplay between t0, tv, and p0, as well as their influence on the effect of vaccination programs. In Fig. 8, we consider two t0 values; t0=30 in panel panelaa and t0=40 in panel panelb.b. In each of these panels, model outcomes corresponding to different levels of vaccination (p0) are compared. When t0=30, panel panelaa shows that the effect of vaccination use on the peak and final sizes are not monotone. That is, the peak and/or final sizes may either decrease or increase with p0. When t0 is increased from 30 to 40, panel panelbb shows that, although the final sizes decreases monotonically when p0 is increased from 0.1 to 0.25, the final epidemic size for the case of p0=0.1 is actually higher than that without vaccination use.
We now examine the simulation results from the model (2.5), which includes also drug treatment, with possible delays in vaccination use. Simulations illustrated in Figs. 9 and and1010 show how the peak and final epidemic sizes may vary with p0 (the levels of vaccination), f0 (drug treatment), and tv (the starting time of vaccination). Figure 9 plots the epidemic curves and final epidemic sizes for the case of t0=35 under various levels of vaccination (p0) without treatment (f0=0). Notice that for t0=35, there is only a single peak in the epidemic curve in the absence of vaccination and treatment (see (seea).a). Simulation results for two different values tv are demonstrated in the middle panel (tv=60) and the bottom panel (tv=100). In each panel, outcomes for three p0 values are compared. A dramatic difference is that for the case of tv=60 the peak and final epidemic sizes exceed those without vaccination use, whereas for the case of tv=100 the peak and final size decrease monotonically with increasing p0. We need to point out that although this set of simulations suggest that an earlier start of vaccination use is not beneficial (due to an increased peak and final sizes), this may not be the cases for other values of t0 (simulations not shown here). That is, whether and how much vaccination can help to reduce the severity of epidemics depends jointly on tv and t0.
Figure 10 illustrates the effect of drug treatment either used alone (i.e., f0 >0 and p0=0) or combined with vaccination (i.e., f0 >0 and p0 >0). For demonstration purposes, it is chosen that t0=55. Figure 10a is for the case of no control (f0=p0=0). By comparing the panels in the middle (p0=0) and the bottom (p0 >0) we observe that although antiviral use may increase the peak and final sizes when used alone, its benefit can be dramatically improved when combined with the use of vaccination.
Another factor that may also affect the final pandemic size is the initial fraction of infected individuals I0/N at t0. Figure 11 shows the outcomes of the model (2.5) for two different initial values of I0/N. Other parameter values used are β=0.35, γ=1/3, f0=0, p0=0.1, t0=30, and tv=0. Again, it shows that the outcomes can be sensitive to the initial values. Particularly, although vaccination reduced the peak and final sizes (see (seeaa and andc)c) for the case of I0/N=10−7, the use of vaccination actually increased the final size for the case of I0/N=10−8 (see (seebb and anddd).
As mentioned earlier, effective control strategies for pandemics should aim at lowering the peak size to lest demand for facilities exceed supply, lowering the final size to reduce morbidity and delaying the peak to provide time for response. However, results presented here show that an increase in control efforts may not always be beneficial for achieving these objectives. Particularly, the outcomes demonstrated in Figs. 7, ,8,8, ,9,9, ,10,10, and and1111 suggest that in some cases, the uses of vaccine and treatment may have a detrimental consequence (e.g., increasing peak and/or final epidemic sizes), and that whether and how much vaccination and treatment can be helpful will depend critically on the time of introduction of infection (t0), the starting time of vaccination use (tv), and the levels of vaccination and treatment (p0 and f0). Thus, it is very important that these control programs are designed appropriately.
A possible explanation for the negative effects of these control measures is the seasonal variation disease transmission. Vaccination/treatment during the spring wave of an epidemic essentially increases susceptibles at the end of the summer. This larger number of susceptibles can, in some cases, act as tinder for a much larger epidemic to occur in the more favorable influenza season of autumn than may have occurred if treatment had not been used in the summer months. Thus, to increase the likelihood of successful control programs, policy decisions should take into consideration the potential impact of these factors on disease transmission dynamics.
In this paper, we use simple SIR-type epidemic models to study the effect of vaccination and antiviral use on the control and prevention of infectious diseases. Model analyses and simulations are carried out using parameter values related to influenza. The models are simple in the sense that they do not include many of the heterogeneities that may affect disease dynamics. These heterogeneities include age-dependent transmission rate, multiple viral strains (including drug-resistant strains), age or spatial variations in contact patterns, etc. We have previously studied the effects of some of these heterogeneities on the disease dynamics (2,11–14). The results presented here suggest that although vaccination and drug treatment will help reduce the disease prevalence in the case of constant transmission rate, their benefits can be compromised if the transmission rate is seasonally forced.
To demonstrate this through modeling, we compared results of epidemic models with and without vaccination and/or treatment, which helped us to examine how the timing and levels of disease control programs may influence the model outcomes in terms of the peak and final sizes of a pandemic. We demonstrated that when the transmission rate β is constant and when vaccination and treatment are available at the beginning of an epidemic, the disease prevalence always decreases with control efforts by increasing the proportion of population vaccinated and/or treated, and by improving the efficacy of vaccines. The main reason for this is that in this case, the control efforts will reduce the reproduction number , which determines the severity of the disease outbreak (see Figs. 3, ,5,5, and and66).
The situation when β(t) incorporates seasonal variations can differ dramatically. Particularly, the model results suggest that the effectiveness of control programs can be very sensitive to the model parameters representing various control measures. We illustrated several scenarios in which the uses of vaccination or treatment may actually increase peak and final epidemic sizes (see Figs. 8, ,9,9, and and10).10). Although more theoretical studies are needed to better understand how periodic forcing may affect the disease dynamics in general, some preliminary observations can be made based on our simulations.
SIR models with periodic transmission rates (β(t)) are extremely non-linear and the final size of the epidemic is strongly dependent upon the initial conditions such as the time of introduction (t0) of the pathogens into the population, the initial proportion of susceptibles and infectious (S0 and I0, see Fig. 11). Particularly, periodic SIR models can produce epidemic curves with more than one peak; depending on values of mode parameters such as t0 and ε, the spring epidemic may not reduce the number of susceptibles sufficiently to prevent the epidemic from once again taking off in the autumn (during a favorable time of year when β(t) begins to rise again).
In summary, our model studies suggest that the patterns of disease transmission (e.g., with or without seasonal variations) can have a significant influence on the effectiveness of control and intervention programs. When a disease exhibits periodic patterns in transmission, decisions of public health policies will be particularly important as to how control measures such as vaccination and drug treatment should be implemented. More importantly, mathematical models can be very helpful for understanding complex disease transmission dynamics and for identifying the optimal control strategies.
We thank the reviewers for helpful comments and suggestions which greatly improved the presentation of this paper. This research is partially supported by NSF grants DMS-0719697 and DMS-1022758.
An erratum to this article can be found at http://dx.doi.org/10.1208/s12248-011-9305-6