Modeling Susceptible Infective Recovered Dynamics and Plague Persistence in California Rodent

Abstract

Plague persists as an enzootic in several very different rodent–flea communities around the world. In California, a diversity of rodent–flea communities maintains the disease, and a single-host reservoir seems unlikely. Logistic regression of plague presence on climate and topographic variables predicts plague in many localities where it is absent. Thus, a dynamic community-based analysis was needed. Deterministic Susceptible Infective Recovered (SIR) models were adapted for plague and analyzed with an eye for insights concerning disease persistence. An R simulation program, Plaguesirs, was developed incorporating multihost and multivector SIR dynamics, demographic and environmental stochasticity, density dependence, and seasonal variation in birth and death. Flea–rodent utilization matrices allowed us to get transmission rates as well as flea carrying capacities. Rodent densities allowed us to estimate host carrying capacities, while maximum birth rates were mainly approximated through an examination of litter phenology and demography. We ran a set of simulations to assess the role of community structure in maintaining plague in a simulated version of Chuchupate campground in Ventura County. Although the actual campground comprises 10 rodent and 19 flea species, we focused on a subset suspected to act as a reservoir community. This included the vole Microtus californicus, the deer mouse Peromyscus maniculatus, the Ceratophyllid fleas Aetheca wagneri and Malareus telchinum, and the Leptopsyllid flea Peromyscopsylla hesperomys. The dynamics of 21 subsets of this community were simulated for 20 years. Single-rodent communities showed much lower disease persistence than two-rodent communities. However, so long as Malareus was present, endemicity was enhanced; removal of the other two fleas slightly increased disease persistence. Two critical features improved disease persistence: (1) host breeding season heterogeneity and (2) host population augmentation (due to two similar host species instead of one). Voles are winter–spring breeders compared to the spring–summer deer mice. While host synchronicity may enhance epidemics, host asynchronicity favors endemics.

Introduction

P lague appeared first in recorded human history in 541 CE at Pelusium on the eastern edge of the Nile delta during the reign at Constantinople of the emperor Justinian (Procopius 551 CE). Its agent, the bacterium Yersinia pestis, shipped around the Mediterranean with rats, causing epidemic waves about every 15 years for 200 years. The 15-year cycles possibly reflect the diminishment of herd immunity in successive human generations. The historical consequences may well include the end of antiquity and the onset of the Middle Ages in Europe (Little 2007, Rosen 2007). Plague then disappeared for 600 years from the historical record, reappearing as the Black Death in the 1340s in Central Asia and Crimea, reaching Italy just in time for the Renaissance. This second Pandemic lasted over 400 years, devastating perhaps half of Europe.

Plague reappeared as the third Pandemic in eastern Asia in the last part of the 19th century, spreading to much of the world. It reached western North America in 1899 at San Francisco and perhaps Los Angeles and Seattle (Adjemian et al. 2007). From these three sites waves of epidemic advance brought plague across western North America to the western edge of the short grass prairie within 40 years. The present distribution of endemic plague is, however, much more limited than its early epidemic scope.

In California, plague seems to be endemic in scattered mid- to high-altitude localities such as Chuchupate campground in Ventura County (1890 m), Nevada County (1951 m), Alturas in Modoc County (1330 m), and San Bruno Mountain in San Mateo County (400 m) (Davis et al. 2002, Vector Borne Disease Section 2006, Adjemian et al. 2008). Statistical analyses (using logistic regression on altitude, climate, etc.) fail to pinpoint specific endemic sites.

Outbreaks and diffusion of plague can be spectacular, rousing a great deal of public and scientific attention. Less spectacular, but more mysterious, are the long periods of endemicity in many parts of the world today, from Kazakhstan to Madagascar to South America to California.

The mathematical threshold R ₀ > 1 for epidemics is well understood, and diffusion requires only some flea or mammal mobility (Bailey 1982, Anderson and May 1991). Mathematical conditions for disease persistence are much trickier to state because they often depend on cyclic, nonlinear dynamics and on demographic and environmental stochasticity (Anderson and May 1991). Moreover, in the case of plague in North America, more than a single host and vector appear to be involved (Gage and Kosoy 2005).

The persistence of sylvatic plague in natural habitats appears to depend in a complicated way on a community of rodents and fleas (Gage and Kosoy 2005). One common theory is that some rodents are reservoirs with a low-mortality persistent infection. Other rodent species do not maintain the bacterium because the disease burns itself out in high-mortality epidemics. It is these epidemics that presumably catch the notice of wildlife disease researchers, and it is these epidemics that raise the likelihood of transmission to human commensals and campfollowers such as roof rats (Rattus rattus). According to this verbal theory, plague outbreaks require a rodent community with at least one reservoir (or “maintenance”) species and one “amplifying” species.

