SIAM News Blog

New Results in Stochastic Analysis Using Dynamical Systems Theory

By Lora Billings and Eric Forgoston

With its long history in a variety of fields, stochastic analysis continues to challenge modelers. Noise plays an important functional role in changing a system’s dynamics in applications across mathematical finance, physics, and biology [7]. Researchers in multiple disciplines have developed numerous approaches and notations to understand the dynamics beyond numerical simulation. Adding to the complexity, one can model noise as an external source or a collection of random internal processes, and few have considered incorporating either concept into models with time-dependent forcing. Inspired by the work in [2, 4], we find that classical tools from dynamical systems and perturbation theory provide an intuitive way to analyze stochastic epidemiological models with seasonal forcing and predict noise-induced rare events.

Seasonal forcing is important in epidemic modeling because contact rates and other parameters may vary in time in diseases like measles, mumps, and chickenpox. This is evident from the periodic and aperiodic recurrence of such conditions, as observed in pre-vaccine United States data and English and Welsh data [1]. Our recent work considers a general class of dynamical models with internal noise that symbolizes demographic stochasticity, or random interactions of individuals [3]. Our model can represent any group of well-mixed individuals from the nano- to macro-scale. We study disease dynamics in a human population that requires seasonality to accurately capture the force of infection, and focus on internal random dynamics that can induce a rare, spontaneous disease eradication.

To understand random dynamics, one can use a master equation that describes the probability of observing a given number of individuals in a population at a certain time. As the master equation generally cannot be solved analytically, we use a Wentzel-Kramers-Brillouin-Jeffreys (WKBJ), or eikonal, approximation to the master equation. This ansatz approximates the probability density function for the system dynamics and leads to the development of a Hamiltonian system. The Hamiltonian’s dimensions are twice the dimensions of the original system’s mean-field equation due to the conjugate momenta variables, which capture the effects of noise. The benefit is that the transformed system is deterministic and lends itself to dynamical systems theory. For instance, standard methods can identify steady states, which correspond to the metastable states in the original stochastic system.

As an example, consider a time-dependent susceptible-infected-susceptible (SIS) model consisting of a finite population of \(N\) individuals divided into two groups: susceptibles and infectives. The flux terms for these groups are as follows. An individual is born susceptible, assuming a constant birth rate. Through contact modeled by a mass action approximation, he/she may become ill and be classified as infectious. We suppose that the contact rate varies throughout the year in a periodic fashion, approximated by the time-dependent function \(\beta(t)=\beta_0(1+\delta\cos(2\pi t))\), where \(t\) is time, \(\beta_0\) is the average contact rate, and \(\delta\) is a small parameter. Assuming an average recovery rate, the individual recovers and returns to the susceptible group. While removal by natural death is possible from both groups, we presume the absence of disease-related deaths in this model.

Figure 1. A realization of the stochastic SIS model with seasonality exhibiting an extinction event. We plot the number of infectives (I) in time (t) in blue. Note how the graph follows the red curve, which represents the endemic limit cycle given by the solution to the mean field equation (p=0, no noise). Image courtesy of [3].

Figure 1 offers a numerical realization of this stochastic time-dependent SIS model for parameters that support the disease’s persistence. We use a Monte Carlo algorithm to evolve the population in time [5]. In Figure 1, the population persists near what appears to be an oscillating steady state and eventually drops to zero via spontaneous extinction. Our goal is to develop a method to analytically predict the average time of occurrence for an extinction event. Numerous simulations reveal the existence of many possible escape paths, with extinction taking place at different times. Aggregation of simulation data demonstrates a path along which extinction is most likely to occur. This optimal path coincides with the manifold connecting steady states in the corresponding Hamiltonian system. Identifying the action along the optimal path lets researchers determine the mean time of escape from a steady state or the average time to an extinction event.

The analytical approach uses the WKBJ approximation \(\rho = e^{-N{\mathcal S}(i,t)}\) for the probability distribution in the master equation, for which \(\mathcal {S}\) is the action and the state variable \(i\) is the scaled infective population. The leading order of a Taylor series expansion of the scaled master equation results in a Hamilton-Jacobi equation with an effective Hamiltonian, for which we analyze the zero-energy solutions. The WKBJ approximation introduces a conjugate momentum variable \(p=\partial {\mathcal S}(i,t)/\partial i\), which doubles the system’s dimension but reveals additional dynamics that stem from stochasticity.

