Extinction appears ubiquitously in many fields, including chemical reactions, population biology, evolution and epidemiology. Even though extinction as a random process is a rare event, its occurrence is observed in large finite populations. Extinction occurs when fluctuations owing to random transitions act as an effective force that drives one or more components or species to vanish. Although there are many random paths to an extinct state, there is an optimal path that maximizes the probability to extinction. In this paper, we show that the optimal path is associated with the dynamical systems idea of having maximum sensitive dependence to initial conditions. Using the equivalence between the sensitive dependence and the path to extinction, we show that the dynamical systems picture of extinction evolves naturally towards the optimal path in several stochastic models of epidemics.
Determining the conditions for epidemic extinction is an important public health problem. Global eradication of an infectious disease has rarely been achieved, but it continues to be a public health goal for polio  and many other diseases, including childhood diseases. More commonly, disease extinction, or fade out, is local and may be followed by a reintroduction of the disease from other regions [2,3]. Extinction may also occur for individual strains of a multi-strain disease , such as influenza or dengue fever. Since extinction occurs in finite populations, it depends critically on local community size [5–8]. Moreover, it is important to know how parameters affect the chance of extinction for predicting the dynamics of outbreaks and for developing control strategies to promote epidemic extinction . The determination of extinction risk is also of interest in related fields, such as evolution and ecology. For example, in the neutral theory of ecology, bio-diversity arises from the interplay between the introduction and extinction of species [10,11].
In general, extinction occurs in discrete, finite populations undergoing stochastic effects owing to random transitions or perturbations. The origins of stochasticity may be internal to the system or may arise from the external environment. Small population size, low contact frequency for frequency-dependent transmission, competition for resources and evolutionary pressure , as well as heterogeneity in populations and transmission , may all be determining factors for extinction to occur.
Extinction risk is affected by the nature and strength of the noise , as well as other factors, including outbreak amplitude  and seasonal phase occurrence . For large populations, the intensity of internal population noise is generally small. However, a rare, large fluctuation can occur with non-zero probability and the system may be able to reach the extinct state. For respiratory diseases such as severe acute respiratory syndrome (SARS), super-spreading may account for these large stochastic fluctuations . Since the extinct state is absorbing owing to effective stochastic forces, eventual extinction is guaranteed when there is no source of reintroduction [18–20]. However, because fade outs are usually rare events in large populations, typical timescales for extinction may be extremely long.
Stochastic population models of finite populations, which include extinction processes, are effectively described using the master equation formalism. Stochastic master equations are commonly used in statistical physics when dealing with chemical reaction processes  and predict probabilities of rare events occurring at certain times. For many problems involving extinction in large populations, if the probability distribution of the population is approximately stationary, the probability of extinction is a function that decreases exponentially with increasing population size. The exponent in this function scales as a deterministic quantity called the action . It can be shown that a trajectory that brings the system to extinction is very likely to lie along a most probable path called the optimal extinction trajectory or optimal path. It is a remarkable property that a deterministic quantity such as the action can predict the probability of extinction, which is inherently a stochastic process [9,23].
Locating the optimal path is desirable because the quantity of interest, the extinction rate, depends on the probability to traverse this path, and the effect of a control strategy on extinction rate can be determined by its effect on the optimal path . By employing an optimal path formalism, we convert the stochastic problem to a mechanistic dynamical systems problem. In contrast to approaches based on diffusive processes, which are valid only in the limit of large system sizes [24,25], this dynamical systems approach can give accurate estimates for the extinction time even for small populations if the action is sufficiently large. Additionally, unlike other methods that are used to estimate lifetimes, this approach enables one both to estimate lifetimes and to draw conclusions about the path taken to extinction. This more detailed understanding of how extinction occurs may lead to new stochastic control strategies .
In this paper, we show that locating the optimal extinction trajectory can be achieved naturally by evolving a dynamical system that converges to the optimal path. The method is based on computing finite-time Lyapunov exponents (FTLE), which have previously been used to find coherent structures in fluid flows [26–31]. The FTLE provides a measure of how sensitively the system's future behaviour depends on its current state. We argue that the system displays maximum sensitivity near the optimal extinction trajectory, which enables us to dynamically evolve towards the optimal escape trajectory using FTLE calculations. For several models of epidemics that contain internal or external noise sources, we illustrate the power of our method to locate optimal extinction trajectories. Although our examples are taken from infectious disease models, the approach is general and is applicable to any extinction process or escape process.
2. Stochastic modelling in the large population limit
To introduce our idea of dynamically constructing an optimal path to extinction in stochastic systems, we show its application to a stochastic susceptible–infectious–recovered (SIR) epidemiological model. Details of the SIR model can be found in the electronic supplementary material, §1. Figure 1 shows the probability density of extinction prehistory in the susceptible–infectious (SI) plane. The probability density was numerically computed using 20 000 stochastic SIR trajectories that ended in extinction. Trajectories are aligned at their extinction point. From the extinction point, the prehistory of each trajectory up to the last outbreak of infection is considered. Small fluctuations in the infectious population are not considered in identifying the last outbreak. In this way, we restrict the analysis to the interval between the last large outbreak of infection and the extinction point. The resulting (S,I) pairs of susceptible and infectious individuals are then binned and plotted in the SI plane . The resulting discrete density has been colour coded so that the brighter regions correspond to higher density of trajectories. The figure shows that, among all the paths that the stochastic system can take to reach the extinct state, there is one path that has the highest probability of occurring. This is the optimal path to extinction. One can see that the optimal path to extinction lies on the peak of the probability density of the extinction prehistory. It should be noted that extinction for the stochastic SIR model has been studied previously .
The optimal path can be obtained using methods of statistical physics. In figure 1, the numerical prediction of the entire optimal trajectory for the stochastic SIR system has been overlaid on the probability density of extinction prehistory that was found using stochastic simulation. The trajectory spirals away from the endemic state, with larger and larger oscillations until it hits the extinct state. The agreement between the stochastically simulated optimal path to extinction and the predicted optimal path is excellent.
The curve of figure 1 is obtained by recasting the stochastic problem in a deterministic form. The evolution of the probability of finding a stochastic system in a given state X at a given time t is described by the master equation . Solving the master equation would provide a complete description of the time evolution of the stochastic system, but in general it is a difficult task to obtain explicit solutions for the master equation. Thus, one generally resorts to approximations to the solution; i.e. one considers an ansatz for the probability density. In this case, since extinction of finite populations is a rare event, we will be interested in examining events that occur in the tail of the probability distribution. Therefore, the distribution is assumed to take the form: 2.1where ρ(X, t) is the probability density of finding the system in the state X at time t, N is the size of the population, x = X/N is the normalized state (e.g. in an epidemic model, the fraction of the population in the various compartments), and 𝒮 is a deterministic state function known in classical physics as the action. Equation (2.1) describes the relationship between the action and the probability density and is based on an assumption for how the probability scales with the population size. The action is the negative of the natural log of the stationary probability distribution divided by the population size. Therefore, the probability (if we normalize the population) is roughly given by the exponential of the action. Intuitively, equation (2.1) expresses the assumption that the probability of occurrence of extreme events, such as extinction, lies in the tails of the probability distribution, which falls steeply away from the steady state.
This approximation leads to a conserved quantity that is called the Hamiltonian . From the Hamiltonian, one can find a set of conservative ordinary differential equations (ODEs) that are known as Hamilton's equations. These ODEs describe the time evolution of the system in terms of state variables x, which are associated with the population in each compartment. For the SIR example, x is the vector 〈S,I,R〉. In addition to the state variables, the equations contain conjugate momenta variables, px. The conjugate momenta, or noise, account for the uncertainty associated with the system being in a given state at a given time owing to the stochastic interactions among the individuals of the population. These ODEs can be constructed from information in the master equation about the possible transitions and transition rates in the system. Details can be found in appendix A.
Integration of the ODEs with the appropriate boundary conditions will then give the optimal evolution of the system under the influence of the noise. Boundary conditions are chosen to be fixed points of the system. A typical case is shown schematically in figure 2a. Deterministically, the endemic state is attracting and the extinct state repelling. However, introducing stochasticity allows the system to leave the deterministic manifold along an unstable direction of the endemic state, corresponding to non-zero noise. Stochasticity leads to an additional extinct state that arises owing to the general non-Gaussian nature of the noise. For the extinction process of figure 1, boundary conditions were the system leaving the endemic steady state and asymptotically approaching the stochastic extinct state.
In general, the optimal extinction path is an unstable dynamical object, and this reflects extinction as a rare event. This has led many authors to consider how extinction rates scale with respect to a parameter close to a bifurcation point [9,33,36,37], where the dynamics are very slow. For an epidemic model, this means that the reproductive rate R0 should be greater than but very close to 1. However, most real diseases have R0 larger than 1.5, which translates into a faster growth rate from the extinct state. In general, in order to obtain analytical scaling results, one must obtain the ODEs for the optimal path either analytically (using the classical theory of large fluctuations mentioned within this section and described in detail in appendix A) or numerically (using shooting methods for boundary value problems). This task may be impossible or extremely cumbersome, especially when the system is far from the bifurcation point. In the following section, we demonstrate how to evolve naturally to the optimal path to extinction using a dynamical systems approach.
3. Finite-time Lyapunov exponents
We consider a velocity field that is defined over a finite time interval and is given by Hamilton's equations of motion. Such a velocity field, when starting from an initial condition, has a unique solution. The continuous dynamical system has quantities, known as Lyapunov exponents, which are associated with the trajectory of the system in an infinite time limit, and which measure the average growth rates of the linearized dynamics about the trajectory. To find the FTLE, one computes the Lyapunov exponents on a restricted finite time interval. For each initial condition, the exponents provide a measure of their sensitivity to small perturbations. Therefore, the FTLE is a measure of the local sensitivity to initial data. Details regarding the FTLE can be found in the electronic supplementary material, §2.
The FTLE measurements can be shown to exhibit ‘ridges’ of local maxima. The ridges of the field indicate the location of attracting (backward time FTLE field) and repelling (forward time FTLE field) structures. In two-dimensional space, the ridge is a curve, which locally maximizes the FTLE field so that transverse to the ridge one finds the FTLE to be a local maximum. What is remarkable is that the FTLE ridges correspond to the optimal path trajectories, which we heuristically argue in the electronic supplementary material, §3. The basic idea is that since the optimal path is inherently unstable, the FTLE shows that, locally, the path is also the most sensitive to initial data. Figure 2b shows a schematic that demonstrates why the optimal path has a local maximum to sensitivity. If one chooses an initial point on either side of the path near the endemic state, the two trajectories will separate exponentially in time. This is due to the fact that both extinct and endemic states are unstable, and the connecting trajectory defining the path is unstable as well. Any initial points starting near the optimal path will leave the neighbourhood in short time.
4. Evolving towards the optimal path using finite-time Lyapunov exponents
We now apply our theory of dynamical sensitivity to the problem of locating optimal paths to extinction for several examples. We consider the case of internal fluctuations, where the noise is not known a priori, as well as the case of external noise. In each case, the interaction of the noise and state of the system begins by finding the equations of motion that describe the unstable flow. These equations of motion are then used to compute the ridges corresponding to maximum FTLE, which in turn correspond to the optimal extinction paths .
4.1. Extinction in a branching–annihilation process
For an example of a system with internal fluctuations, which has an analytical solution, consider extinction in the stochastic branching–annihilation process 4.1where λ and μ > 0 are constant reaction rates [39,40]. Equation (4.1) is a single species birth–death process and can be thought of as a simplified form of the Verhulst logistic model for population growth . The mean field equation for the average number of individuals n in the infinite population limit is given by . The stochastic process given by equation (4.1) contains intrinsic noise that arises from the randomness of the reactions and the fact that the population consists of discrete individuals. This intrinsic noise can generate a rare sequence of events that causes the system to evolve to the extinct state. The probability Pn(t) to observe, at time t, n individuals is governed by the master equation 4.2
The Hamiltonian associated with this system is 4.3where q is a conjugate coordinate related to n through a transformation , and p plays the role of the momentum. The equations of motion are given by 4.4aand 4.4bThe mean field is retrieved in equation (4.4a) when p = 0 (no fluctuations or noise). The Hamiltonian has three zero-energy curves. The first is the mean-field zero-energy line p = 0 (no fluctuations), which contains two unstable points h1 = (p,q) = (0,λ/μ) and h0 = (p,q) = (0,0). The second is the extinction line q = 0 (trivial solution), which contains another unstable point h2 = (p,q) = (−1,0). The third zero-energy curve is non-trivial and is 4.5
The segment of the curve given by equation (4.5), which lies between −1 ≤ p ≤ 0 corresponds to a (heteroclinic) trajectory that exits, at t = −∞, the point h1 along its unstable manifold and enters, at t = ∞, the point h2 along its stable manifold. This trajectory is the optimal path to extinction and describes the most probable sequence of events, which evolves the system from a quasi-stationary state to extinction .
To show that the FTLE evolves to the optimal path, we calculate the FTLE field using the system of Hamilton's equations given by equations (4.4a,b). Figure 3a shows both the forward and backward FTLE plot computed for the finite time T = 6, with λ = 2.0 and μ = 0.5. In this example, as well as the following two examples, T was chosen to be sufficiently large so that one obtains a measurable exponential separation of trajectories. In figure 3a, one can see that the optimal path to extinction is given by the ridge associated with the maximum FTLE. In fact, by overlaying the forward and backward FTLE fields, one can see all three zero-energy curves including the optimal path to extinction. Also shown in figure 3a are the analytical solutions to the three zero-energy curves given by p = 0, q = 0 and equation (4.5). There is excellent agreement between the analytical solutions of all three curves and the ridges that are found through numerical computation of the FTLE flow fields.
It is possible to compute analytically the action along the optimal path for a range of λ/μ values. Using equation (4.5), it is easy to show that the action 𝒮 is 4.6
It is clear from equation (4.6) that the action scales linearly with λ/μ.
4.2. Susceptible–infectious–susceptible epidemic model: external fluctuations
We now consider the well-known problem of extinction in a susceptible–infectious–susceptible (SIS) epidemiological model, which is a core model for the basis of many recurrent epidemics. The SIS model is described by the following system of equations: 4.7aand 4.7bwhere μ denotes a constant birth and death rate, β represents the contact rate and γ denotes the recovery rate. Assuming the total population size is constant and can be normalized to S + I = 1, then equations (4.7a,b) can be rewritten as the following one-dimensional process for the fraction of infectious individuals in the population: 4.8
The stochastic version of equation (4.8) is given as 4.9where η(t) is uncorrelated Gaussian noise with zero mean and σ is the standard deviation of the noise intensity. The noise models random migration to and from another population [15,36].
Equation (4.9) has two equilibrium points given by I = 0 (disease-free state) and I = 1 − (μ + γ)/β (endemic state). Using the Euler–Lagrange equation of motion  from the Lagrangian determined by equation (4.9) () along with the boundary conditions given by the extinct and endemic states, one finds that the optimal path to extinction (as well as its counterpart path from the disease-free state to the endemic state) is given by .
As in the first example, one can numerically compute the optimal path to extinction using the FTLE. Figure 3b shows the forward and backward FTLE plot computed for T = 5, with β = 5.0 and κ = μ + γ = 1.0. Note that we can consider the combination μ + γ since the Lagrangian depends only on the combination rather than on either parameter separately. In figure 3b, one can see that the optimal path from the endemic state to the disease-free state is given by the ridge associated with the local maximum FTLE. Also shown in figure 3b is the counterpart optimal path from the disease-free state to the endemic state (found by computing the backward FTLE field). In addition, the agreement with the analytical prediction is excellent, as shown in figure 3b.
If one solves equation (4.9) for I(t) and substitutes both I(t) and into the Lagrangian, then one can analytically find an expression for the action along the optimal path. The expression for the action is a function of κ and the reproductive number R0 and is given by 4.10where R0 = β/κ.
4.3. Susceptible–infectious–susceptible epidemic model: internal fluctuations
We next consider the one-dimensional stochastic version of the SIS epidemic model for a finite population with internal fluctuations using the transition rates found in the electronic supplementary material, §4. Using the formalism of Gang , one then has the following Hamiltonian associated with this model: 4.11where I is the fraction of infectious individuals and p is the momentum. Internal fluctuations arise from the random interactions of the population. Although there is no analytical solution for the optimal path to extinction, we can once again determine the optimal path by computing the FTLE flow field associated with this system. Figure 4 shows the forward FTLE plot computed using Hamilton's equation of motion for T = 10, with β = 2.0 and κ = μ + γ = 1.0, and as in previous examples, the optimal path to extinction from the endemic state to the disease-free state is apparent. Note that the non-zero momentum corresponding to the extinct state qualitatively agrees with similar boundary conditions found in Dykman et al.  and is associated with non-zero probability flux.
Once again, it is possible to compute the action along the optimal path for a range of values of the reproductive number R0. In contrast to the prior two examples, here the action must be computed numerically. Moreover, even the optimal path must be found numerically using the FTLE plot generated for each value of R0.
Given an R0, we computed the forward FTLE flow field using a grid resolution of 0.005 in both position and momentum space. To find the optimal path corresponding to the ridge of maximal FTLE values, we let the deterministic, endemic steady state be the starting point z0 of the path. Then a second point z1 on the path was found by locating the maximum of the FTLE values in an arc of radius 20 grid points spanning nearly π radians for negative momentum values and originating at z0. Subsequent points zn+1 were found by taking the maximum FTLE value on an arc of radius 10 grid points spanning nearly π radians originating at the most recent point zn and centred around the vector zn − zn−1 until the extinction line was reached. The complete optimal path was estimated from the zn using cubic spline interpolation. Once the optimal path had been found, the action was computed by numerically integrating the Lagrangian along this path.
Figure 5a shows the numerically computed action versus reproductive number over the range 1.1 ≤ R0 ≤ 20. The inset of figure 5a shows the portion of figure 5a from 1.1 ≤ R0 ≤ 2. Also shown in the inset of figure 5a is an analytical, asymptotic scaling result that is valid for values of R0 close to unity. As can be see in figure 5a, there is a good agreement between the numerically computed action and the analytically computed action near R0 = 1. Details of the derivation of the analytical scaling law can be found in the electronic supplementary material, §5.
Figure 5b shows the numerically simulated mean extinction time versus reproductive number for the one-dimensional stochastic SIS model for a finite population (see §4.3). Also shown in figure 5b is the analytical prediction found using the previously mentioned scaling law derived for values of R0 close to unity. As one can see, the agreement is excellent.
In all of the examples of §4, we have shown that the optimal path to extinction is equated with having a (locally) maximal sensitivity to initial conditions. Even though many possible paths to extinction exist, the dynamical systems approach converges to the path that maximizes extinction. The parameter values chosen for the three examples are such that the extinct and endemic states are far away from each other. Therefore, in general, there will be no possible approximate analytical treatment as performed in Dykman et al. . In addition, we have shown how to constructively compute numerically the action for a wide range of reproductive numbers. Our method allows for the computation of extinction times and can be extended to high-dimensional problems.
Because the method is general, and unifies dynamical systems theory with the probability of extinction, we expect that any system found in other fields can be understood using this approach. Indeed, in problems of general extinction, it is now possible to evolve naturally to the optimal path, and thus predict the path that maximizes the probability to extinction. Future work in this area will include improved sampling methods to find the optimal path more efficiently in higher dimensional models. Specific applications of optimal path location in the future will include spatio-temporal extinction processes such as those that occur in pre-vaccine measles [7,43] and multi-strain extinction in diseases such as dengue .
Finally, the optimal path method may lead to novel monitoring and control strategies. In one biological application, knowledge of the most probable path to extinction may help with monitoring an environment and with providing an estimate of the likelihood of extinction based on where the population lies relative to the path. It is known that once a trajectory has left the neighbourhood of the endemic state, most paths to extinction occur near the optimal path. This phenomenon can be seen in figure 1, where the optimal path lies on the peak of the probability density of extinction prehistory. Therefore, the optimal path provides a good location to monitor the system for possible extinction events. Furthermore, although the momenta (noise effects) are not directly observable, it may be possible to infer them dynamically  from time series data of observable quantities in conjunction with the equations for the time evolution of position and momentum variables. Knowledge of where a system lies in position-momentum space could provide an estimate for how quickly a desired epidemic extinction could occur or could provide the risk of extinction for a species one wishes to conserve.
In yet another application, knowledge of the optimal path to extinction has the advantage that, in the presence of noise (that is estimated from data) and a known population of infectious individuals, it may be possible to develop better vaccine controls that reduce the time to extinction. Figure 2a shows a schematic of the path to extinction, where the extinct state is a saddle with stable and unstable directions. For many epidemic models, the extinct state can be shown to have a similar geometry of stable and unstable directions. An approach to the extinct state on the optimal path will lead to the fastest time to extinction. Moreover, since the extinct state has a saddle structure in the presence of noise, it may be controlled with projection methods [46,47] or probabilistic techniques .
One can consider instituting a method of parameter control, where the parameters could be vaccine levels or treatment of infectious individuals. Combined with the monitoring techniques mentioned above, the control method will allow one to move an existing state that deviates from the optimal path closer to the optimal path as time evolves. By adjusting the parameters, we may target the stable directions of the path when we are close to epidemic extinction . Comparing observations with model predictions of the optimal path allow us to use controls to minimize the time to extinction.
The authors gratefully acknowledge support from the Office of Naval Research, the Air Force Office of Scientific Research, the Army Research Office and the National Institutes of Health. I.B.S and L.B.S. are supported by Award Number R01GM090204 from the National Institute Of General Medical Sciences. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institute of General Medical Sciences or the National Institutes of Health. We also gratefully acknowledge M. Dykman for helpful discussions.
Appendix A. Theory of large fluctuations
Letting N denote the population size, the state variables X ∈ ℝn of the system describe the components of a population, while the random state transitions, which govern the dynamics are described by the transition rates W(X,r), where r ∈ ℝn is an increment in the change of X. In the large population limit without any fluctuations, the mean field equations are given by the system .
Consideration of the net change in increments of the state of the system results in the following master equation for the probability density ρ(X,t) of finding the system in state X at time t: A 1
The solution to the master equation (A 1) has a peak around a stable steady state in the limit of large population (N ≫ 1) with width ∝ N1/2 [9,23]. However, since we are interested in the probability of extinction, we will consider the tail of the distribution, which gives the probability of having a small number of individuals. The tail of the distribution can be obtained by employing the ansatz given by equation (2.1), which is an assumption for how the probability density scales with population size [9,22,49]. Equation (2.1) also implies that a maximum in the extinction probability can be found minimizing the action over a set of extinction paths starting from the stable endemic state. The assumption of this functional form for ρ allows the action to be derived from properties of the master equation.
The density, ρ(X,t), can be found by substituting the ansatz given by equation (2.1) into the master equation. The resulting equations for the action are given by the Hamilton–Jacobi equation for a Hamiltonian, H, given by ∂𝒮/∂t + H(x,∂𝒮/∂x;t) = 0, where we have normalized the population components and rates, respectively, as x = X/N, w(x,r) = W(X,r)/N.
Following the classical mechanics convention, define a conjugate momentum to the state space, x, by letting p = ∂𝒮/∂x and where H(x,p;t) is the classical Hamiltonian . The Hamiltonian, H, depends both on the state of the system, x, as well as the momentum, p, which provides an effective force owing to stochastic fluctuations on the system. The Hamiltonian equations of motion provide trajectories in time of x(t) and p(t), and as such, describe a set of paths going from an initial point at time ti to some final point at time tf in (x,p) space. For a given path, the action is given by , and as such determines the exponent of the probability of observing that path. (It should be noted that instead of using the Hamiltonian representation, one could use the Lagrangian representation , which results in an equivalent solution.)
The action reveals much information about the probability evolution of the system, from scaling near bifurcation points in non-Gaussian processes to rates of extinction as a function of epidemiological parameters [9,49]. As already stated, in order to maximize the probability of extinction, one must minimize the action 𝒮. The minimizing formulation entails finding the solution to the Hamilton–Jacobi equation, which means that one must solve the 2n-dimensional system of Hamilton's equations [, ] for x and p, where the Hamiltonian is explicitly given as A 2
The appropriate boundary conditions of the system are such that a solution starts at a non-zero state, such as an endemic state, and asymptotically approaches one or more zero components of the state vector, representing a disease-free state. Therefore, a trajectory that is a solution to the two-point boundary value problem determines a path, which in turn yields the probability of going from the initial state to the final state. The optimal path to extinction is the path that minimizes the action in either the Hamiltonian or Lagrangian representation.
We compute the trajectory satisfying the Hamiltonian system that has as its asymptotic limits in time the endemic state as t → −∞ and the extinct state as t →+∞. The momentum p represents the force of the fluctuations on the population, and this momentum changes the stability of the equilibrium points. Both the endemic and extinct states have attracting and repelling directions for p ≠ 0, as shown schematically in figure 2.
Given optimal path trajectories (x(t), p(t)), the action with correct limits is found by .
- Received March 22, 2011.
- Accepted April 21, 2011.
- This journal is © 2011 The Royal Society