Abstract
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 , defined as the expected number of secondary cases resulting from a single infective case in an otherwise susceptible population. Supercritical diseases, those with , tend to stabilize around a positive number of infective cases that can persist for very long times, while in subcritical cases () 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 , the quasistationary distribution of number of infective cases in simple transmission models is often approximately geometric, with probability of infective cases proportional to [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 , where is the Hamiltonian. The Hamiltonian ODE system includes the deterministic limit of the stochastic model as an invariant subsystem within the equipotential () 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 is labeled and the rate of decrease from size is labeled . Writing for the probability that the size is at time , the change in probability over time is governed by a master equation: Taking and for all , the dynamics of the master equation is confined to nonnegative values of . 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 become unboundedly large, the ratio remains finite. For example, in a system with finite population size , we can use , as we will see below. Then letting , we obtain a transformed master equation where and . Let the functions and be smooth functions of for each , with a smooth limit as .
Additionally, let be a probability density function that is smooth in and , such that . Following Hu [18], this allows construction of a Kramers-Moyal expansion of the dynamics, by substituting and Taylor expanding the master equation around so that it is expressed using only values at :
To derive a partial differential equation in the large-system limit, we rewrite the density as an exponential expression: Assume that the function can be expanded in powers of on , and that the terms of that expansion other than vanish asymptotically as approaches infinity. This ansatz, known as the WKB approximation [18, 26], makes it possible to generate a partial differential equation in .
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 :
2.1. The Associated Hamiltonian System
Because the right hand side of (8) contains only first partial derivatives of , 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 , the scaled population size, and a conjugate variable , which takes the place of 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 versus and , which can be used to analyze the behavior of 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 , are identified with curves on the plane on which .
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 .
The flat subspace is always a solution set for 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 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 is an absorbing state, a common factor of can be taken out of and , 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 and infective cases 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 , which is proportional to a quadratic susceptible-infective contact rate, and decreases at a per capita constant rate , with , and total population held fixed. Thus it is the number of infective cases, , that is the stochastically varying state variable of the model. Infective cases are added by transmission events, at rate , where is the transmission rate per susceptible-infective pair [1]. Cases return to the susceptible class at rate , 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 is the basic reproduction number [28].
Using system size , the analysis we have presented for birth-death systems applies to the SIS model, with Hamiltonian where is the infective fraction of the population.
3.1. The Supercritical Case
In the supercritical () case, the SIS process is attracted to a positive, or endemic, equilibrium value , 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 , 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.
Stationary solutions of the PDE correspond to solutions of on this plane, when is interpreted as . The Hamiltonian factors into three parts: which directly identifies the three solution curves of 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 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 and the endemic equilibrium at . They are located at the points where the horizontal axis intersects the other two solution curves. A third fixed point, at , also corresponds to the disease-free state () but is at the intersection of solution curves away from the horizontal axis.
The nontrivial solution curve (17) corresponds to the stationary solution of 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 that solves Changing variables to and integrating produce a closed-form solution, This provides a closed-form solution for the quasistationary probability density: The constant is determined by the constraint that .
In supercritical models in general, the equipotential surfaces (solutions of ) 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- SIS system, in which extinction is a rare event given large .
It provides an approximation for the time to extinction in the stochastic dynamics. The function is the action of classical mechanics. The most probable path to extinction can be obtained by maximizing the function , which produces the equipotential surfaces . The path is explicitly calculated by integrating along the 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, relaxes to zero for all initial conditions . The master equation solution also relaxes to , with probability mass declining to zero at all other values of [2]. In this case, the quasistationary distribution is not stationary even in the large- limit due to the deterministic attraction of the origin. The WKB hypothesis that the probability current near the absorbing state vanishes when the system size 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 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 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 , is now above .
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 at which the curve crossing the horizontal axis is vertical and cannot be translated to values of as a function of . To study the quasistationary distribution of this system requires further analysis.
Any smooth initial distribution can be mapped onto a curve in the plane on which at every value of , where is defined by as above. This curve for an example initial distribution is plotted in Figure 3.
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 is visible, and that relation between and provides information about the form of the function .
In terms of Hamiltonian dynamics, the function 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 directly when integrating the Hamiltonian dynamics numerically, by extending the dynamical system to include as a state variable:Integrating this system, with initial conditions at selected points of the initial curve, then yields values of explicitly for positive .
4.2. Evolution of the Subcritical System from Initial Conditions
As time passes, each point of the -versus- 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 and the gray curve that extends into the first quadrant.
From the moving points of this curve, a plot of versus can be constructed, or of versus , at each time . Figure 5 presents this plot of versus in time. The peak of the probability density moves asymptotically toward , and there is a declining tail to the right of the peak.
A number of features of the evolution of versus 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 is read as , it follows that horizontal axis, where , corresponds to the extrema of the potential function with respect to . In the case pictured in these figures, the only extremum is a minimum of , which is a maximum of . 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 values for which a curve in the - plane is below the horizontal axis are regions where and equivalently on which is increasing in , and regions where the curve is above the axis are where is decreasing in . Near the vertical axis, the -versus- curve diverges to . The fact that , representing , becomes negatively infinite there strongly suggests that is divergent to at , and so , at least in cases like the one illustrated in which is zero in the initial conditions.
If the Hamilton-Jacobi PDE (8) is used to approximate any finite- system, by grouping the probability density into bins of width , the result will be that probability mass accumulates in the bin that includes , and all the other bins contain a tail that is decreasing in , and whose total mass declines asymptotically to zero as .
Figure 4 demonstrates that, in the long term, the -versus- curve becomes asymptotically close to the union of the vertical axis below the positive- equilibrium and the nontrivial curve (17) at and above that equilibrium. We conclude that as the probability density accumulates near , the shape of the tail of the density on approaches a function described by the diagonal curve, which is the nontrivial solution (17) of . That tail defines the conditional distribution of given , 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 curve (17). This is solved in the same way as in the supercritical case: where .
While in the supercritical case this density function has a mode at the endemic value , in this case the density is greatest at (), as the function is monotonic decreasing on the interval .
Changing variables back to the number of infective cases, , the quasistationary approximation becomes using the appropriate normalizing factor for this discrete probability mass function.
This quasistationary approximation is closely related to the classical approximation 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 , 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 , with the probabilities of successive values of having ratio .
Thus the approximating geometric distribution has the form The geometric distribution is characterized by the constant slope of its logarithm:
Comparing to our approximation , the slope of is not constant:
However, near , the nonconstant term is approximately zero, and the slope of the logarithm is approximately , with the consequence that the distribution is approximately geometric with the desired ratio when .
Since the ratio is smaller than one when and thus its logarithm is negative, it follows that the probability mass function decreases to zero more rapidly than the geometric function does as 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.
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 transformed to the phase plane representation defined above, . We assume a population size 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 .
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 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.
Appendix
Comparison to Poisson Birth-Death Process
The close approximation to the geometric by the SIS and other transmission models when infective counts are small can also be explained by comparing the transmission model to a Poisson birth-death process, that is, a process in which depletion of the susceptible class is not accounted for [1]. The quasistationary limit of this process is the Yaglom limit of its associated branching process [39, 40] and is exactly geometric.
Here we use the above Hamiltonian phase plane analysis to approximate the quasistationary limit of a Poisson birth-death process with the same basic reproduction number , for comparison to the above results.
In the Poisson birth-death process, the birth and death rates are the same as in the SIS model, except for the absence of the nonlinear factor in the birth rate: It follows that the resulting Hamiltonian also is the same except for the absence of that factor:
As in the SIS case, the Hamiltonian factors into three parts, corresponding to three intersecting components of the solution set of . Again, two components are the vertical and horizontal axis. The third, nontrivial solution in this case is a horizontal line rather than a rising curve. Unlike the SIS case, here the nontrivial curve does not intersect the horizontal subsystem (except in the critical case, which we will leave aside in this discussion).
In the supercritical case (Figure 8), the dynamics on the horizontal axis (the deterministic subsystem) is similar to the SIS model in that positive values rise away from zero, but with the difference that they increase to infinity rather than to a finite endemic equilibrium value. In the subcritical case (Figure 9), in which the birth-death process tends to extinction, the behavior of the Hamiltonian system in the half-plane is qualitatively the same as for the subcritical SIS. The only difference is in the form of the nontrivial curve that specifies the quasistationary distribution.
The quasistationary curve is specified by the equation , or . Substituting for produces the quasistationary solution for the action, , and for the probability density, . This is equivalent to the geometric distribution with parameter , which is well known to be the quasistationary distribution of this process.
We note that we can relate the geometric approximation to the geometry of the Hamiltonian phase plane. As discussed above, the geometric approximation to the quasistationary distribution is characterized by the logarithmic slope . The density function in terms of the fraction is , and the action is defined by , or . Combining, so that But note that, on the phase plane, the vertical coordinate is identified with ; so it follows that the parameter of the geometric distribution approximating the quasistationary distribution near is revealed by the value of the nontrivial solution for at , that is, the intercept where the limiting curve describing the quasistationary distribution crosses the vertical axis. That intercept is equal to , where is the parameter of the approximating geometric distribution.
We note that the horizontal line found in the Poisson birth-death model’s phase plane and the SIS process’s limiting curve converge on the same value at , confirming visually that this birth-death process is a good approximation for the SIS process when its infective count is much smaller than .
Conflicts of Interest
The authors declare they each have no conflicts of interest.
Acknowledgments
This study was supported by a Models of Infectious Disease Agent Study (MIDAS) grant from the US NIH/NIGMS to the University of California, San Francisco (U01GM087728), by US NEI R01-EY025350 and by a Research to Prevent Blindness Award. Ira B. Schwartz was supported by NRL base funding (N0001414WX00023) and Office of Naval Research (N0001414WX20610).