The advantage of recasting the problem in terms of the conjugate momentum is that it reveals a higher-dimensional Hamiltonian system with steady states classified as saddles and centers. This allows a trajectory to escape from the deterministically stable endemic limit cycle through the expanded \(p\) space (see Figure 2). Specifically, there exists a heteroclinic trajectory connecting it to the fluctuational extinction state for which \(i=0\) and \(p \ne 0\).

Figure 2. The phase space for Hamilton’s equations with forcing, which corresponds to the stochastic SIS model with seasonality. Note the location of the endemic and fluctuational extinction states and the connecting manifold. Approximation of the optimal path found with the iterative action minimizing method numerical scheme appears in red. The implicit ordinary differential equation solution is plotted in black. Image courtesy of [3].
If one rewrites the Hamiltonian using a linear expansion in \(\delta\), the Poincaré-Birkhoff fixed point theorem guarantees existence of the perturbed hyperbolic points. The presence of a corresponding heteroclinic trajectory depends on the optimal correction to the action, or the minimization of \(\mathcal {S}\) with respect to initial time \(t_0\). Finding the zeros of the derivative of \(\mathcal {S}\) yields an integral of a Poisson bracket of the unperturbed Hamiltonian and the time-dependent Hamiltonian evaluated on the unperturbed optimal path; this integral is the perturbed problem’s Melnikov function. A sufficient condition for the perturbed optimal path’s existence is the Melnikov function’s possession of simple zeros. We can numerically evaluate this equation to identify a periodic set of zeros, which determines the phase for the optimal correction to the action.

While we cannot explicitly solve the time-dependent Hamilton’s equations, we can numerically approximate the heteroclinic manifold via the optimal correction — or appropriate phase defined by initial time \(t_0\). One can use the implicit MATLAB solver ode15i to do so. We compared the results with those from the iterative action minimizing method [6], a numerical scheme based on Newton’s Method. Figure 2 depicts agreement of the time-dependent optimal paths resulting from the two different numerical methods. One can now use the paths to analyze the extinction event in Figure 1. Additional detail for this example is available in our published work [3].

Much work remains pertaining to the analysis of rare noise-induced transition events in time-dependent models. For epidemic models, such as the aforementioned SIS model, knowledge of the optimal path enables discovery of the mean time to extinction from the endemic state, which is verifiable by stochastic realizations. Determining the mean time to extinction allows researchers to quantify the effect of various control schemes—including vaccination pulses and treatment—on infectious disease. More generally, one can use classical dynamical systems theory to understand noise-induced transitions in a wide variety of stochastic problems, from biology to physics. Some examples include genetic switching, pest eradication, and bit error rate in communication systems. As scientists increasingly discover more systems that exhibit rare events, we look forward to the implementation of old and new techniques to understand the effect of stochasticity and control measures.

The authors presented this research at the 2017 SIAM Conference on Applied Dynamical Systems, which took place in Snowbird, Utah.

[1] Anderson, R.M., May, R.M., & Anderson, B. (1992). Infectious Diseases of Humans: Dynamics and Control. New York, NY: Oxford University Press.
[2] Assaf, M., Kamenev, A., & Meerson, B. (2008). Population extinction in a time-modulated environment. Phys. Rev. E, 78(4), 041123.
[3] Billings, L., & Forgoston, E. (2018). Seasonal forcing in stochastic epidemiology models. Ricerc. Mat., 67(1), 27-47.
[4] Dykman, M.I., Rabitz, H., Smelyanskiy, V.N., & Vugmeister, B.E. (1997). Resonant directed diffusion in nonadiabatically driven systems. Phys. Rev. Lett., 79, 118-1181.
[5] Gillespie, D.T. (1976). A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys., 22(4), 403-434.
[6] Lindley, B.S., & Schwartz, I.B. (2013). An iterative action minimizing method for computing optimal paths in stochastic dynamical systems. Phys. D: Nonlin. Phenom., 255, 22-30.
[7] Tsimring, L.S. (2014). Noise in biology. Rep. Prog. Phys., 77(2), 026601.

Lora Billings is dean of the College of Science and Mathematics at Montclair State University. Her research focuses on applied deterministic and stochastic dynamical systems that model applications in epidemiology, physics, and ecology. Eric Forgoston is an associate professor of applied mathematics at Montclair State University. His research involves the study of complex physical and biological phenomena, including material transport in the ocean, the outbreak and extinction of infectious diseases, behavior of biological and robotic swarms, food web dynamics in ecological systems, and the stability of fluid flows.

blog comments powered by Disqus