In California, there are about 94 rodent species (Zeiner et al. 1988) and 27 suspected plague-carrying flea species (Hubbard 1968). In a typical 1-km² area there may be a dozen rodent species and a dozen rodent-based flea species, fewer than half of which are competent to transmit the disease. At plague-ridden Chuchupate campground, 19 flea species were taken from 10 rodent species and a rabbit species. Of these, a woodrat, ground squirrel, chipmunk, and four flea species seemed especially susceptible to plague (Davis et al. 2002).

Several researchers have suggested that if we better understood rodent and flea distribution and phenology, we could better predict plague outbreaks (Eskey and Haas 1940). This paper attempts to model plague dynamics as a function of rodent and flea community characteristics. Seasonality is incorporated into the model by allowing birth and death rates to vary seasonally for hosts and vectors.

Bailey (1982) developed a “general” theory of host–vector disease dynamics based on the (Kermack and McKendrick 1927) model, but its generality extends only to one host and one vector species and lacks host mortality, fecundity, and seasonality. To examine the dynamics of plague in realistic communities, we need a community SIR model.

This paper develops a more detailed, realistic SIR model and its implementation as R code for the simulation function Plaguesirs. We do this in a series of steps, starting from a simple one-host SIR model, adding more reality along the way. At each step, we look for insight into plague persistence and extinction. In culmination, the community host–vector SIR model tracks a vector of h (rodent) hosts and a vector of v (flea) vectors. This model allows us to examine community structure effects on plague persistence. We do this by examining a critical subset of the Chuchupate campground. Running thousands of simulations of this community and its subsets, we conclude that community structure can have substantial effects on plague persistence.

Materials and Methods

The SIR dynamics are written out as a system of differential equations building on the work of Kermack and McKendrick, Bailey, and Anderson & May (Anderson and May 1991) and others. The system is simulated in the open source statistical language R based on S developed at Bell Labs (Chambers and Hastie 1992). R was also used to do the statistical analyses and create graphic plots.

The present version of the simulation Plaguesirs is based on an earlier version described more fully in (Foley et al. 2007). We have improved the demographic stochasticity approach and made other improvements. The R program simsirstoch is a simple one-host simulation of SIR with demography and environmental stochasticity. Demographic stochasticity was implemented using binomial random variables for transitions and Poisson random variables for births.

We used a variety of approaches to obtain plausible parameter values we need: a vector-to-host transmission matrix β, a host-to-vector transmission matrix β_V, a vector of host recovery rates γ, and so forth. We develop a plausible way, discussed later, to decompose these transmission matrices into quantities that are easier to estimate.

Demographic parameters including carrying capacities and environmental stochasticity are estimated using the approaches of Foley (1994, 2000). To adapt a yearly population multiplication rate R to a daily birth rate r _m, we use the formula R = exp(r _m T _b), where T _b is the length of the breeding season in days. For example, if woodrats can maximally produce three successful litters of 1.5 females in a 180 day breeding season, then R = 4.5 (ignoring survivorship), and r _m is ln(4.5)/180 = 0.00836 during the season.

Although we intend this models to be as general as possible, we often refer to hosts as rodents and vectors as fleas, partly to anchor the abstract arguments and partly to deal with the confusions inevitable in studying vectors of vectors.

Results

Kermack–McKendrick single-host SIR model with demography

We can gain some insight into the dynamics of plague by examining a single-host SIR model first developed by Kermack and McKendrick (1927). This will also serve as an introduction to the community model.

A population of hosts of size N can be decomposed into three groups of individuals, susceptibles, S, infectious, I, and resistant, R. Thus an SIR model with some demography can be written out like this: \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} N = S + I + R . \tag{1} \end{align*} \end{document} \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align} \begin{split} \frac {dS} {dt} & = bN - \beta SY - d_S S \\ \frac {dI} {dt} & = \beta SY - \gamma I - d_I I \\ \frac {dR}{dt} & = \gamma I - d_R R \end{split}\tag{2} \end{align} \end{document}

The transmission rate β is tricky to estimate in practice. The recovery rate γ is the reciprocal of the mean time to recovery. In rodent hosts γ can vary from 0.1 to 0.001 per day. Birth rates b and death rates d also vary widely among rodents.

This system of equations is well studied, but can only be applied to host–vector systems if we let the transmission term β depend on flea density, disease state, and behavior. This is a plausible approach for future modeling.

Epidemics depend upon the growth of I individuals. This can only occur in a naive host population when \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} R_0 = N \frac {\beta} {\gamma + d_I} > 1. \tag {3} \end{align*} \end{document}

When there is some herd immunity, an epidemic is harder to start since the initial number of susceptibles S ₀ must satisfy the same equation, replacing N with S ₀. It is unfortunate that there is no standard R ₀ notation for the general case of partial herd immunity.

Although deterministic conditions for endemicity do not make for a good fit to data, this is a good place to start. To do so, we set the system to equilibrium by requiring I > 1 and make a couple of temporary assumptions: d = d _S = d _R; the overall birth rate b compensates for the overall death rate and that all the rates are constant. At equilibrium we get \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align} \begin{split} S^{*} & = \frac{\gamma + d_I} {\beta} \\ R^{*} & = \frac{\gamma} {d} I \\ I^{*} & = \frac{N - \frac{\gamma + d_I} {\beta}} {1 + \frac{\gamma} {d}} \end{split}\tag {4} \end{align} \end{document}

