Hamiltonian Analysis of Subcritical Stochastic Epidemic Dynamics.

We extend a technique of approximation of the long-term behavior of a supercritical stochastic epidemic model, using the WKB approximation and a Hamiltonian phase space, to the subcritical case. The limiting behavior of the model and approximation are qualitatively different in the subcritical case, requiring a novel analysis of the limiting behavior of the Hamiltonian system away from its deterministic subsystem. This yields a novel, general technique of approximation of the quasistationary distribution of stochastic epidemic and birth-death models and may lead to techniques for analysis of these models beyond the quasistationary distribution. For a classic SIS model, the approximation found for the quasistationary distribution is very similar to published approximations but not identical. For a birth-death process without depletion of susceptibles, the approximation is exact. Dynamics on the phase plane similar to those predicted by the Hamiltonian analysis are demonstrated in cross-sectional data from trachoma treatment trials in Ethiopia, in which declining prevalences are consistent with subcritical epidemic dynamics.

1. Introduction

Stochastic models are a common tool in epidemiological research, where public health interventions aim at the reduction of fluctuating counts of infected or infective individuals [1], and models are used in explaining, predicting, and responding to acute and chronic diseases of public health significance.

A fundamental result is the presence of a critical value of the basic reproduction number R0, defined as the expected number of secondary cases resulting from a single infective case in an otherwise susceptible population. Supercritical diseases, those with R0 > 1, tend to stabilize around a positive number of infective cases that can persist for very long times, while in subcritical cases (R0 < 1) the infective count declines to zero on a relatively short timescale. In either case, the long-term, stationary probability distribution of number of infective cases is trivial, as all epidemics in finite population stochastic transmission models must eventually die out due to chance fluctuations, but the quasistationary distribution—the distribution conditional on nonextinction of the disease—can be very informative about the behavior of the system within finite time intervals.

When R0 < 1, the quasistationary distribution of number of infective cases in simple transmission models is often approximately geometric, with probability of I infective cases proportional to (R0) [2, 3]. Prevalences consistent with the geometric distribution, when analyzed statistically across multiple locations simultaneously, have been observed in trachoma elimination trials at times in which the disease's dynamics are subcritical [4-6].

Such statistics of case count distributions observed in multiple communities at a single time may be able to help provide an assessment of the dynamics of a disease, possibly of its basic reproductive number and, hence, of the future time course of the disease. An approximately geometric distribution of prevalences also implies that there will be more high-prevalence communities than there would be in a lighter-tailed distribution, even when the mean prevalence is low and declining. This suggests that an exceptionally high-prevalence community may be simply a statistical outlier, which can be expected to regress to the mean without intervention, rather than a “transmission hotspot” calling for intensified intervention [5].

While the quasistationary distribution of a specific stochastic model can be calculated as an eigenvector of a Markov transition matrix, since the equations for the entries of that vector cannot be solved explicitly for even very simple models, research has focused on approximations ([2, 7–9], e.g.). Barbour and Pollett [10] established that the quasistationary distribution is a fixed point of a given map defined on probability mass functions, allowing efficient approximation techniques [11]. The fixed point of that map can also be found using a “ratio of means” approach built on waiting times rather than transition rates [12] that can aid in calculation. Quasistationary approximations for diffusion processes and branching processes are also well developed and are the subject of active research and development [3, 11, 13].

In this paper we introduce a method of approximating the quasistationary distribution of a stochastic model in the subcritical regime, using a technique that has been used previously to approximate rare large-deviation events in supercritical dynamics [14-16]. This technique takes a large-population limit of the model dynamics in a way that yields a Hamilton-Jacobi equation, which can be understood by analyzing the geometry of an associated Hamiltonian ODE system.

