PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Eur J Oper Res. Author manuscript; available in PMC 2012 December 16.
Published in final edited form as:
Eur J Oper Res. 2011 December 16; 215(3): 679–687.
doi:  10.1016/j.ejor.2011.07.016
PMCID: PMC3182455
NIHMSID: NIHMS315337

Generalized Markov Models of Infectious Disease Spread: A Novel Framework for Developing Dynamic Health Policies

Abstract

We propose a class of mathematical models for the transmission of infectious diseases in large populations. This class of models, which generalizes the existing discrete-time Markov chain models of infectious diseases, is compatible with efficient dynamic optimization techniques to assist real-time selection and modification of public health interventions in response to evolving epidemiological situations and changing availability of information and medical resources. While retaining the strength of existing classes of mathematical models in their ability to represent the within-host natural history of disease and between-host transmission dynamics, the proposed models possess two advantages over previous models: (1) these models can be used to generate optimal dynamic health policies for controlling spreads of infectious diseases, and (2) these models are able to approximate the spread of the disease in relatively large populations with a limited state space size and computation time.

Keywords: Infectious Disease Models, Dynamic Health Policy, Discrete-Time Markov Chain, Dynamic Programming, Epidemiology

1 Introduction

The appearance of novel human pathogens (e.g. H1N1 and H5N1 influenza, SARS) and the persistent circulation of infectious diseases in communities (e.g. HIV and tuberculosis), have stimulated efforts to develop dynamic health policies for controlling the spread of infectious diseases. Dynamic health policies make real-time recommendations in response to changing population characteristics (e.g. disease prevalence, proportion of individuals that are immune), disease characteristics (e.g. infectivity, antimicrobial resistance), and resource constraints (e.g. vaccines, antimicrobial drugs, personnel, and budget) (Wallinga et al., 2010; Merl et al., 2009; Ludkovski and Niemi, 2010).

Most existing approaches for identifying optimal policies for infectious disease control use mathematical or simulation models of disease spread as a basis for comparing the performance of a number of pre-determined health policies (Dimitrov et al., 2009; Goldstein et al., 2010; Halloran et al., 2008). Although these approaches allow for projection of the potential impact of different interventions, they are not generally structured to assist dynamic decision making as real-time data on disease spread become available during an epidemic.

In contrast, a few studies have proposed new methods to find exact or approximate optimal dynamic health policies for epidemics. Wallinga et al. (2010) developed a framework for finding the optimal allocation schemes as new observations accrue in the initial phase of an emerging epidemic. These recommendations, however, are dependent on estimates of R0 which poses inherent drawbacks in dynamic decision making in stochastic environments (refer to Larson (2007) for detailed discussion on the limitations of R0). Lefevre (1981) used a continuous-time Markov decision model, Merl et al. (2009) developed a statistical framework and Ludkovski and Niemi (2010) developed a simulation-based model for dynamic determination of optimal policies for emerging epidemics.

Since the transmission of infectious diseases is a stochastic process, optimal dynamic health policies for limiting disease spread can potentially be determined through dynamic programming techniques (Bertsekas, 2005) (or approximate dynamic programming (Powell, 2007)). These methods have proven to be efficient and effective for dynamic decision making in a wide variety of areas (e.g., medical treatment optimization (Schaefer et al., 2005), economics (Van and Dana, 2003), operations research (Winston, 2003; Bertsekas, 2005; Powell, 2007)). Yet, the use of these techniques for assisting the selection of infectious disease control strategies is limited (to the best of our knowledge, Lefevre (1981), Ludkovski and Niemi (2010) and Ge et al. (2010) are the only examples). In part, this may reflect the failure of existing models of infectious diseases to satisfy the requirements of dynamic optimization techniques.

There are two main features of infectious disease transmission processes that pose challenges for the use of dynamic optimization methods:

  1. Prohibitively large state space: The state of an epidemic is usually described by the number of individuals in each disease compartment (e.g. susceptible, infectious, recovered). Therefore, as the population size grows, the size of the state space increases rapidly. For instance, in a closed population of size N, the size of the state space of a simple epidemic susceptible-infectious-removed (SIR) model will be N(N + 1)/2. Such enormous state spaces mean dynamic programming methods to lose their efficiency very rapidly.
  2. Unobservability of state: The within-host natural history of an infectious disease can be complicated and variable; this, coupled with limited availability or access to diagnostic tests results in uncertainty about the true state of the epidemic at any point in time. For example, for infectious diseases with a longer period of incubation than latency, infectious individuals may be asymptomatic and hence unlikely to be diagnosed for a variable amount of time. Although the number of symptomatic infectious individuals may be observable, asymptomic infectious individuals will be hard to detect. This means that the overall state of the epidemic cannot be measured accurately. Therefore, in order to utilize (approximate) dynamic programming, the model of disease spread should be structured such that a probability belief can be formed about the state of disease spread using real-time data.

In this paper, we propose a class of mathematical models that retains the strength of existing modeling approaches (e.g. representation of the within-host natural history of disease and between-host transmission dynamics) while permitting the efficient use of dynamic programming techniques to develop optimal dynamic health policies for disease control. The state space in this class of models can be reduced by state aggregation while maintaining a desired level of accuracy for the disease spread model. Also, the proposed class of models provides a structure for formation of probability beliefs about the actual state of disease spread based on observed data.

In §2, we review a number of available infectious disease models and highlight the manner in which these models fail to meet the requirements of dynamic optimization. In §3, we describe a procedure to construct the proposed class of discrete-time Markov models for infectious disease transmission. Section 4 demonstrates how the proposed framework can be employed in constructing Markov models for susceptible-infectious-susceptible (SIS) and SIR models as well as their derivatives. Section 5 discusses a method to reduce the state space of Markov models of infectious diseases through state aggregation. Finally, using the proposed class of models, we present an illustrative SIR model for an influenza outbreak in an English boarding school.

2 Literature Review

Although infectious diseases generally spread in a stochastic fashion, deterministic models are commonly used as tools for studying epidemic behavior (Anderson and May, 1992; Hethcote, 2000). These deterministic models have been very useful in understanding the dynamics of infectious disease, estimating important epidemiologic parameters (e.g. basic reproductive numbers), and determining targets for disease control (e.g. critical proportions of the population to immunize). However, if the aim of a model is to help develop a dynamic health policy (that is a policy that can recommend switching interventions based on real-time observations about the epidemic state), it is not useful to consider models that produce epidemics with deterministic trajectories. If the epidemic trajectory could be known with certainty, it is possible to determine an optimal series of interventions at baseline and thus the motivation for dynamic decision-making is lost. Therefore, finding optimal dynamic health policies requires the use of stochastic models of infectious disease spread.