To keep an endemic alive I must be at least one. Demographic (or other forms of ) stochasticity may require us to maintain a higher value I _e. Thus we get this threshold for endemicity: \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} N_e = \frac {\gamma + d_I} {\beta} + \left(1 + \frac {\gamma} {d} \right) I_e \tag {5} \end{align*} \end{document}

Note that γ is likely to be only on the order of 10 times d, so that the endemic threshold for N may not be much higher than the epidemic threshold N. Yet documented endemic plague is relatively rare in contrast with regions of epidemic success. This suggests that endemic fadeout is likely due to stochastic variation of I(t) around the equilibrium value I*.

We include stochasticity into the analysis for a better fit to data. Simulations done using our simple R program simsirstoch() show that I(t) is approximately normal around I*, that is \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} I (t)\, \sim \, N (I^{*}, \sigma_{I}^{2}) \tag{6} \end{align*} \end{document}

where σ_I ² rises quickly with environmental stochasticity in the overall population growth. These random I(t) values are not independent. The extinction of the disease in the host population will occur when the random I(t) drops to 0. The most important quantities for predicting the time to extinction of the disease are then I*/σ_I and the effective time between independent I(t) values. High I* and low stochasticity forestall endemic fadeout.

Real rodent populations show annual cycles often with greatest birth rates in spring/summer, and lowest populations in winter. Such populations exhaust most of their susceptibles early in the year, and the deterministic equilibrium value I* itself tends to decline until the spring birth pulse. Epidemic fadeout may only be avoided by a burst of new susceptibles.

This single-host SIR model assigns β a big job. Yersinia transmission from I individuals to S individuals is mediated by flea vectors. β thus varies with flea densities, temperatures, and other environmental factors, increasing fluctuations in I(t). On the other hand, the internal dynamics of vector disease dynamics may slow down the epidemic and thus increase the chances of endemicity. Some subpopulations of fleas or their larvae may be an effective long-term refuge for Yersinia. There is also likely to be a latent period before an infected flea is infectious.

If we cannot neglect flea dynamics, we must resort to a more complex model.

A modified Bailey single host–single vector SIR model

Bailey worked out a few more implications of Kermack and McKendrick's early host–vector SIR model for malaria (Bailey 1982). Here we apply these models to plague.

Similar to the SIR compartments of the host N, a flea vector population of size V can be separated into three groups, though the immunes, Z, may be rare or lacking in nature. \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} V = X + Y + Z. \tag{7} \end{align*} \end{document}

In a short time interval, individuals move from one host category to another by

flea-to-host transmission at a rate βSY,

host recovery at a rate γI,

host birth at a rate bN.

Similar rates apply to fleas, so that we can write the dynamics out as a system of differential equations. \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align} \begin{split} \frac{dS}{dt} &= bN - \beta SY - d_S S \quad \frac{dX}{dt} = b_V V - \beta_V XI - d_X X \\ \frac{dI}{dt} &= \beta SY - \gamma I - d_I I \quad \ \frac{dY}{dt} = \beta_V XI - \gamma_V Y - d_Y Y \\ \frac{dR} {dt} & = \gamma I - d_R R \quad\quad\quad\quad \frac {dZ}{dt} = \gamma_V Y - d_Z Z \end{split}\tag{8} \end{align} \end{document}

Host death rates are likely to be the same for S and R individuals. Further, resistant individuals Z apparently are rare so infectious individuals get recycled as susceptibles X. This allows a few simplifications. \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align} \begin{split} \frac{dS}{dt} & = bN - \beta SY - d_S S \quad \frac{dX}{dt} = b_V V + \gamma_V Y - \beta_V XI - d_X X \\ \frac{dI}{dt} &= \beta SY - \gamma I - d_I I \quad \ \frac{dY}{dt} = \beta_V XI - \gamma_V Y - d_Y Y\end{split} \tag{9} \end{align} \end{document} \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} \frac{dR}{dt} = \gamma I - d_S R \end{align*} \end{document}

With a little simplification, neglecting demography, Bailey was able to get an epidemic threshold \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} R_0 = \frac{NV \beta \beta_V} {\gamma \gamma_V} \ge 1 \tag {10} \end{align*} \end{document}

Note the importance of both rodent and flea population sizes in producing an epidemic. Protracted recovery times for rodents or fleas can increase the chance of epidemics. Since both β and β_V depend on flea attack rates, this component of transmission has a squared effect on epidemics.

Conditions for endemicity are complicated by the consideration that the disease can be maintained (temporarily) by either host or vector. The rough approaches of the last section should still apply. Besides temporal variations in birth and death rates, the potentially complicating effects of the entire rodent–flea community will cause fluctuations in I(t) and Y(t) that the single host model neglects.

Community SIR model