This Hamiltonian approach to stochastic mechanics, innovated by Graham and Tél [17] for diffusion equations and extended by Hu [18] to master equations, has primarily been used to study stationary solutions of the limiting stochastic process, by locating special solutions of the Hamiltonian ODE system, characterized by H = 0, where H is the Hamiltonian. The Hamiltonian ODE system includes the deterministic limit of the stochastic model as an invariant subsystem within the equipotential (H = 0) set, and at each limit set of the deterministic system, the equipotential set extends outwards into the nondeterministic regions of the Hamiltonian system's phase space. Those extensions reveal quantitative information about the system's stochastic behavior near attractors. Thus they are used to analyze stationary probability densities associated with attractors and other limit sets of the deterministic system and the frequencies and paths of rare escape events from one attractor to another [15, 19–22]. This geometric structure, which encodes characteristics of the deterministic limit of the stochastic system and the probability distribution of deviations from the deterministic limit, is strange in comparison to the structures seen in Hamiltonian systems from physics and is much less well understood.

Here we investigate the use of structures within the equipotential set, but at a distance from the deterministic subsystem, to analyze a stochastic model's behavior. We identify such a structure far from the deterministic subsystem with the quasistationary behavior of an epidemic model, in contrast to the use of structures intersecting the deterministic subsystem to analyze stationary behavior.

2. Limiting Behavior of Birth-Death Process

Many models of stochastic epidemic dynamics, biological population dynamics more generally, and branching processes are included in the category of birth-death processes. Here we apply the analysis of Hu [18] to this class of processes, and below we will apply it to specific example models.

A stochastic birth-death process models the size of a single population, altered by events in which the size either increases by one or decreases by one. The rate of increase from size k is labeled B(k) and the rate of decrease from size k is labeled D(k). Writing P(k, t) for the probability that the size is k at time t, the change in probability over time is governed by a master equation: Taking D(0) = 0 and B(−1)P(−1, t) = 0 for all t, the dynamics of the master equation is confined to nonnegative values of k. In order to take a large-system-size limit, let Ω be a measure of system size such as, for example, a maximum population size, such that, as we consider increasingly large birth-death systems in which both Ω and k become unboundedly large, the ratio k/Ω remains finite. For example, in a system with finite population size N, we can use Ω = N, as we will see below. Then letting x = k/Ω, we obtain a transformed master equation where b(x) = (1/Ω)B(Ωx) and d(x) = (1/Ω)D(Ωx). Let the functions b and d be smooth functions of x for each Ω, with a smooth limit as Ω → ∞.

Additionally, let ϕ(x, t) be a probability density function that is smooth in x and t, such that ϕ(k/Ω, t) = ΩP(k/Ω, t). Following Hu [18], this allows construction of a Kramers-Moyal expansion of the dynamics, by substituting and Taylor expanding the master equation around x so that it is expressed using only values at x:

To derive a partial differential equation in the large-system limit, we rewrite the density as an exponential expression: Assume that the function U can be expanded in powers of Ω on 0 < x < 1, and that the terms of that expansion other than u(x, t) vanish asymptotically as Ω approaches infinity. This ansatz, known as the WKB approximation [18, 26], makes it possible to generate a partial differential equation in u.

With these assumptions, derivatives of products of ϕ take on a simplified form, Substituting, the expansion of (3) to first order is Thus, in the large size limit, (3) becomes a partial differential equation for u:

2.1. The Associated Hamiltonian System

Because the right hand side of (8) contains only first partial derivatives of u, it has the form of a Hamilton-Jacobi equation of classical mechanics [27], with the consequence that it can be analyzed using characteristic curves described by an associated system of ordinary differential equations [18]. This analysis is based on the Hamiltonian function

From that Hamiltonian a two-dimensional dynamical system can be written, whose state variables are x, the scaled population size, and a conjugate variable p, which takes the place of ∂u/∂x in the Hamiltonian. The associated Hamiltonian dynamical system is Trajectories of this system do not correspond to realizations of the stochastic birth-death process but rather trace out curves along the surface of u versus x and t, which can be used to analyze the behavior of u over time.