Many stochastic models of infectious diseases utilize non-negative integer-valued Markov processes in continuous or discrete time. To date, most stochastic models of infectious diseases have been based on a continuous-time Markov chain (Jacquez and O’Neill, 1991; Jacquez and Simon, 1993; Nåsell, 2002; Keeling and Ross, 2008). In these Markov models, the state of the process is defined as the number of individuals that are susceptible, infected, etc. Therefore, for population of any reasonable size, the number of Kolmogorov’s differential equations describing the disease spread is prohibitively large. For example, in a population of N individuals, to model SIS dynamics by continuous-time Markov chain, N + 1 differential equations are needed, and to model SIR dynamics, (N + 1)(N + 2)/2 differential equations are required (Keeling and Ross, 2008). In addition, dynamic optimization over a continuous-time model can be rather challenging, resulting in policies which are not convenient to implement in practice, and limited to very simple models of infectious diseases (Lefevre, 1981).

Most discrete-time Markov models assume that the time step is sufficiently small so that only one change in state is possible during the time step. A change may be a birth or death of a susceptible or infected individual, recovery of an infected individual, an infection of a susceptible individual, etc. (Allen, 1994; Allen and Burgin, 2000; Castillo-Chavez and Yakubu, 2001). As such, the transition probabilities obtained from these models simply approximate the transition probabilities in a continuous-time Markov jump process (Allen and Burgin, 2000), and cannot be used to build the transition probability matrix of the associated discrete-time Markov chain, which is required for dynamic optimization methods.

There are important historical examples of discrete-time Markov chain models for infectious diseases. Reed-Frost and Greenwood models are probably the best-known discrete-time stochastic epidemic models (Abbey, 1952; Greenwood, 1931). In these Markov models, the state of the disease spread is defined as the number of individuals that are susceptible, infected, etc. Therefore, similar to the case of continuous-time Markov process, as the population size grows, the size of the probability matrixes becomes prohibitively large. In addition to requiring a large state space, these discrete-time Markov chains have generally been used to describe the spread of a pathogen with an infectious period that is relatively short in comparison with the latent period (Daley and Gani, 1999). This assumption is violated for many diseases with complex natural histories.

3 Discrete-Time Markov Models for Spread of an Infectious Disease

In this section, we discuss the construction of a generalized class of discrete-time Markov models of infectious diseases that fulfill the requirements of dynamic programming techniques. We then show how the parameters of the proposed class of models can be estimated and belief states can be formed using the available data.

3.1 Model Construction

To illustrate the steps required to construct discrete-time Markov models of infectious disease spread, we consider a hypothetical infectious disease with a natural history that can be adequately summarized with M serial classes (see Figure 1).

Figure 1
A population with M serial classes

At a given time t, we denote the number of individuals in class Ci by XCi(t), for i [set membership] {1, 2, …, M}. As with any infectious disease model, we must first determine the set of classes that can adequately describe the disease dynamics of interest. Once this set of classes is identified, we can write a dynamics state equation, which specifies the relationship among the numbers of individuals in each class (i.e., XCi(t)’s). This fundamental equation is used to determine the disease states that must be included in the Markov process. If we assume that in Figure 1 the population size is fixed and equal to N, then the dynamics state equation for this disease will be:

i=1MXCi(t)=N.
(1)

By Eq.1, the disease state is fully identified if we know M − 1 variables of {XC1(t), XC2(t), …, XCM(t)}, since the variable excluded is determined by the values of the others.

The state of the system changes as events occur, such as births or deaths of susceptibles, transmission episodes, recoveries or deaths of infectives, etc. We call these events the dynamics driving events. In Figure 1, the driving event from class Ci to Ci+1 is denoted by Ci+1(t), which is a non-negative discrete random variable representing the number of transitions from class Ci to Ci+1during the interval [t, t + Δt]. For example, the driving event from the class susceptible S to infective I, is the random variable I(t) which represents the number of new infections occurring during the interval [t, t + Δt]. Let PCi(t)(·) denote the probability mass function for the random variable Ci(t); that is PCi(t)(c) = Pr{Ci(t) = c}, i [set membership] {2, 3, …, M}.

The random variables C2(t), C3(t), …, CM (t) are assumed to be independently distributed and to be determined only by the state of the disease at time t, i.e., {XC1(t), XC2(t), …, XCM(t)}. This assumption implies that transitions from states occur independently of each other; for example, over a given period, the number of births, the number of new infections and the number of recovered from infection all occur independently. In other words, information about the number of new infections over a given period does not change the probability distribution for the number of births or deaths over the same period.

The second step in constructing the discrete-time Markov model is to find the probability distribution of each of the driving event random variables conditional on the state of the disease at time t: PCi(t)(·|XC1(t), XC2(t), …, XCM(t)). Since the random variables C2(t), C3(t), …, CM (t) are independently distributed, the joint probability mass function of the random variable (C2(t), C3(t), …, CM (t)) will be:

P(C2(t),C3(t),,CM(t))((c2,c3,,cM)XC1(t),XC2(t),,XCM(t))=i=2MPCi(t)(ciXC1(t),,XCM(t)).
(2)

As we will see later, it is usually straightforward to find the probability function of each driving event random variable, and in consequence the joint probability mass function (2).

To construct the discrete-time Markov model, we then calculate the transition probability matrix whose rows represent the probability distribution of the disease state at time t + Δt, given the state of the disease at time t:

Pr{(XC1(t+Δt),XC2(t+Δt),,XCM(t+Δt))=(x1,x2,,xM)XC1(t),XC2(t),,XCM(t)}.
(3)

In calculating the transition probabilities (3), we first need to find a way to relate the joint random variable (XC1(tt), XC2(tt), …, XCM(tt)) to the driving event random variables (C2(t), C3(t), …, CM (t)) whose joint probability mass function was calculated in (2); and second, we need to determine the support of the joint random variable1 (XC1(tt), XC2(tt), …, XCM(tt)) conditional on the disease state at time t, {XC1(t), XC2(t), …, XCM(t)}. Let ΩX(t) denote this support. These two conditions can be accomplished by two sets of constraints: dynamics driving constraints and dynamics feasibility constraints.