A rodent host guild and a flea vector guild make up the subset of a community that concerns us. Consider a vector of rodent population sizes \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} N = (N_{1}, N_{2}, \ldots , N_{h}) \end{align*} \end{document} and a vector of flea vectors \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} V = ( V_{1}, V_{2}, \ldots , V_{v}) \end{align*} \end{document} . Within the text we will often write out vectors horizontally, but as is standard, they are column vectors.

We use vector and matrix notation to write out our community model. Thus, the species-specific birth rates form diagonal matrices. \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} bN = \left[{\begin{matrix}b_1 & 0 & & 0 \\ 0 & b_2 & & 0 \\ \\ 0 & 0 & & b_h\end{matrix}} \right] \left[{\begin{matrix}N_1 \\ N_2 \\ \\ N_h\end{matrix}} \right] = \left[{\begin{matrix}b_1 N_1 \\ b_2 N_2 \\ \\ b_h N_h\end{matrix}} \right] \tag{11} \end{align*} \end{document}

Transmission to new hosts can take place from several possible flea species. Let β₁₂ be the transmission rate from vector 2 to host 1. Then, we can model the entire host transmission process according to Equation 12. \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} S1 \beta Y &= \left[{\begin{matrix}S_1 \\ S_2 \\ \vdots \\ S_h\end{matrix}} \right] \left[{\begin{matrix}1 & 1 & \ldots & 1\end{matrix}} \right] \left[{\begin{matrix}\beta_{11} & \beta_{12} & & \beta_{1v} \\ \beta_{21} & \beta_{22} & & \beta_{11} \\ & & \ddots & \\ \beta_{h1} & \beta_{h2} & & \beta_{hv}\end{matrix}} \right] \left[{\begin{matrix}Y_1 \\ Y_2 \\ \vdots \\ Y_v\end{matrix}}\right] \\&= \left[{\begin{matrix}\beta_{11} S_1 Y_1 + \beta_{12} S_1 Y_2 + \ldots + \beta_{1v} S_1 Y_v \\ \beta_{21} S_2 Y_1 + \beta_{22} S_2 Y_2 + \ldots + \beta_{2v} S_2 Y_v \\ \\ \beta_{h1} S_h Y_1 + \beta_{h2} S_h Y_2 + \ldots + \beta_{hv} S_h Y_v\end{matrix}} \right] \tag{12} \end{align*} \end{document}

Similarly, d _S, d _I, and γ _I are h*h matrices, and b _V, γ _V, d _X, and d _Y are v*v matrices. β_V is a v*h matrix giving the transmission rates from host species j to vector species i.

After fleas contract plague, there often is a substantial latency period during which they are not infective (but see Eisen et al. 2008). This leads in most cases to the added complication of another vector state E_v . Our new community SIR model includes this complication. \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} \frac{dS}{dt} & = bN - S1 \beta Y - d_S S\quad \frac{dX}{dt} = b_V V + \gamma_V Y - X1 \beta_V I - d_X X \\ \frac {dI}{dt} & = S1 \beta Y - \gamma I - d_I I \quad \frac{dE_V}{dt} = X1 \beta_V I - \psi E_V - d_X E_V \\ \frac{dR}{dt}& = \gamma I - d_S R \quad\quad \quad\qquad \frac{dY} {dt} = \psi E_V - \gamma_V Y - d_Y Y \\ \tag {13} \end{align*} \end{document}

The rate at which latency shifts to infectivity is ψ. The mean time in latency is then 1/ψ. We have assumed that the death rate of latents is the same as the susceptibles.

Translation to a simulation program Plaguesirs()

The system of Equation 13 is 3h + 3v dimensional and nonlinear. Although some insights can be gleaned from it by qualitative analysis, for many purposes a computer simulation is needed. This forces a series of compromises, since the simulation uses time discretely, while the original equations are set in continuous time. The accuracy of the simulation's ability to mimic the differential equation 13 is improved as the time chunks get smaller. Running the time in day increments seems a reasonable compromise that is easy to interpret. More complicated approximation approaches are discussed at length in (Keeling and Rohani 2008).

Differential equations can handle a transmission model such as βSY, but in discrete time, dangers arise. As infected vector populations build up so that βY > 1, it becomes unfortunately possible to produce more newly infected than the original susceptibles. To avoid this problem for high populations, while retaining the standard SIR model dynamics for small populations, Plaguesirs uses a Nicholson-Bailey type approach to model the surviving susceptible hosts. The rate at which one host individual contacts infective vectors is called the force of infection. The force of infection on individuals of species i is \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} \hbox{forceH}_i = \sum\nolimits_j {\beta_{ij} Y_j}.\tag{14} \end{align*} \end{document}

So, if infection is a Poisson process, only a certain fraction of the susceptibles can be expected to avoid infection completely \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} S_i \leftarrow S_i \exp ( -\hbox{forceH}_i). \tag{15} \end{align*} \end{document}

Overall we need \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} S \leftarrow S \exp (- \beta Y) \tag{16} \end{align*} \end{document}

which is a bit more elegant than the ugly matrix equation for transmission in Equation 18.