Thus we can gain information about birth-death processes in the large size limit by using this associated system to analyze the Hamilton-Jacobi equation (8). Stationary solutions of the master equation, characterized by the equilibrium condition dϕ(x, t)/dt = 0, are identified with curves on the (x, p) plane on which H(x, p) = 0.

In the case of this one-dimensional system, though not in the general master equation case, the Hamiltonian has two factors, which contribute two solution sets to the solution of H = 0.

The flat subspace p = 0 is always a solution set for H = 0 in Hamiltonian systems constructed from master equations in this way [18]. The dynamics within this set are the dynamics of the ODE approximation to the stochastic dynamics, and fixed points and other limit sets of the Hamiltonian system located in this set correspond to fixed points and other limit sets of this deterministic subsystem. Other solutions to the equation H = 0 pass transversely through those limit sets and can reveal information about the stochastic behavior of the master equation system, as we will see in the treatment of the supercritical SIS model, below.

In the birth-death systems we consider here, in which k = 0 is an absorbing state, a common factor of x can be taken out of b(x) and d(x), allowing us to describe three components of the solution set in all.

3. The SIS Model

The SIS (susceptible-infective-susceptible) model provides a simple representation of infectious disease processes in the absence of immunity [28]. Classically, this model describes the number of susceptibles S and infective cases I in a population of fixed size, where increase in the infective class is driven by infective-susceptible contact events, and infective cases return to the susceptible class at a rate independent of contact with others. SIS models have been used to describe a range of diseases, including trachoma [29] and sexually transmitted infections [30]. In population biology, a model identical in form to this one is known as a stochastic logistic model [31].

In the basic SIS model, the infective class increases at a rate βS(I/N), which is proportional to a quadratic susceptible-infective contact rate, and decreases at a per capita constant rate γI, with S = N − I, and total population N held fixed. Thus it is the number of infective cases, I, that is the stochastically varying state variable of the model. Infective cases are added by transmission events, at rate β(S/N)I, where β is the transmission rate per susceptible-infective pair [1]. Cases return to the susceptible class at rate γI, where γ is the per capita removal rate. The parameters can be combined into one nondimensional value by rescaling the time variable by a factor of γ, after which the birth and death rates are where R0 = β/γ is the basic reproduction number [28].

Using system size Ω = N, the analysis we have presented for birth-death systems applies to the SIS model, with Hamiltonian where x = I/N is the infective fraction of the population.

3.1. The Supercritical Case

In the supercritical (R0 > 1) case, the SIS process is attracted to a positive, or endemic, equilibrium value x = 1 − 1/R0, at which the birth and death rates are equal. The probability density of the fraction infective case concentrates around that value. On very long time scales, however, in finite systems, stochastic fluctuation will bring the fraction infective case to zero, which is an absorbing state from which the epidemic cannot return. Thus the stationary distribution of the process is a point mass at x = 0, and the density function concentrated around the endemic equilibrium, while it is a stationary distribution in the infinite-size limit and is the quasistationary distribution in the finite cases.

The Hamiltonian analysis of the supercritical SIS model has been treated exactly elsewhere [16, 20]. The phase plane of the Hamiltonian system is shown in Figure 1.

Figure 1