The set of dynamics driving constraints summarizes the relationships among the driving events during interval [t, tt] and the state of the disease at time t and tt. For the disease in Figure 1 with a fixed population size, the dynamics driving constraints are:

XC1(t+Δt)=XC1(t)C2(t),
(4)

XCi(t+Δt)=XCi(t)+Ci(t)Ci+1(t),fori{2,3,,M},
(5)

XCM(t+Δt)=XCM(t)+CM(t).
(6)

These equations specify the new values of (XC1(t + Δt), XC2(t + Δt), …, XCM(t + Δt)), if during the period [t, t + Δt] the driving events take value of C2(t), C3(t), …, CM (t). Adding up the equations (4)(6) results in i=1MXCi(t+Δt)=i=1MXCi(t), which must equal N by the dynamics state equation (1).

Therefore, constraint (6) is redundant and can be dropped. By solving the set of equations (4)(5) for Ci(t), i [set membership] {2, 3, …, M 1}, we obtain:

Ci(t)=j=1i1(XCj(t)XCj(t+Δt)),fori{2,3,,M}.
(7)

Given that the random variables Ci(t), i [set membership] {2, 3, …, M} are nonnegative, and cannot exceed the number of individuals in class Ci−1, the set of constraints (7) can now be used to derive the dynamics feasibility constraints, which specify the feasible ranges for the state of the disease at time t + Δt, given {XC1(t), XC2(t), …, XCM(t)}. Assuming a fixed population size, the dynamics feasibility constraints for Figure 1 will be:

0j=1i1(XCj(t)XCj(t+Δt))XCi1(t),fori{2,3,,M}.
(8)

Therefore, in finding the transition probabilities (3), the support of the random variable (XC1(tt), XC2(t + Δt), …, XCM(t + Δt)) conditional on {XC1(t), XC2(t), …, XCM(t)} is determined by the set of dynamics feasibility constraints (8):

ΩX(t)={(x1,x2,,xM)NM0j=1i1(XCj(t)xj)XCi1(t),i{2,3,,M}}
(9)

and the relationship between the random variables (XC1(t + Δt), XC2(t + Δt), …, XCM(t + Δt)) and (C2(t), C3(t), …, CM (t)) are established by the set of dynamics driving constraints (7).

Thus, by using the probability mass function (2) and the set of dynamics driving equations (7), the transition probabilities (3) can be calculated by:

Pr{(XC1(t+Δt),,XCM(t+Δt))=(x1,,xM)XC1(t),,XCM(t)}={P(C2(t),C3(t),,CM(t))((j=11(XCj(t)xj),,j=1M1(XCj(t)xj)XC1(t),,XCM(t)),if(x1,,xM)ΩX(t),0,otherwise.
(10)

To summarize, the discrete-time Markov model can be constructed through the following steps:

  • Step 1: Define the classes and form the dynamics state equation (Eq.1).
  • Step 2: Find the joint probability distribution of the driving events (Eq.2).
  • Step 3: Form the dynamics driving and feasibility constraints (constraints (7) and (8)).
  • Step 4: Calculate the transition probability matrix for the Markov chain {(XC1(t), …, XCM(t)): t = 1, 2, …} (Eq.10).

3.2 Parameter Estimation

Data accrued over the course of an epidemic permits estimation (and updating) of model parameters. In this subsection, we discuss how a maximum likelihood method can be used to estimate the parameters of our proposed class of infectious disease models.

Suppose that we want to estimate the parameters ω = (ω1, ω2ωZ) of the Markov chain {(XC1(t), XC2(t), …, XCM(t)): t = 1, 2, …}. Depending on disease dynamics, a variety of observations may be gathered during each period. Among the M − 1 driving events {C2, C3, …, CM}, let us assume that driving events with index k [set membership] O [subset or is implied by] {2, 3, …, M} can be observed during period [t, t + Δt] (e.g. birth of a susceptible or the number of new infections). To simplify notation, without loss of generality, we assume throughout this section that each period’s length equals a single unit of time, i.e. Δt = 1. Let the random vector D(t) = (Ck(t)), k [set membership] O denote the vector of the observable driving events during period [t, t + 1] whose observed values are stored in vector Dt = (Ĉk(t)), where Ĉk(t) denotes the observed value of the event Ck, k [set membership] O, during period [t, t+1].

Suppose that the initial state of disease spread is observable and denoted by X^(1)=(x11,x21,xM1), and that during T periods {1, 2, …, T} observations D = (D1, D2, …, DT) are gathered. The likelihood of observing D = (D1, D2, …, DT) is then equal to:

L(ω1,,ω2;D^1,D^2,,D^T,X^(1))=Pr{D(1)=D^1,D(2)=D^2,,D(T)=D^Tω,X^(1)}
(11)

Now, the model parameters ω = (ω1, ω2ωZ) can be estimated by maximizing the likelihood function (11) over parameters ω = (ω1, ω2ωZ). The likelihood function (11) is equivalent to:

L(ω1,ω2;D^1,D^2,,D^T,X^(1))=Pr{D(T)=D^TD(1)=D^1,D(2)=D^2,,D(T1)=D^T1;ω,X^(1)}×Pr{D(T1)=D^T1D(1)=D^1,D(2)=D^2,,D(T2)=D^T2;ω,X^(1)}××Pr{D(2)=D^2D(1)=D^1;ω,X^(1)}×Pr{D(1)=D^1ω,X^(1)}.
(12)

In order to calculate the likelihood function (12) we define the following notation:

  • X(t) = (XC1(t), XC2(t), …, XCM(t)): random vector representing the state of the disease spread at time t.
  • πt(x1, x2, … xM; ω, X (1)): probability that at time t the state of the disease spread is X(t) = (x1, x2, … xM) given the model parameters ω = (ω1, ω2ωZ) and the initial state X^(1)=(x11,x21,,xM1).
  • Ot: set of all states that the disease spread may visit at time t given the observations D = (D1, D2, …, Dt−1).

Now given the observations D = (D1, D2, …, Dt−1), the term

Pr{D(t)=D^tD(1)=D^1,D(2)=D^2,,D(t1)=D^t1;ω,X^(1)}

in likelihood function (12) is calculated as follows:

Pr{D(t)=D^tD(1)=D^1,D(2)=D^2,,D(t1)=D^t1;ω,X^(1)}=(x1,x2,xM)Ω^tPr{D(t)=D^tX(t)=(x1,x2,xM)}πt(x1,x2,xM;ω,X^(1))(x1,x2,xM)Ω^tπt(x1,x2,xM;ω,X^(1)).
(13)

It only remains to establish a procedure to update the set Ot as observations accrue. For time t = 1, the set O1 has only one member: the initial state of the epidemic, X^(1)=(x11,x21,xM1). Therefore:

Ω^1={(x1,x2,xM)NM(x1,x2,xM)=(x11,x21,xM1)}.

For time t > 1, the set Ot can be determined recursively by using the set Ot−1 and the observation Dt−1 = (Ĉk(t − 1)), k [set membership] O, gathered during period [t − 1, t]. To this end, we establish how the set of observations Dt−1 = (Ĉk(t − 1)), k [set membership] O, during period [t − 1, t], restricts the possible value of the disease state at time t. This can be achieved by using dynamics driving constraints (7) and dynamic feasibility constraints (8). Let us denote members of set Ot−1 by (y1, y2, …, yM) [set membership] Ot−1 and members of set Ot by (x1, x2, …, xM) [set membership] Ot. The set Ot is now determined by:

Ω^t={(x1,x2,,xM)NM(y1,y2,,yM)Ω^t1,0j=1i1(yjxj)yj1,i{2,3,,M},andC^i(t1)=j=1i1(yjxj),iO}.

We note that for some special cases, the number of individuals in some classes (compartments) may also be observable. For example, if all individuals remain symptomatic throughout their infectious period, the prevalence of disease can be measured at any time. While in this section we assumed that the only observable data on the disease spread come from a set of driving events, the likelihood function (11) and analysis presented above can be modified when a set of disease classes (compartments) is also observable.

3.3 Calculation of Belief States

In the previous subsection, we demonstrated how model parameters can be updated by using the observations D = (D1, D2, …, DT) gathered during the first T periods. To make a decision at decision epoch T + 1, dynamic programming methods require the calculation of a belief state at time T +1, given the model parameters ω = (ω1, ω2ωZ), the observations D = (D1, D2, …, DT), and the initial state of disease spread X(1). This belief state can be calculated as below:

Pr{X(T+1)=(x1,x2,xM);ω,D^,X^(1))={πt(x1,x2,xM;ω,X^(1))(y1,y2,yM)Ω^Tπt(y1,y2,yM;ω,X^(1)),if(x1,x2,xM)Ω^T,0,otherwise,

where πt(x1, x2, … xM; ω, X (1)) is defined in the previous subsection.

4 Discrete-Time Markov Models for SIS, SIR and Their Derivatives

Discrete-time Markov chain models are typically used for pathogens with relatively short (and fixed) durations of infectiousness (Daley and Gani, 1999; Larson, 2007); in these models, the time step (Δt) is usually set to the duration of the infectious period. In §4.1 we develop a discrete-time Markov SIS model incorporating this convention to demonstrate that this model is a special case (a member) of our proposed class of infectious disease models. For the purpose of dynamic decision making, however, we require a time step Δt that represents the interval between the two consecutive decision epochs, which is determined by the decision maker and will depend on the context of the problem being studied. Hence, models for disease transmission that will support dynamic decision-making should be flexible enough to incorporate different Δt, varying over a reasonable range where transition probabilities are still reasonably accurate. Accordingly, in subsequent subsections the time steps Δt can be set by the decision maker.

To denote the classes, we use conventional notation: Upon birth, an individual enters into the susceptible class S. If contact between an susceptible and an infective occurs and results in the transmission of infection, the susceptible moves to exposed class E; individuals in this class are in a latent period and are infected but not yet infectious. When the latent period ends, the individual enters the infective class I and is capable of transmitting the infection. If infection results in permanent immunity, an individual cured from infection enters the recovered class R, otherwise the individual moves back to the susceptible class S.

For each models presented, we assume that the members of population become infected only through contact with other infectious members, and that contacts during the interval [t, t + Δt] occur according to a homogenous Poisson process, with rate λ Delta;t. We initially assume a closed population size of size N. In §4.3, we discuss extensions allowing for birth and death.

4.1 SIS with short, fixed infectious period

SIS models are simple models in which individuals move from the susceptible class to the infective class and then back to the susceptible class upon recovery; there is no immunity conferred by a previous infection (Hethcote, 2008). We assume that the infectious period is fixed at length Δt; that is, an infective at time t remains infectious over the interval [t, tt] but will be diagnosed and effectively treated at the end of this period. Such person will reenter the population as susceptible at time t + Δt. The SIS model consists of two classes: the susceptibles C1 = S, and the infectives C2 = I. Let XS(t) denote the number of susceptibles at time t, XI(t) denote the number of infectives at time t, I(t) denote the number of new infections during the period [t, t + Δt], and S(t) denote the number of infectives who are treated during period [t, t + Δt] and return to the susceptible class by time t + Δt.

Step 1 - dynamics state equation

For a population with fixed size N, the dynamics state equation will be XS(t)+ XI(t) = N, and hence only a single class is needed to define in the Markov model. We select XS(t) as the state of the Markov model.

Step 2 - distribution of the driving events

In Step 2, we find the probability distribution of the driving event I(t) conditional on the state of the disease at time t, i.e., PI(t)(i|XS(t)). To this end, let us first introduce the following notation:

  • r(t): probability that a susceptible person becomes infected upon contact with an infectious individual.
  • β(t): probability that the next interaction of a random susceptible person is with an infectious person.
  • q(t): overall probability that a susceptible person becomes infected.

Both r(t) and β(t) can be decision variables, affected by “hygienic interventions” (reducing the chance of transmission given contact between infectious and susceptible individuals) and “social distancing” (reducing the likelihood of contact between susceptible and infectious individuals), respectively. When social distancing has not been used and mixing is homogenous β(t) is equal to β(t) = 1 − XS(t)/N. To calculate the probability q(t), first note that a random susceptible person will contact n individuals during the interval Δt, where n has Poisson distribution with rate λΔt. The number of infective persons, j, among these n individuals has binomial distribution (n, β(t)); and finally, the susceptible person may become infected if the contact with at least one of these j individuals results in infection (so the probability that the susceptible person becomes infected is one minus the probability that none of the interactions with infectives results in infection). Hence,

q(t)=n=0eλΔt(λΔt)nn!(j=0n(nj)β(t)j(1β(t))nj(1(1r(t))j))=1n=0eλΔt(λΔt)nn!(j=0n(nj)β(t)j(1β(t))nj(1r(t))j).
(14)

The expression j=0n(nj)β(t)j(1β(t))nj(1r(t))j is the z-transform of the binomial distribution (n, β(t)) for z = 1 − r(t), and is equal to [1 − β(t)+β(t)(1 − r(t))]n. Therefore, Eq.14 results in:

q(t)=1n=0eλΔt(λΔt)nn!(1β(t)r(t))n.
(15)

In Eq.15, the expression n=0eλΔt(λΔt)nn!(1β(t)r(t))n is the z-transform of the Poisson distribution with rate λΔt for z = 1 − β(t)r(t), and hence, Eq.15 results in:

q(t)=1eλΔtβ(t)r(t).
(16)

When the probability that a random person from XS(t) susceptible becomes infected is q(t), then the number of new infections during the interval [t, t + Δt] will have binomial distribution (XS(t), q(t)); hence the probability mass function for the driving event I(t) will be:

PI(t)(iXS(t))={(XS(t)i)q(t)i(1q(t))XS(t)i,for0iXS(t),0,otherwise.
(17)

Note that for an SIS model, the state XS(t) = N is an absorbing state since no infectives remain in the population. For this state, the probability function (17) is defined as:

PI(t)(iXS(t)=N)={1,fori=0,0,otherwise.

Next, we find the probability mass function for the driving event S(t). Since by assumption, all the infectives at time t, XI(t) will be diagnosed and treated (and hence become susceptible again) by time t + Δt, the random variable S(t) takes a single value XI(t) = NXS(t):

PS(t)(sXS(t))={1,fors=NXS(t),0,otherwise.
(18)

Step 3 - dynamics driving and feasibility constraints

In Step 3, we form the dynamics driving and feasibility constraints. For an SIS model, the number of susceptibles at time t + Δt, XS(t + Δt), is equal to the number of susceptibles at time t, XS(t), minus the number of new infections during period [t, t + Δt], I(t), plus the number of infectives who are treated and return to the susceptible class by time t, S(t); that is, XS(t + Δt) = XS(t) − I(t) + S(t). Therefore, the dynamics driving and feasibility constraints will be:

I(t)=XS(t)XS(t+Δt)+S(t),and
(19)

0XS(t)XS(t+Δt)+S(t)XS(t).
(20)

Since, by assumption, all the infectives at time t will be removed by time t + Δt, we have S(t) = XI(t) = NXS(t). Substituting S(t) = NXS(t) in constraint (19)–(20) results in2:

I(t)=NXS(t+Δt),and
(21)

0XS(t+Δt)N.
(22)

To calculate the probability support (9), by probability function (17) we have 0 ≤ I(t) ≤ XS(t), which combined with Eq.21 and inequality (22) leads to NXS(t) ≤ XS(t + Δt) ≤ N. Therefore, the probability support (9) will be:

ΩXS(t)={xNNXS(t)xN.}
(23)

Step 4 - transition probabilities of Markov chain {XS(t): t = 1, 2, …}

The transition probability matrix of the Markov chain {XS(t): t = 1, 2, …} is given by (10):

P{XS(t+Δt)=xXS(t)}={PI(t)(NxXS(t)),ifNXS(t)xN,0,otherwise,
(24)

where the probability function PI(t)(·|XS(t)) is given by Eq. 17.

4.2 SIR

If individuals that recover from infections acquire permanent immunity, then an SIR model can be used to represent disease spread (Hethcote, 2008). Here we assume that the infectious period is exponentially distributed with mean duration 1/μI. Let XS(t) denote the number of susceptibles, XI(t) denote the number of infectives and XR(t) denote the number of recovereds at time t. Also let I(t) denote the number of incident infections during the period [t, t + Δt] and R(t) denote the number of recovered from infection during the period [t, t + Δt].

Step 1 - dynamics state equation

In SIR model for a population of fixed size N, the dynamics state equation is: XS(t) + XI(t) + XR(t) = N; hence two classes are sufficient to construct the Markov model. We select (XS(t), XI(t)) as the state of the Markov model.

Step 2 - distribution of the driving events

The probability distribution of the driving event I(t) conditional on the state (XS(t), XI(t)) is the same as the probability distribution (17) for the SIS model of §4.1, except that in calculating q(t), β(t) should be set to β(t) = XI(t)/N:

PI(t)(iXS(t),XI(t))={(XS(t)i)q(t)i(1q(t))XS(t)i,for0iXS(t),0,otherwise.
(25)

To calculate the probability distribution of the driving event R(t) conditional on the state(XS(t), XI(t)), let us first define one additional parameter:

  • ρ(t): probability that an infective at time t recovers by time t + Δt.

Since we assume that the duration of infectiousness is exponentially distributed with a mean length of 1/μI, ρ(t) can be calculated as: ρ(t) = 1 − eμIΔt.

The number of infectives at time t who recover during the period [t, t + Δt] will have binomial distribution with probability of success ρ(t) and the total number of trials XI(t):

PR(t)(rXS(t),XI(t))={(XI(t)r)ρ(t)r(1ρ(t))XI(t)r,for0rXI(t),0,otherwise.
(26)

Note that for an SIR model, the states (XS(t), XI(t)) = (xS, 0), 0 ≤ xSN, are absorbing states since no infectives remain in the population. For these states, the probability function (25) is defined as:

PI(t)(i(XS(t),XI(t))=(xS,0))={1,fori=0,0,otherwise,

and the probability function (26) is defined as:

PR(t)(r(XS(t),XI(t))=(xS,0))={1,forr=0,0,otherwise.

Step 3 - dynamics driving and feasibility constraints

The dynamics driving equations of SIR model can be easily obtained by Eq.7:

I(t)=XS(t)XS(t+Δt),and
(27)

R(t)=XS(t)XS(t+Δt)+XI(t)XI(t+Δt).
(28)

The dynamics feasibility constraints for the SIR model can be obtained by Eq.8:

0XS(t)XS(t+Δt)XS(t),and
(29)

0XS(t)XS(t+Δt)+XI(t)XI(t+Δt)XI(t).
(30)

Constraints (29)–(30) are equivalent to:

0XS(t+Δt)XS(t),andXS(t)XS(t+Δt)+XI(t+Δt)XS(t)+XI(t).

Therefore, the probability support (9) will be:

Ω(XS(t),XI(t))={(xS,xI)N20xSXS(t),XS(t)xS+xIXS(t),XI(t)}.
(31)

Step 4 - transition probabilities of the Markov chain {(XS(t), XI(t)): t = 1, 2, …}

By dynamics driving equations (27)(28) and probability support (31), the transition probability Pr{(XS(t + Δt), XI(t + Δt)) = (xS, xI)|XS(t), XI(t)} will be equal to:

Pr{(XS(t+Δt),XI(t+Δt))=(xS,xI)XS(t),XI(t)}={Pr{(I(t),R(t))=(XS(t)xS,XS(t)+XI(t)xSxI)XS(t),XI(t)},for0xSXS(t),XS(t)xS+xIXS(t)+XI(t),0,otherwise.
(32)

By the assumption that I(t) and R(t) are independent, the probability function (32) results in:

Pr{(XS(t+Δt),XI(t+Δt))=(xS,xI)XS(t),XI(t)}={PI(t)(XS(t)xSXS(t))PR(t)(XS(t)+XI(t)xSxI)XS(t),XI(t)),for0xSXS(t),XS(t)xS+xIXS(t)+XI(t),0,otherwise.
(33)

4.3 Other Infectious Disease Models: Incorporating Classes with Non-Exponential Stay Time, Branching, Birth and Death

The proposed models can be easily extended to capture the spread of diseases with more complex natural histories. For example, for diseases such as tuberculosis that are characterized by latent (non-infectious) periods after infection, models will often also include an exposed class E, to which infected individuals progress prior to becoming infectious. These SEIR models include four serial classes and can be easily constructed through the five steps outlined in §3.

If birth is assumed to occur according to Poison distribution with rate λB, it can be incorporated into the disease dynamics by adding a dummy class at the beginning of class sequence. Death can be modeled using a similar approach, though several death classes may be needed to reflect each of the classes from which individuals die.

Individuals may also be expected to have different rates of exit from particular classes. For example, consider a disease in which αk, k [set membership] {1, …, K} portion of population may remain in the latent period according to exponential distribution with an average duration of 1/μk, where k=1Kαk=1 (See Figure 2).

Figure 2
Branching in Infectious Disease Models

In Figure 2, the joint probability function of the driving events (I1(t), …, IK(t)) is calculated in exactly the same way as for a disease with serial classes (refer to §3 and Figure 1). The joint probability function of the driving events (L1(t), …, LK(t)) can also be easily calculated noting that when L(t) individuals leave the susceptible class, then the driving events (L1(t), …, LK(t)) follow a multinomial distribution with L(t) number of trials and probability of success (α1, …, αK).

Finally, the assumption of exponential stay time in a class might be violated for some classes; as an alternative the modeler may choose Gamma or empirical distribution (Wearing et al., 2005). A Gamma distribution with scale parameter μ and an integer shape parameter K, can be modeled by adding K serial classes each with stay time exponentially distributed with rate μ. An empirical distribution may be modeled by adding a number of parallel classes (branches) to the model as in Figure 2.

5 Reducing the State Space through State Aggregation

For a population of size N, the transition probability matrix of the Markov chain {(XC1(t), …, XCM(t)): t = 1, 2, …} can be of size (N + 1)M−1 × (N + 1)M−1, which for large populations causes computational problems. Therefore, as the final step in constructing the discrete-time Markov models for a large population, we propose an approach for approximating the potentially enormous Markov chain {(XC1(t), …, XCM(t)): t = 1, 2, …} with the Markov chain {(ΘC1(t), …, ΘCM(t)): t = 1, 2, …}, where ΘCi(t), i [set membership] {1, 2, …, M} is the proportion of the population in class Ci at time t, and can only take a limited number of values from the set { θ1Ci,θ2Ci,,θdCiCi}. To properly determine the set { θ1Ci,θ2Ci,,θdCiCi}, let the control points {0=b0Ci,b1Ci,b2Ci,,bdCiCi=1},0b1Ci<b2Ci<<bdCiCi=1, divide the interval [0, 1] into dCi regions; the possible values of ΘCi(t), i.e., { θ1Ci,θ2Ci,,θdCiCi}, can now be determined by, θjCi=(bj1Ci+bjCi)/2 for j [set membership] {1, …, dCi}.

The transition probabilities for the Markov chain {(ΘC1(t), …, ΘCM(t)): t = 1, 2, …} can now be calculated by:

Pr{(ΘC1(t+Δt),,ΘCM(t+Δt))=(θk1C1,,θkMCM)(ΘC1(t),,ΘCM(t))}=(x1,,xM)ΩΘ(t)P(XC1(t+Δt),,XCM(t+Δt))((x1,,xM)NΘC1(t),,NΘCM(t))
(34)

where the probability mass function of (XC1(t + Δt), …, XCM(t + Δt)) conditional on the current state (XC1(t), …, XCM(t)) = ([left floor]NΘC1(t)[right floor], …, [left floor]NΘCM(t)[right floor]) is given by (10), and

ΩΘ(t)={(x1,,xM)NMNbk11C1x1Nbk1C1,,NbkM1CMxMNbkMCM}.

While the probability for some state transitions calculated by (34) may be positive, in some situations, the probability of moving from state (ΘC1(t), …, ΘCM(t)) to (ΘC1(tt), …, ΘCM(tt)) should be set to zero in the transition probability matrix of the Markov chain {(ΘC1(t), …, ΘCM(t)): t = 1, 2, …}. This situation may occur when transitions from state (ΘC1(t), …, ΘCM(t)) to (ΘC1(t + Δt), …, ΘCM(t + Δt)) are not possible (for instance, in a SIR model with closed population, the number of susceptibles at time t + Δt cannot exceed the number of susceptible at time t) but the probability support ΩX(t) in (10) and the set ΩΘ(t) in (34) intersect, i.e., ΩX(t) ∩ ΩΘ(t)[empty], which may result in a positive transition probability in (34). This can arise as a result of approximating the original Markov chain {(XC1(t), …, XCM(t)): t = 1, 2, …} with the reduced Markov chain {(ΘC1(t), …, ΘCM(t)): t = 1, 2, …} and, if not corrected, may lead to an invalid transition matrix. Accordingly, if the probability transition from one state to another should be set to zero in the transition matrix, then this should be ensured after calculating the transition matrix of the Markov chain {(ΘC1(t), …, ΘCM(t)): t = 1, 2, …}; in the attached e-companion, we describe the further simple steps to construct a valid transition probability matrix for these cases using the transition probability function (34).

6 An Illustrative Example: SIR Model

In this section, we use the SIR model developed in §4.2 to capture an influenza outbreak in an English boarding school reported in (Anonymous, 1978) and recently used by Merl et al. (2009) and Ludkovski and Niemi (2010). The population consisted of N = 763 students and the infection was believed to be introduced by one student returning from Asia. The situation satisfies many requirements of a simple SIR model, particularly since no specific intervention was employed during the outbreak.

Given the population size N = 763, the SIR Markov model in §4.2 will have N(N + 1)/2 = 291,466 states. Here we show how the state aggregation method described in §5 can be used to approximate this Markov chain. We also investigate the effect of such approximation on the power of the model in fitting the observed data. The influenza spread model described in §4.2 has three parameters: the contact rate λ, infectivity r, and the mean recovery time 1/μ. Throughout this section, we use maximum likelihood estimation (as discussed in §3.2) to find these parameters.

To construct the approximate Markov chain {(ΘS(t), ΘI(t)), t = 1, 2, …}, we must specify the control points { 0=b0Ci,b1Ci,b2Ci,,bdCiCi=1}, Ci [set membership] {S, I}, dividing the interval [0, 1] for the proportion of susceptibles and infectives. The simplest way to discretize the proportion of susceptibles and infectives is to distribute the control points uniformly over the interval [0, 1]. Figure 3(a) displays the number of students infected and the expected prevalence of influenza provided by the model over 13 days. In this figure, the proportions of susceptibles and infectives have been discretized by control points {0, 0.05, 0.10, 0.15, …, 0.95, 1}. The Markov chain generated has only 210 states (a 99.93% reduction in the state space compared with the model without approximation).

Figure 3
Influenza outbreak in the English boarding school

Although the expected number of infectives provided by the Markov model in Figure 3(a) fits the data reasonably well, the variance of this model output is substantial. A tighter fit can be obtained by selecting a finer grid for model approximation. Figure 3(b) displays the number of students infected and the expected prevalence of influenza provided by the model when the proportions of susceptibles and infectives have been discretized by control points {0, 0.02, 0.04, 0.06, …, 0.98, 1}. The generated Markov chain has 1275 states (a 99.56% reduction in the state space). The model’s parameters are estimated to be λ = 15 contacts per day, r = 11.5%, and 1/μ = 1.8 days.

Note that in both Figure 3(a) and Figure 3(b), the model is not able to perfectly capture the observed prevalence at time t = 1. This occurs because the number of infectives at time t = 1, XI(1) = 0, is too low to be captured by the approximate Markov chain {(ΘS(t), ΘI(t)): t = 0, 1, …}. The model starts fitting the data after time t = 1.

We note that the grid size in this illustrative example is arbitrary and was chosen solely to demonstrate the effect on model outputs. Choice of grid size will ideally reflect the intended use of the model. For example, to predict the magnitude of an outbreak, grid size may be selected to minimize the variance of the prediction, whereas, to determine health policies for controlling an outbreak, grid size may be selected to minimize the expected decision errors caused by approximation.

Table 1 summarizes the effect of grid choice on parameter estimation (see §3.2) and the predictive power of the model (see §3.3) when fitted to an initial set of observations. In this example, parameter estimates are not very sensitive to grid choice; however, choosing a finer grid significantly lowers the standard deviations of predicted values (i.e. XI(5) and XI(6) when the model is fit to the initial 4 and 5 observations, respectively). This can also be observed in Figure 3.

Table 1
Parameter estimation using an initial set of observations with different grids

Finally, to demonstrate the computational burden of the proposed SIR model, we compare the times needed to built the model for different population sizes in Figure 4. The SIR models here use the same parameter values as the model of Figure 3(b) with a relatively fine grid {0, 0.025, 0.05, 0.075, …, 0.975, 1} for both the number of susceptibles and infectives.

Figure 4
Computation time for the SIR model with 1275 states and grid {0, 0.025, 0.05, …, 0.975, 1}. The model is implemented in Microsoft Excel 2007, coded in Visual Basic for Applications (VBA), and run on a computer with a 2.80 GHz CPU.

7 Conclusion

Over the past several years, substantial effort has been devoted to developing methods that facilitate real-time decision making over the course of an epidemic as new information become available. Dynamic programming (Winston, 2003; Bertsekas, 2005), a technique that aims to aid decision making in changing environments, has been successfully applied to problems with similar characteristics (Schaefer et al., 2005; Van and Dana, 2003; Winston, 2003). To date, dynamic programming techniques have not been exploited in efforts to aid policy choices for the control of emerging or persistent infectious diseases. One reason for this may be that currently available infectious disease models do not generally satisfy the requirements of these techniques, as discussed in §2.

In this paper, we proposed a class of models for the spread of infectious diseases in relatively large populations. This class of models generalizes the discrete-time Markov chain models in the work of Reed-Frost (Abbey, 1952) and Greenwood (Greenwood, 1931) that have been used to describe the spread of pathogens with infectious periods that are relatively short in comparison with their latent periods (Daley and Gani, 1999). The proposed models possess two advantages over existing models of infectious diseases: (1) these models can be effectively used by dynamic optimization methods to select optimal dynamic health policies, and (2) they are able to approximate the spread of the disease in relatively large populations with a limited state space size, and reasonable degree of accuracy and computation time. In a related paper, we have demonstrated how the proposed class of models can be employed by an MDP to generate stationary optimal health policies for controlling an simplified Influenza epidemic (Yaesoubi and Cohen, 2010).

Although the primarily motivation in developing the proposed class of Markov models was to use them in generating dynamic health policies, the framework developed here extends the application of stochastic mathematical models to populations with relatively large size. The current stochastic mathematical models of infectious diseases have been limited to small and moderate size populations (Keeling and Ross, 2008), and hence, to date, computer simulation has been the primarily modeling approach for capturing the stochastic dynamics of infectious disease spreads in large populations.

The primarily disadvantage of using a computer simulation model is the need for a large number of replications to determine the stochastic behavior of the disease spread. The proposed class of Markov models here provides a mathematical framework to capture the stochastic dynamics of infectious disease spreads in relatively large population while allowing the modeler to incorporate the desired level of complexity in the representation of the within-host natural history of disease and between-host transmission dynamics. If, however, disease spread has negligible stochasticity, deterministic models would be preferred to the stochastic framework proposed in this paper.

The main disadvantage of our proposed framework is the need to use grids to reduce the state size of the Markov model for modeling disease spread in large populations. Methods to determine the optimize grid size to minimize approximation errors and developing alternative approximation techniques are promising topics for future research.

Research Highlights

  • We consider dynamic programming for the optimal control of infectious spreads.
  • The major limitations of existing infectious disease models are discussed.
  • We propose a class of models which can be employed by DP or approximate DP.
  • We demonstrate the ability of these models to fit data from an emerging epidemic.

Supplementary Material

01

Acknowledgments

The authors would like to thank Marc Lipsitch for his comments and suggestions. The work is supported by NIH grants DP2OD006663 and U54GM088558. The content is solely the responsibility of the authors and does not necessarily represent the official views of the Office of the Director of the National Institute of Health and the National Institute of General Medical Sciences of the National Institutes of Health.

Footnotes

1The support of a distribution is the smallest closed set whose complement has probability zero.

2By an alternative reasoning, since all the infectives at time t will be removed at time t + 1, infectives at time t + Δt are only those who were infected during the period [t, t + Δt]; that is, NXS(t + Δt) = I(t), which is exactly the same as constraint (21).

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Contributor Information

Reza Yaesoubi, Harvard School of Public Health - Department of Epidemiology, 677 Huntington Ave., Boston, MA 02115, U.S.A.

Ted Cohen, Brigham and Women’s Hospital - Division of Global Health Equity, Harvard School of Public Health - Department of Epidemiology, 641 Huntington Ave., Boston, MA 02115, U.S.A.

References

  • Abbey Helen. An examination of the reed-frost theory of epidemics. Human Biology. 1952;24(3):201–233. [PubMed]
  • Allen Linda JS. Some discrete-time SI, SIR, and SIS epidemic models. Mathematical biosciences. 1994;124(1):83–105. [PubMed]
  • Allen Linda JS, Burgin Amy M. Comparison of deterministic and stochastic sis and sir models in discrete time. Mathematical Biosciences. 2000;163:1–33. [PubMed]
  • Anderson Roy M, May Robert M. Infectious Diseases of Humans: Dynamics and Control. Oxford University Press; 1992.
  • Influenza in a boarding school. British Medical Journal. 1978:587. [PubMed]
  • Bertsekas Dimitri P. Dynamic Programming and Optimal Control. 3. Vol. 1. Athena Scientific; 2005.
  • Castillo-Chavez Carlos, Yakubu Abdul-Aziz. Discrete-time SIS models with complex dynamics. Nonlinear Analysis-Theory Methods and Applications. 2001;47(7):4753–4762.
  • Daley Daryl J, Gani Joseph M. Epidemic Modelling: An Introduction. Cambridge University Press; Cambridge; New York: 1999.
  • Dimitrov N, Goll S, Meyers LA, Pourbohloul B, Hupert N. Optimizing tactics for use of the US antiviral strategic national stockpile for pandemic (H1N1) Influenza, 2009. PLoS Curr Influenza 2009 [PMC free article] [PubMed]
  • Ge L, Kristensen AR, Mourits MC, Huirne RB. A new decision support framework for managing foot-and-mouth disease epidemics. Annals of Operations Research. 2010:1–14.
  • Goldstein E, Apolloni A, Lewis B, Miller JC, Macauley M, Eubank S, Lipsitch M, Wallinga J. Distribution of vaccine/antivirals and the ‘least spread line’ in a stratified population. Journal of the Royal Society Interface. 2010;7(46):755–764. [PMC free article] [PubMed]
  • Greenwood M. On the statistical measure of infectiousness. The Journal of Hygiene. 1931;31(3):336–351. [PMC free article] [PubMed]
  • Halloran ME, Ferguson NM, Eubank S, Longini IM, Cummings DAT, Lewis B, Xu S, Fraser C, Vullikanti A, Germann TC, et al. Modeling targeted layered containment of an influenza pandemic in the United States. Proceedings of the National Academy of Sciences. 2008;105(12):4639–4644. [PubMed]
  • Hethcote Herbert W. Mathematical Understanding of Infectious Disease Dynamics, chap. The basic epidemiology models: models, expressions for R0, parameter estimation, and applications. World Scientific Publishing Company; 2008. pp. 1–61.
  • Hethcote HW. The mathematics of infectious diseases. SIAM review. 2000;42(4):599–653.
  • Jacquez John A, O’Neill Philip. Reproduction numbers and thresholds in stochastic epidemic models i. homogeneous populations. Mathematical Biosciences. 1991;107(2):161–186. [PubMed]
  • Jacquez John A, Simon CP. The stochastic SI model with recruitment and deaths I. Comparison with the closed SIS model. Mathematical Biosciences. 1993;117(1–2):77–125. [PubMed]
  • Keeling MJ, Ross JV. On methods for studying stochastic disease dynamics. J R Soc Interface. 2008;5:171–181. [PMC free article] [PubMed]
  • Larson Richard C. Simple models of influenza progression within a heterogeneous population. Operations Research. 2007;55(3):399412.
  • Lefevre C. Optimal control of a birth and death epidemic process. Operations Research. 1981;29(5):971–982. [PubMed]
  • Ludkovski Michael, Niemi Jarad. Optimal dynamic policies for influenza management. Statistical Communications in Infectious Diseases. 2010;2(1)
  • Merl D, Johnson LR, Gramacy RB, Mangel M. A statistical framework for the adaptive management of epidemiological interventions. PloS One. (4) 2009;(6):e5087. [PMC free article] [PubMed]
  • Nåsell Ingemar. Stochastic models of some endemic infections. Mathematical Biosciences. 2002;179(1):1–19. [PubMed]
  • Powell WB. Approximate Dynamic Programming: Solving the curses of dimensionality. Wiley-Interscience; 2007.
  • Schaefer A, Bailey M, Shechter S, Roberts M. Modeling medical treatment using markov decision processes. Vol. 23. Springer; New York: 2005. Operations Research and Health Care -A Handbook of Methods and Applications; pp. 593–612.
  • Van Cuong Le, Dana Rose-Anne. Dynamic Programming in Economics. Kluwer Academic Publishers; Boston: 2003.
  • Wallinga J, van Boven M, Lipsitch M. Optimizing infectious disease interventions during an emerging epidemic. Proceedings of the National Academy of Sciences. 2010;107(2):923–928. [PubMed]
  • Wearing HJ, Rohani P, Keeling MJ. Appropriate models for the management of infectious diseases. PLoS Medicine. 2005;2(7):621. [PMC free article] [PubMed]
  • Winston Wayne L. Operations Research: Applications and Algorithms. 4. Duxbury Press; 2003.
  • Yaesoubi Reza, Cohen Ted. Dynamic Health Policies for Controlling the Spread of Emerging Infections: Influenza as an Example. 2010. Submitted. [PMC free article] [PubMed]