Modeling host birth and death

Rodent (and other host) populations manifest features of high growth rate, annual population fluctuations, density dependence, and seasonality. In some cases they also show population cycles (Turchin 1993, 2003). The simplest way to capture the most indispensable features (expect seasonality and cycles) uses a stochastic Ricker equation. \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} N ( t + 1 ) = N (t) e^{r(t)} = N (t) e^{r_{m}(1 - N/K) + \varepsilon_t} \tag{17} \end{align*} \end{document}

where ɛ₁∼N(0, v _r) is the effect of environmental stochasticity at time t (Foley 1997). The growth parameters r _m, K, and v _r can be estimated with some uncertainty from time series, but Equation 17 does not handle seasonality well because death rates [which are implicitly incorporated in r(t)] are often different for infected and uninfected hosts. This is implied in Equation 13, where d _s need not equal d _I.

Therefore, we need to retain carrying capacity and environmental stochasticity while teasing apart b(t) and d(t) from r(t). If death rate is constant over the SIR classes, this leads to the model \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} N(t + 1) = N (t) e^{r(t)} = N (t) e^{b(t) - d(t)} = N (t) e^{b_m(1- N(t)/K) + \varepsilon_t - d(t)}. \tag{18} \end{align*} \end{document}

If environmental stochasticity is low (v _r∼0) and seasonality can be neglected, then the populations size N reaches the equilibrium \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} N^* = K(1 - db_m). \tag{19} \end{align*} \end{document}

If d _i is much smaller than b _m, K is the carrying capacity according to our usual intuitions, but if mortality caused by the disease is high enough, the species can go locally extinct.

Death rates and birth rates can be approximated by a linear combination of cosine and/or sine functions (Chatfield 1999) (Keeling and Gilligan 2000), for example \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} d(t) = A \cos(t - t_{peak}) + {\bar d} \tag{20} \end{align*} \end{document}

where d(t) is the death rate at a time t counted in radians, that is t = 2π*real time, A is the amplitude of the death rate, and t _peak is the time of peak mortality. Birth rates in rodents typically are seasonal in several pulses during the favorable season. We implement this process by defining a birth window for each rodent species. Outside the window, b _m is 0.

Flea demography and seasonality is handled similarly in Plaguesirs. An additional flea problem is the population dependency of vectors on their hosts.

Central role of the host–vector count table

The natural way to organize community information about host–vector relationships is a table showing the observed counts of each vector species on each host species (Linsdale and Davis 1956) (Davis et al. 2002). With random host sampling, it is possible to estimate relative host abundances, relative vector abundances and host–vector relationships.

A host–vector count table provides useful data for constructing the transmission matrices β and β_V, and for parameterizing the densitterizing the density dependence model for the vectors. We can also adjust each row of the table to reflect the relative abundance of hosts, if we have such supplementary data.

The carrying capacity of a vector species j depends on the abundance of each host species i, together with some measure k_ij of its ability to use that host species as a resource. In the absence of any competitive fleas the carrying capacity K_j of vector j is given by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} K_{v_j} = \sum\nolimits_i k_{ij} N_i = \sum\nolimits_i u_{ij} q_i w_j N_i \tag{21} \end{align*} \end{document}

Transmission rates and density dependence are all affected by the utilization rate or preference of a particular vector species j for a particular host species i. We call this preference u_ij . It is zero if vector j does not utilize host i. An example u matrix for rodent rows Microtus californicus and Peromyscus maniculatus is given below. The column fleas are Aetheca wagneri, Malareus telchinum, and Peromyscopsylla hesperomys, respectively. \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} u = \left[{\begin{matrix}0.0065 & 0.838 & 0.0 \\ 0.9935 & 0.162 & 1.0\end{matrix}} \right]\tag{22} \end{align*} \end{document}

The carrying capacity offered by host i to vector j, k_ij , depends on: (1) the vector's host utilization rate u_j , (2) the size of the host individuals or some other index of the quality of resources a host provides q_i , and (3) the diminutiveness of the vector w_j . w_j is close to unity for most flea species, but may rise greatly above 1 for small fleas.

To transmit a disease from an infective vector j to a susceptible host i, the vector needs to (1) find a host (proportional to aSY), (2) be competent to transmit the disease (with vector competency c_vj ), and (3) want to use the host (with utilization rate u_ij ). The vector-to-host matrix β has entries that look like this. \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} \beta_{ij} = ac_{vj} u_{ij} \tag{23} \end{align*} \end{document}

To transmit a disease from a host j to a vector i, similar reasoning applies, and the host-to-vector matrix β_v has entries of this form. \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} \beta_{v_{ij}} = ac_j u_{ji}. \tag{24} \end{align*} \end{document}

Community structure affects plague persistence