Phase plane of the Hamiltonian dynamical system ((R0 = 2). Arrows depict the flow of the dynamics of x and p. The three invariant curves of the dynamics (solution curves of H = 0) are shown in gray: the two axes of the space and one nontrivial curve. The nontrivial curve corresponds to the quasistationary solution of the stochastic SIS model, as discussed in the text.

Stationary solutions of the PDE correspond to solutions of H(x, p) = 0 on this plane, when p is interpreted as ∂u/∂x. The Hamiltonian factors into three parts: which directly identifies the three solution curves of H = 0 in the plane: two trivial solutions,and one nontrivial solution, shown in Figure 1. These curves are trajectories of the Hamiltonian dynamical system (11).

The horizontal axis of the phase plane, which is the p = 0 solution, is isomorphic to the deterministic SIS system. Two of the fixed points of the Hamiltonian system are the fixed points of that deterministic system—the disease-free equilibrium at (0,0) and the endemic equilibrium at (1 − 1/R0, 0). They are located at the points where the horizontal axis intersects the other two solution curves. A third fixed point, at (0, −ln⁡R0), also corresponds to the disease-free state (x = 0) but is at the intersection of solution curves away from the horizontal axis.

The nontrivial solution curve (17) corresponds to the stationary solution of u(x) on which probability concentrates around the endemic equilibrium, and the fixed points on it describe the probability density at the endemic and disease-free equilibria. That solution is a function u(x) that solves Changing variables to s = 1 − x and integrating produce a closed-form solution, This provides a closed-form solution for the quasistationary probability density: The constant C1 is determined by the constraint that ∫01ϕ(s)ds = 1.

In supercritical models in general, the equipotential surfaces (solutions of H = 0) near the nontrivial solution of the deterministic subsystem describe the behavior of the probability distribution of rare events, which are located in the tail of the stationary distribution.

The above stationary solution approximates the quasistationary density in the finite-N SIS system, in which extinction is a rare event given large N.

It provides an approximation for the time to extinction in the stochastic dynamics. The function u is the action of classical mechanics. The most probable path to extinction can be obtained by maximizing the function u(x), which produces the equipotential surfaces H = 0. The path is explicitly calculated by integrating along the H = 0 curves, both in this SIS case and in more complex models (e.g., [16]).

4. Subcritical Dynamics

In the deterministic SIS system in the subcritical case, x relaxes to zero for all initial conditions 0 ≤ x ≤ 1. The master equation solution also relaxes to x = 0, with probability mass declining to zero at all other values of x [2]. In this case, the quasistationary distribution is not stationary even in the large-N limit due to the deterministic attraction of the origin. The WKB hypothesis that the probability current near the absorbing state x = 0 vanishes when the system size N grows without bound is not satisfied, and we do not use the stationary behavior of the PDE (which relaxes to a point mass) to analyze the quasistationary behavior of the master equations. Instead we use the transient behavior of the PDE to identify the equilibrium structure in the Hamiltonian phase plane that describes the master equation's quasistationary solution.

4.1. Using the Phase Plane to Analyze Dynamics of the Hamilton-Jacobi Equation

In the Hamiltonian phase plane for the subcritical model, the same three solution curves for H = 0 are present as in the supercritical case, but they fall in different places on the phase plane, as shown in Figure 2. In this case, the point of intersection of the nontrivial curve (17) and the horizontal axis is shifted to the left of the origin. The endemic equilibrium represented by that point is lost in a transcritical bifurcation when R0 declines below 1, and the origin becomes the attracting solution for the stochastic SIS system. The intercept where the nontrivial curve (17) meets the vertical axis, at p = −ln⁡R0, is now above p = 0.

Figure 2

Phase plane of Hamiltonian dynamical system for subcritical SIS system (R0 = 0.5). Flow is represented by arrows and the three invariant curves of the dynamics (solution curves of H = 0) are shown in gray, as in Figure 1. In this case, the nontrivial curve is shifted to a different position, and its intersections with the axes are located above and to the left of the origin, where in the supercritical case (Figure 1) they are below and to the right of the origin. This leads to qualitatively different dynamics, requiring a different analysis to explain the quasistationary behavior of the model.

Because of this bifurcation, in the subcritical case we cannot apply the analysis used for the supercritical case, as the system is drawn to a singular value of x at which the H = 0 curve crossing the horizontal axis is vertical and cannot be translated to values of ∂u/∂x as a function of x. To study the quasistationary distribution of this system requires further analysis.

Any smooth initial distribution ϕ(x) can be mapped onto a curve in the (x, p) plane on which p = ∂u/∂x at every value of x, where u is defined by ϕ(x) = Ne− as above. This curve for an example initial distribution is plotted in Figure 3.

Figure 3

Initial condition for the subcritical SIS system on the Hamiltonian phase plane, represented by a curve of p values as a function of x. In this and following figures, the initial condition used is a β distribution with α = β = 2 (i.e., ϕ0(x) = 6x(1 − x)) and using N = 100, transformed to a curve in the x-p plane using the relations u(x) = −ln⁡(ϕ0(x)/N)/N and p = ∂u/∂x.

Integrating points of this curve forward along trajectories of this system produces a geometric representation of the time evolution of the system as a moving curve in the phase plane, on which the changing shape of ∂u/∂x is visible, and that relation between ∂u/∂x and x provides information about the form of the function u(x).

In terms of Hamiltonian dynamics, the function u(x, t) is the action of the system, a scalar quantity that can be evaluated by integrating along its trajectories:

For convenience, it is possible to calculate u directly when integrating the Hamiltonian dynamics numerically, by extending the dynamical system to include u as a state variable:Integrating this system, with initial conditions u(x, 0) = u0(x) at selected points of the initial curve, then yields values of u(x, t) explicitly for positive t.

4.2. Evolution of the Subcritical System from Initial Conditions

As time passes, each point of the p-versus-x curve moves on the phase plane according to the Hamiltonian dynamics. Their evolution stretches and translates the curve across the phase plane, as shown in Figure 4. While any given point may move in somewhat strange ways, including many that tend to infinity in the upper right direction, the curve moves smoothly to the left, approaching the vertical line x = 0 and the gray curve that extends into the first quadrant.

Figure 4

Transient dynamics of the subcritical SIS system on the Hamiltonian phase plane, evolving from the initial condition depicted in Figure 3 (red) toward later states (yellow, green, and blue), as each point of the initial curve moves according to the Hamiltonian dynamics (11).

From the moving points (x, p, u) of this curve, a plot of u versus x can be constructed, or of ϕ = Ne− versus x, at each time t. Figure 5 presents this plot of ϕ versus x in time. The peak of the probability density moves asymptotically toward x = 0, and there is a declining tail to the right of the peak.

Figure 5

Transient dynamics of probability density in the subcritical SIS system, displayed as ϕ(x, t) = Ce− versus x using the same data points as in Figure 4, with N = 100. Each curve is normalized to total probability one. The quasistationary distribution (23) is plotted in gray.

A number of features of the evolution of u(x, t) versus x are visible in this view of the dynamics. As discussed above, the dynamics on the horizontal axis of the phase plane is identical to the usual deterministic ODE for the SIS system. When p is read as ∂u/∂x, it follows that horizontal axis, where p = 0, corresponds to the extrema of the potential function u(x, t) with respect to x. In the case pictured in these figures, the only extremum is a minimum of u(x, t), which is a maximum of ϕ(x, t). This implies that the maximum point of the probability density function ϕ, which is the mode of the probability distribution, in the large-system approximation we are using (8), moves in exact accordance with the deterministic SIS dynamics.

Regions of x values for which a curve in the x-p plane is below the horizontal axis are regions where ∂u/∂x < 0 and equivalently on which ϕ(x, t) is increasing in x, and regions where the curve is above the axis are where ϕ(x, t) is decreasing in x. Near the vertical axis, the p-versus-x curve diverges to p = −∞. The fact that p, representing ∂u/∂x, becomes negatively infinite there strongly suggests that u(x) is divergent to +∞ at x = 0, and so limϕ(x, t) = 0, at least in cases like the one illustrated in which ϕ(0) is zero in the initial conditions.

If the Hamilton-Jacobi PDE (8) is used to approximate any finite-N system, by grouping the probability density into bins of width 1/N, the result will be that probability mass accumulates in the bin that includes x = 0, and all the other bins contain a tail that is decreasing in x, and whose total mass declines asymptotically to zero as t → ∞.

Figure 4 demonstrates that, in the long term, the p-versus-x curve becomes asymptotically close to the union of the vertical axis below the positive-p equilibrium and the nontrivial H = 0 curve (17) at and above that equilibrium. We conclude that as the probability density accumulates near x = 0, the shape of the tail of the density on x > 0 approaches a function described by the diagonal curve, which is the nontrivial solution (17) of H = 0. That tail defines the conditional distribution of x given x > 0, and therefore the limiting curve (17) should provide an approximation for the quasistationary distribution of the SIS master equations.

4.3. Explicit Approximation for the Quasistationary Distribution

From the above analysis we conclude that the quasistationary probability density function of the master equation system (1) is approximated by the density function represented by the nontrivial H = 0 curve (17). This is solved in the same way as in the supercritical case: where s = 1 − x.

While in the supercritical case this density function has a mode at the endemic value s = 1/R0, in this case the density is greatest at x = 0 (s = 1), as the function is monotonic decreasing on the interval 0 < x < 1.

Changing variables back to the number of infective cases, I = Nx = N(1 − s), the quasistationary approximation becomes using the appropriate normalizing factor C2 for this discrete probability mass function.

This quasistationary approximation is closely related to the classical approximation p(1) of Kryscio et al. [2, 8] (see also [32]): their approximation, when transformed using Stirling's approximation for factorials, yields the approximation we have derived: (where C3, C4 are normalizing constants).

Previous approximations and numeric evaluation have established [2, 7, 8] that the quasistationary distribution of the subcritical SIS system is approximately geometric near I = 0, with the probabilities of successive values of I having ratio R0.

Thus the approximating geometric distribution has the form The geometric distribution is characterized by the constant slope of its logarithm:

Comparing to our approximation p, the slope of ln⁡p is not constant:

However, near I = 0, the nonconstant term is approximately zero, and the slope of the logarithm is approximately ln⁡R0, with the consequence that the distribution is approximately geometric with the desired ratio when I ≪ N.

Since the ratio (N − I)/N is smaller than one when 0 < I < N and thus its logarithm is negative, it follows that the probability mass function p decreases to zero more rapidly than the geometric function Γ does as I increases.

In an appendix we compare the SIS process to a birth-death process that has the transmission and removal rates of the SIS model without the effect of depletion of susceptibles and whose quasistationary distribution is exactly the geometric distribution that approximates the above distribution. The phase plane analysis of the birth-death process provides visual evidence that the parameter characterizing the approximating geometric distribution by its rate of decay is determined by the intercept where the nontrivial curve (17) crosses the vertical axis.

5. Application of SIS Model Analysis to Trachoma Case Counts

Trachoma is a common subclinical childhood infection in certain regions of the less-developed world. Repeated infection results in scarring of the eyelid and trichiasis (turning inward of the eyelashes, so that the eyelids scrape against the cornea). Millions of cases of blindness have resulted. The causative agent, Chlamydia trachomatis, can be cleared with high efficacy with a single dose of azithromycin [33]. The World Health Organization currently recommends annual mass treatment in affected communities as a public health control measure [33, 34].

During a clinical trial of timing of mass administration of azithromycin in the Amhara Region of Ethiopia [23, 34, 35], village-level prevalence data were collected. At baseline the probability distribution of village-level prevalences, omitting zero values, had a mean of 0.39 (range 0.08–0.62) (Figure 6, top plot). After the initiation of mass treatment at or exceeding recommended WHO levels, the mean prevalence declined, and the distributions became indistinguishable from exponential [5] (Figure 6, subsequent plots). This finding is consistent with the approximately exponential distributions predicted by simple epidemic models, as discussed above. The matter is of more than theoretical interest, as mentioned in our introduction: the long tail of the exponential distribution implies that, during an elimination campaign, some communities may have unexpectedly large prevalence and appear to be outliers when in fact they are entirely consistent with the variation expected.

Figure 6

Changing trachoma prevalence at baseline and at 6-month intervals during the TANA trial of mass administration of azithromycin [23]. As the trial progresses, the prevalences become smaller and become more closely approximated by the exponential [5]. (Individual village prevalences are shown in tick marks on the horizontal axis. Curves result from beta distribution kernel density smoothing [24], with smoothing parameter determined from leave-one-out cross-validation [25].)

The SIS model has been used in practice to assess treatment frequency needed for elimination of trachoma [4, 29, 36].

Figure 7 displays these probability density functions ϕ(x) transformed to the phase plane representation defined above, p(x) = −(d/dx)ln⁡(ϕ(x)/N)/N. We assume a population size N = 100 per village, which is approximately the number of children at risk in one of these villages [23]. In this plot, the same motion from lower right to upper left is visible, with convergence to the vertical axis and possibly to a curve leaving that axis in the positive quadrant. More abundant data may permit location of such a limiting curve that would intersect the vertical axis in this representation of the data. That curve would provide an estimate of the quasistationary behavior of the disease, and its intercept would provide an estimate of the disease's R0.

Figure 7

Phase plane representation of changing trachoma prevalence data from TANA trial shown in Figure 6. Each curve on the plot corresponds to one of the distributions shown in Figure 6, transformed to the x-p plane as in earlier figures (see text for details). Over time the curves shift upward and to the left, moving close to the vertical axis for smaller values of p and diverging from it at larger values of p, similar to the motion seen in the Hamiltonian analysis of the SIS model (Figure 4). Each curve in this figure is restricted to the range of the nonzero prevalence values.

6. Summary

Hamiltonian structures describing master equation and diffusion equation systems are the subject of ongoing exploration in stochastic processes research, where the solution sets of H = 0 near the deterministic subspace are used to model quasistationary behaviors and rare transition events, such as switching between states or noise-induced extinctions. We have presented an application of these structures far away from the deterministic subsystem, to approximate the probability distribution of a process near an absorbing singular point, where the WKB hypothesis does not hold and transient dynamics of the limiting PDE rather than its large-time limit behavior must be used to identify the structure corresponding to the quasistationary probability distribution of the finite-size system.

Quasistationary solutions in epidemic models can generally not be solved exactly, so approximation techniques are crucial in analysis of these processes. We present an alternative approach to this approximation problem, which may be extensible to other similar model settings and whose full usefulness is yet to be discovered. The WKB approximation and the Hamiltonian and Lagrangian techniques of analysis that it makes available are powerful and flexible and may have applications in subcritical disease settings that go well beyond the quasistationary distribution.

Our exploration of cross-sectional prevalence data from trachoma trials, when the prevalence distributions are represented as curves on the Hamiltonian phase plane, reveals a pattern of motion consistent with the motion on the phase plane predicted by this analysis for a subcritical transmission model. Thus it is consistent, at least qualitatively, with a hypothesis that trachoma transmission in that trial setting is in fact subcritical and stochastic. This analysis fails to disconfirm that hypothesis, though other explanations are possible. In epidemiological settings where more data are available, it may become possible to observe an upper limiting curve in such a plot as well as the convergence to the vertical axis. By revealing an emerging shape of the tail of the prevalence distribution, information about that curve could contribute to description of the quasistationary behavior of the disease. Such information also may contribute to an estimate of its basic reproduction number, arrived at independently of any estimate based on temporal change in prevalences.

Beyond the one-variable birth-death models that we have analyzed, the techniques that we explore here for study of quasistationary dynamics may be of use with models with more stages of disease progression or differing transition rates, multitype models, models with patch or network structure (cf. [37]), and other cases that are more complex than the simple models presented here. In population biology, the SIS model we have discussed is also known as a stochastic logistic model [38], and this analysis has promise for population biology models that are similar but not identical to this model. While the primary goal in conservation biology is to preserve the populations in question, rather than to eradicate them as in epidemiology, declining populations are clearly of interest and the models in use may benefit from a similar analysis. This analysis may be of use in other applications as well, where quasistationary dynamics near an absorbing state is of interest.