We constructed a subset of the rodent–flea of Chuchupate campground (Ventura County, California) to dissect out the community features that encourage plague persistence. A reconstruction of the entire rodent–flea is parameterized and discussed in an earlier paper (Foley et al. 2007). Although the actual campground comprises 10 rodent species and 19 flea species, we focused on a subset that we suspected to act as a reservoir community. This included the vole Microtus californicus, the deer mouse Peromyscus maniculatus, the Ceratophyllid fleas Ae. wagneri and M. telchinum, and the Leptopsyllid flea P. hesperomys with the u matrix given in Equation 22. These species were chosen because they have been reported as likely reservoirs and are not reported to experience high mortality (and thus likely epidemic fade-out), but have been reported to be culture positive in plague-endemic areas in the absence of morbidity, which is a very unusual occurrence (Stoenner et al. 1959, Marchette et al. 1962, Kartman et al. 1966). In a recent study at Chuchupate, these rodents also were reported as reservoirs, with three additional species (the dusky-footed woodrat, the Merriam's chipmunk, and the California ground squirrel), all considered to be “susceptible,” that is, not reservoirs (Davis et al. 2002).

The dynamics of 21 distinct subsets of this community were simulated for 20 years with Plaguesirs. The full five-species community became endemic. We then slightly modified parameter values to a mean disease persistence of 10 years for the five-species community. Using these same parameter values, smaller communities varied in their capacity to maintain plague as can be seen in Table 1.

Table 1.

Simulated Plague Persistence Times T _e (Measured in Days) in Subsets of a Rodent–Flea Community at Chuchupate Campground, Ventura County, California

Host breeding seasons different: voles, winter-spring; deer mice, spring-summer
		Vectors
		1, 2, 3	1, 2	1, 3	2, 3	1	2	3
Hosts	1, 2	3652 (82)	3850 (84)	908 (63)	3603 (83)	997 (69)	4813 (86)	732 (37)
	1	830 (68)	786 (76)	36 (4)	833 (69)	39 (4)	918 (73)	39 (4)
	2	848 (59)	821 (56)	864 (53)	890 (62)	836 (58)	87 (10)	807 (50)

Host breeding season identical: voles and deer mice, both spring-summer
		Vectors
		1, 2, 3	1, 2	1, 3	2, 3	1	2	3
Hosts	1, 2	1886 (146)	1593 (132)	1052 (71)	1967 (159)	848 (53)	1512 (121)	705 (41)
	1	962 (64)	803 (62)	50 (4)	829 (62)	49 (5)	857 (85)	44 (4)
	2	849 (48)	680 (38)	918 (50)	858 (61)	792 (58)	121 (13)	788 (44)

The values in the cells are mean values from either 100 or 800 runs. In parentheses are standard errors for the means. Host 1 is the vole Microtus californicus; host 2 is the deer mouse Peromyscus maniculatus. Vector 1 is the flea Aetheca wagneri; vector 2 is Malareus telchinum; vector 3 is Peromyscaopsylla hesperomys.

Communities consisting of a single-rodent species showed much lower endemicity than the fuller community. Although rodent diversity enhanced persistence, this was not so for fleas. So long as Malareus was present, endemicity was enhanced. In fact removal of the other two fleas slightly increased disease persistence. The other two fleas greatly favor deer mice, whereas Malareus is more evenhanded.

Two distinct community features contribute to enhanced disease persistence in the larger host community: breeding season extension and host population augmentation. Voles are winter–spring breeders compared to the spring–summer deer mice. Voles save the disease from extinction in the winter and deer mice in the summer, as can be seen in Figure 1. With our parameters, this adds 2.4 months per year of disease persistence security.

FIG. 1.

Typical 2-year output from Plaguesirs with parameter values used to produce Table 1. Microtus and Aetheca population sizes are shown with solid lines. Broader dashes belong to Peromyscus and Malareus. Peromyscopsylla is shown with shorter dashes.

The two-host system has the added advantage (from the bacterial point of view) of doubling the potential number of individual hosts. By itself this effect should diminish the effect of demographic stochasticity and increase persistence times. We examined the host augmentation effect by artificially shifting the voles to the same breeding season as the deer mice, March to September. Disease persistence times remained the same for the single-host communities as can be seen in Table 1. The two-host communities with Malareus dropped their persistence times from about 10 years to about 5. Thus, in our specific community, the host augmentation effect and the breeding season extension effect contributed about equally to the increased endemicity of plague in the two-host community.

Discussion

Disease persistence in rodent–flea communities

The SIR models we have examined provide several insights into plague persistence. The one-host deterministic SIR model predicts endemic thresholds that are only a little higher than epidemic thresholds. Our simulation simsirstoch shows that the addition of even a little environmental stochasticity leads to greatly increased disease extinction rates if the deterministic equilibrium I* is close to 1. If a host population has a seasonal low point in I*, then disease extinction rates will be roughly proportional to the length of the bad season. The deterministic host–vector SIR model has an epidemic threshold analogous to the simpler SIR model, but both host and vector numbers matter. Endemic thresholds should also be analogous.

Our culminating model is a community host–vector SIR model that we translated into a simulation models Plaguesirs. To investigate community effects on plague persistence, we used Plaguesirs to simulate a two-rodent/three-flea subset of the Chuchupate rodent–flea community.

Our simulated subset of the Chuchupate campground rodent–flea community is a simplification that captures only a small part of the heterogeneity that promotes plague persistence in real California communities. Even in this simplified system, we found that a multiple-rodent host community can increase disease persistence due to a simple augmentation of available host individuals and by extending the effective host breeding season. The host augmentation effect diminishes endemic vulnerability to demographic stochasticity (endemic fadeout), since the deterministic values of I* are raised higher above the threshold of 1. The window in which the pathogen becomes vulnerable to extinction is shortened as the host breeding season is extended. In California, plague appears to be endemic in regions of high topographic relief. Such landscapes have high habitat diversity, even on a local scale, increasing the opportunities for rodents to diverge in breeding phenology. On a larger scale, montane landscapes create even more heterogeneity in host phenology and herd immunity. Thus, rodent metapopulations and metacommunities may act as reservoirs as has been argued, for very different landscapes by (Keeling and Gilligan 2000, Davis et al. 2007).

High flea diversity did not translate into high plague persistence in our model community. Ae. wagneri and P. hesperomys have a strong affinity for deer mice, while M. telchinum leans toward voles, but encounters many deer mice. Plague persists most strongly in the absence of the two deer mouse specialists, and we are not sure why. Higher vector loads may lead to a hotter early epidemic with greater likelihood of epidemic fadeout.

Parameter estimation problems

The models we have presented each give some insight into the persistence of plague in natural populations. We have attempted to build increasingly accurate models, throwing a great deal of natural history into the mix. Unfortunately, the more precise and detailed the model, the more parameters need to be estimated. Our earlier paper deals with the estimation problems in greater detail (Foley et al. 2007), but our results suggest a few strategic approaches.

Field epidemiologists and land managers may want to look for evidence of rodent host heterogeneities that increase the persistence of plague, such as breeding season differences and high seasonal mortality differences. There appears to be local variation in mortality attributable to disease, such that plague endemic rodents often show lower mortality attributable to disease (Quan and Kartman 1962, Hudson et al. 1964). Therefore, it is dangerous to generalize our findings across various locations without making local estimates of parameters. The parameter estimation problems loom especially large for the thought experiment we would most like to realize: using Plaguesirs to predict plague presence and plague absence across the California landscape.

Three problems for future modeling

Plaguesirs includes intraspecific competition in the density-dependent birth rates of both rodents and fleas; however, no interspecific competition is modeled. Various fudges could be used, but the problem is a deep one. How vectors divide up the resources of an individual host is hard to know, but how they partition all the resources provided by a host community is a very multidimensional problem, one that might drive even Santa Rosalia to drink (Hutchinson 1959).

Our flea populations are distributed among host species in an empirically measurable way. Unfortunately, we do not consider the fidelity of individuals to their individual host. This could be especially problematic when mammal hosts are not dying, because flea host-fidelity is strongest then (Gage and Kosoy 2005). To fix this with our model may require a more compartmented rodent-oriented SIR model with fleas attached to categories of hosts, rather than having their own differential equations. It is our intuition that such models may be (counterintuitively!) actually easier to model analytically at least approximately.

Our original interest was to see if community SIR models could improve the local persistence of plague, and we have presented evidence that this is so. However, if we want to study more realistic models in a landscape, then we need to convert our models to metapopulation models. At first glance this seems like adding one more layer of confusion and aggravation, but we have developed similar models for Hendra virus in Australian Flying Foxes (Plowright et al., unpublished data), and the results should be worth the effort. In a best-case scenario, a metacommunity SIR model be easier to analyze and even connect better to the real world.

Footnotes

Acknowledgments

We would like especially to thank members of Janet Foley's lab at UCDavis for energetic discussion. Encouragement and funds to attend the 2008 Fort Collins Plague Ecology Symposium are also greatly appreciated.

Disclosure Statement

Neither author has any conflict of interest that could influence this research or the writing of this article.

References

Adjemian

, Adjemian

, Foley

, Chomel

et al. Evidence of multiple zoonotic agents in a wild rodent community in the eastern Sierra Nevada. J Wildl Dis, 2008; 44:737–742.

Adjemian

, Foley

, Gage

, Foley

. Initiation and spread of traveling waves of plague, Yersinia pestis, in the western United States. Am J Trop Med Hyg, 2007; 76:365–375.

Anderson

, May

. Infectious Diseases of Humans. Oxford: Oxford Univ. Press, 1991.

Bailey

. The Biomathematics of Malaria. London: C. Griffen, 1982.

Chambers

, Hastie

. Statistical Models in S. Wadsworth & Brooks. California: Pacific Grove, 1992.

Chatfield

. The Analysis of Time Series: An Introduction. London, Boca Raton, Florida: Chapman and Hall/CRC, 1999.

Davis

, Klassovskiy

, Ageyev

, Suleimenov

et al. Plague metapopulation dynamics in a natural reservoir: the burrow system as the unit of study. Epidemiol Infect, 2007; 135:740–748.

Davis

, Smith

, Madon

, Sitko-Cleugh

. Flea, rodent, and plague ecology at Chuchupate Campground, Ventura County, California. J Vector Ecol, 2002; 27:107–127.

Eisen

, Holmes

, Schotthoefer

, Vetter

et al. Demonstration of early-phase transmission of Yersinia pestis by the mouse flea, Aetheca wagneri (Siphonaptera: Ceratophylidae), and implications for the role of deer mice as enzootic reservoirs. J Med Entomol, 2008; 45:1160–1164.

10.

Eskey

, Haas

. Plague in the western part of the United States. Public Health Bull, 1940; 254:83.

11.

Foley

. Predicting extinction times from environmental stochasticity and carrying capacity. Conserv Biol, 1994; 8:124–137.

12.

Foley

. Extinction Models for Local Populations. 215–246 Hanski

, Gilpin

. Metapopulation Biology. San Diego: Acad. Press, 1997.

13.

Foley

. Problems in extinction model selection and parameter estimation. Environ Manage, 2000; 26:S55–S73.

14.

Foley

, Zipser

, Chomel

, Girvetz

, Foley

. Modeling plague persistence in host-vector communities in California. J Wildl Dis, 2007; 43:408–424.

15.

Gage

, Kosoy

. Natural history of plague: perspectives from more than a century of research. Annu Rev Entomol, 2005; 50:505–528.

16.

Hubbard

. Fleas of Western North America, Their Relation to the Public Health. New York: Hafner Pub. Co., 1968.

17.

Hudson

, Quan

, Goldenberg

. Serum antibody responses in a population of Microtus californicus and associated rodent species during and after Pasteurella pestis epizootics in the San Francisco bay area. Zoonoses Res, 1964; 3:15–29.

18.

Hutchinson

. Homage to Santa Rosalia or why are there so many kinds of animals? Am Nat, 1959; 93:145–159.

19.

Kartman

, Goldenberg

, Hubbert

. Recent observations on the epidemiology of plague in the United States. Am J Public Health, 1966; 56:1554–1569.

20.

Keeling

, Gilligan

. Bubonic plague: a metapopulation model of a zoonosis. Proc R Soc Lond B, 2000; 267:2219–2230.

21.

Keeling

, Rohan

. Modeling Infectious Diseases. New Jersey: Princeton Univ. Press, 2008.

22.

Kermack

, McKendrick

. Contributions to the mathematical theory of epidemics. Proc R Soc, 1927; 115A:700–721.

23.

Linsdale

, Davis

. Taxonomic appraisal and occurrence of fleas at the Hastings Reservation in Central California. Univ Calif Publ Zool, 1956; 54:293–370.

24.

Little

. Plague and the End of Antiquity. Cambridge: Cambridge University Press, 2007.

25.

Marchette

, Lundgren

, Nicholes

, Bushman

, Vest

. Studies on infectious diseases in wild animals in Utah. II. Susceptibility of wild mammals to experimental plague. Zoonoses Res, 1962; 1:225–250.

26.

Procopius. The Persian War, 551 CE.

27.

Quan

, Kartman

. Ecological studies of wild rodent plague in the San Francisco Bay area of California. VIII. Susceptibility of wild rodents to experimental plague infection. Zoonoses Res, 1962; 1:121–144.

28.

Rosen

. Justinian's Flea: the First Great Plague and the End of the Romen Empire. London: Penguin, 2007.

29.

Stoenner

, Holdenried

, Lackman

, Orsborn

Jr.

The occurrence of Coxiella burnetil, Brucella, and other pathogens among fauna of the Great Salt Lake Desert in Utah. Am J Trop Med Hyg, 1959; 8:590–596.

30.

Turchin

. Chaos and stability in rodent population dynamics evidence from non-linear time-series analysis. Oikos, 1993; 68:167–172.

31.

Turchin

. Complex Population Dynamics: a Theoretical Empirical Synthesis. Princeton, N.J.: Princeton University Press, 2003.

32.

Vector Borne Disease Section. Vector-borne diseases in California. 2004 Annual Report, 58. Infectious Diseases Branch, California Department of Health Services: Sacramento CA, 2006.

33.

Zeiner

, Laudenslayer

, Mayer

, Program

CWHR

. California's Wildlife: Amphibians, Birds, Mammals, Reptiles. Sacramento, CA: California Dept. of Fish and Game, 1988.

Modeling Susceptible Infective Recovered Dynamics and Plague Persistence in California Rodent–Flea Communities

Abstract

Introduction

Materials and Methods

Results

Kermack–McKendrick single-host SIR model with demography

A modified Bailey single host–single vector SIR model

Community SIR model

Translation to a simulation program Plaguesirs()

Modeling host birth and death

Central role of the host–vector count table

Community structure affects plague persistence

Discussion

Disease persistence in rodent–flea communities

Parameter estimation problems

Three problems for future modeling

Footnotes

Acknowledgments

Disclosure Statement

References