Hamiltonian dynamics of the SIS epidemic model with stochastic fluctuations

Empirical records of epidemics reveal that fluctuations are important factors for the spread and prevalence of infectious diseases. The exact manner in which fluctuations affect spreading dynamics remains poorly known. Recent analytical and numerical studies have demonstrated that improved differential equations for mean and variance of infected individuals reproduce certain regimes of the SIS epidemic model. Here, we show they form a dynamical system that follows Hamilton’s equations, which allow us to understand the role of fluctuations and their effects on epidemics. Our findings show the Hamiltonian is a constant of motion for large population sizes. For small populations, finite size effects break the temporal symmetry and induce a power-law decay of the Hamiltonian near the outbreak onset, with a parameter-free exponent. Away from the onset, the Hamiltonian decays exponentially according to a constant relaxation time, which we propose as a metric when fluctuations cannot be neglected.

spatial structures in epidemic models and pathogen variability have been under investigation for some time 20,21 , and they have been linked to chaotic dynamics in local population 22 . Experiments on the effects of migration between metapopulations, i.e. similar populations but spatially separated, subjected to temporal fluctuations have shown that pathogen prevalence is greatly influenced by the nature of the fluctuation 23 , highlighting the interplay between synchronization and pathogen prevalence in epidemics 24 .
Traditionally, the detailed examination of fluctuations-either temporal or spatial-and their effects on system dynamics have been largely described by correlation functions 23,25 . More recently, autocorrelation functions have been used to reveal the nature and general aspects of fluctuations in a simple agent-based epidemic model for a population of size N, in which temporal fluctuations are divided into two broad classes: gaussian and non-gaussian 26 . In the gaussian regime, the prevalence of the disease is well described by its instantaneous average, finite variance, and higher cumulants can be neglected. This is remarkable as it allows one to derive the exact contributions of fluctuations to disease outbreaks in the asymptotic limit N  1. Here, we show that the dynamical equations form a Hamiltonian dynamical system, and the way external noise can be incorporated to model disease outbreaks. This approach allows us to discuss quantitatively the relevant scales of the problem, and interpret the resulting Lagrangian and canonical transformations.

Model
We begin our discussion using the susceptible-infectious-susceptible (SIS) epidemic model. The SIS model describes the dissemination of a single communicable disease in a susceptible population of size N. The transmission of the pathogen occurs when infectious hosts transmit the disease pathogen to healthy susceptible individuals. The infectious period extends throughout the whole course of the disease until the recovery of the patient, warranting a two-stage model: either infected or susceptible. The essence of the model is summarized by inset in Fig. 1.
The traditional formulation of the problem assumes the random-mixing hypothesis (see Introduction) holds for a large population size N1, compromised of statistically equivalent individuals. Under these circumstances, the only relevant variable is the instantaneous density of infected elements ρ(t), which means that fluctuations can be safely neglected. Furthermore, ρ(t) decreases with rate γρ, where γ is the recovery rate. New infections per unit of time (disease incidence) are proportional to αρ(1 − ρ), i.e., they depend on the chance that infected elements interact with susceptible ones, with intensity given by the transmission rate α. This picture provides an interpretation where ρ(t) is continuously exchanged between two compartments, leading a simple description called compartmental equation: dρ(t)/dt = αρ(1 − ρ) − γρ. For the sake of convenience, redefine the timescale as τ ≡ αt (1) fails to reproduce the data (dashed line). The forward-derivative dσ 2 /dτ from data (circles) also agrees with the formula in Eq. (5b) (line). All the lines are drawn using the simulated data for 〈ρ(τ)〉, σ 2 (τ), and Δ 3 (τ).
Clearly, the equilibrium density can either be ρ eq = 0 or ρ eq = ρ 0 . Also, ρ 0 is related to the basic reproduction number R 0 = N(α/γ) which provides an estimate on the number of new infections per generation 27 .
In light of its long age, compartmental equations have met considerable success in predicting the time evolution of disease outbreaks, providing valuable insights for intervention strategies and funding allocation 28 . However, outbreaks that fail to meet the underlying hypotheses (random mixing and large population of statistically equivalent elements) can contradict compartmental equations. These inconsistencies are largely attributed to stochastic effects and their inherent fluctuations 2 .

Improved Compartmental Equations
Stochastic variables are known to cause the emergence of critical phenomena in computer simulations of epidemic models, under certain parameter ranges 29,30 . One key ingredient common to critical phenomena is the scale invariance of fluctuations. This special symmetry remains the foundation of cooperative phenomena and critical phase transitions, whose contributions span over a broad set of research fields such as condensed-matter, quantum field theories, and neuroscience to name a few [31][32][33][34][35][36] . In these systems, fluctuations occur in all sizes and dictate the general behavior of the problem, which prompts for in-depth studies of their effects in epidemic models.
Stochastic epidemic models include non-deterministic events that intrinsically occurs during the course of the disease spreading process. Examples of these events include the transmission of the pathogen, and the elapsed time required for the complete recovery of patients. The inclusion of these effects brings the models closer to more realistic expectations, which in turn can deviate from predictions using compartmental equations [37][38][39] . One prime example is the effects of the absorbing state in the SIS model, which has been examined in detail by Nåsell and others [40][41][42] . New experiments on this subject provide evidence that temporal fluctuations can drastically alter the prevalence of pathogens 23 . Spatial heterogeneity also introduces an extra layer of complexity as it may trap or delay the pathogen transmission 20 . As a result, the requirements of statistical equivalence may not hold for all scales. To deal with this issue, stochastic formulations and numerical simulations have been the default tools to investigate fluctuations in disease outbreaks. In what follows, we outline the general ideas behind the agent-based approach to describe the SIS model. A more detailed account of the results listed here are explained in detail in refs 26,43,44 .
Our discussion assumes the disease spreading follows a Markov chain in discrete time δt. Moreover, δt is such that at most a single recovery or transmission event is likely to occur during the course of its duration. Under these requirements, the master equation of the SIS model in discrete time reads Here, P μ (t) refers to the instantaneous probability to observe the system in the μ-th configuration. Configuration labels follow the binary ruling where A k is the adjacency matrix, n k represents the k-th occupation operator (with eigenvalues n k = 1 if infected, 0 otherwise), and σ + k are the localized spin-1/2 ladder operators that produce the transition S → I. Clearly, k σ − produce the opposite transitions, I → S in relation the k-th agent. As notation, the hat symbol always accompanies operators to quickly distinguish them from numbers.
The master equation Eq.
(2) provides the means to evaluate the time evolution of relevant statistical moments of ρ(t). Notice that the average density of infected agents in the system readŝ k k , which recovers the random mixing hypothesis. In those particular instances, the complete time evolution of the system comprehends a set of hierarchical equations that involves the statistical moments of ρ(t), as shown in refs. 43,45 . More explicitly 45 , the first two equations for instantaneous mean 〈ρ〉 and variance www.nature.com/scientificreports www.nature.com/scientificreports/ where Δ 3 (τ) = 〈ρ 3 (τ)〉 − 〈ρ(τ)〉 3 . These results find excellent agreement with simulated data using an ensemble with 10 6 replicas starting from the same initial condition (see Fig. 1).
Comparing Eqs (1) and (5), the case that considers temporal fluctuations decays faster than the compartmental equation by σ 2 (τ), even in the regime  N 1. Both equations are equivalent whenever σ(τ) becomes irrelevant compared to 〈ρ〉. Therefore, a generalization of compartmental equations for the SIS model is readily available by retaining both mean and variance, neglecting higher statistical moments. Thus, the dynamical system describes a gaussian variable evolving along time. The skewness coefficient vanishes as a direct consequence of this assumption, so that Δ 3 (τ) ≈ 3〈ρ(τ)〉σ 2 (τ). For  N 1, the resulting equations are We emphasize that the variance in Eq. (6a) slows down the growth rate of 〈ρ(τ)〉, recalling the Allee effect 28,46 .
The dynamical system represented by Eq. (6) can also be obtained from a stochastic expansion following the guidelines in ref. 47 . In this case, assume the density ρ = 〈ρ〉 + η is well described by its instantaneous average plus some noise function η(τ). We assume 〈η〉 = 0 and 〈η 2 〉 = σ 2 (τ) for the sake of consistency. The additional requirement η ρ | | 〈 〉  ensures stochastic effects act as perturbations. Expanding Eq. (1) as a Taylor series around ρ = 〈ρ〉, and then taking the ensemble average, produces Taking the ensemble average, one obtains (1/2)(d/dt)lnσ 2 = [ρ 0 − 2〈ρ〉] − 〈η 3 〉/σ 2 . Now, since η mimics a gaussian variable, 〈η 3 〉 = 0 and Eq. (6b) is recovered as well. Equation (6) can be further combined into a single second-order differential equation 26  The constants c 1 and c 2 depend solely on the initial conditions. The special case = c c 2 1 2 recovers the usual solution of Eq. (1). We assumed that fluctuations behave as gaussian fluctuations. While reasonable for various situations, the assumption does not hold for γ/α around unity or small population sizes, according to numerical simulations 26 , in which Eq. (5b) should be used instead of Eq. (6b). For the sake of completeness, there is another solution in which σ 2 0 .

Hamilton Equations
The fact that the dynamical system Eq. (6) can be combined into a single second-order differential equation suggests an interpretation of the epidemic model in terms of Hamilton equations 48 . Hamiltonian systems are ubiquitous in Physics, serving as basis to describe and explain countless physical phenomena. The hallmark of systems are the Hamilton equations: where q(t) and p(t) are conjugated variables, and the Hamiltonian function  encodes some information about the problem-usually associated with energy for conservative systems but not restricted to them. Besides classical mechanics and related areas, quantum field theories and statistical mechanics are deeply intertwined with Hamilton's principle and Liouville theorem. Despite its usefulness in Physics, Hamilton formulation and surrounding principles are rarely used in population dynamics, ecological problems, or epidemic models, where first-order differential equations are dominant. The lack of second-order differential equations in these areas, although not prohibitive, raises questions about the description of the dynamics, as discussed extensively in ref. 49  www.nature.com/scientificreports www.nature.com/scientificreports/ because it means some interactions and stochastic effects acting on the system remain unaccounted. Significant advances in the Hamiltonian formulation of stochastic epidemic models have been obtained using the eikonal approximation, with emphasis on the disease extinction and vaccination 50,51 . Generalizations for heterogeneous networks have also been investigated, providing improved control strategies 52,53 . Even more, the aforementioned Hamiltonian formulation explains other effects, such as the emergence of noise-induced metastable states 54 . The crucial ingredient of these studies concerns the eikonal approximation, in which either the generating function or the density of infected is written as the exponential of the classical action of the system 55 . For populations with fixed size, one variable describes 〈ρ(t)〉, while its conjugate variable encompasses fluctuations 52 . However, the description of epidemics using this conjugate pair remains a complex task because the conjugate variable is not directly related to familiar statistics such as the standard deviation or variance. In what follows, instead of using the eikonal approach, we argue that Hamilton dynamics for the stochastic SIS model can be constructed from 〈ρ(τ)〉 and σ(τ).
In view of the inherent stochasticity behind disease spreading, it seems necessary to determine whether Eq. (6) form a Hamiltonian system or not. A brief inspection shows the pair (〈ρ〉, σ 2 ) does not satisfy the usual Hamilton equations. The solution to this issue is obtained by assuming, instead, that the correct conjugated pair is (〈ρ〉, h(σ 2 )), where h(x) is some analytical function. Inspiration from common pairs of conjugate variables can be used to refine the choice of h(x). For instance, the product 〈ρ〉 × h(σ 2 ) should be dimensionless, in close analogy the scalar product between position and wave vectors. One possible candidate is h(x) = x −1/2 , which entails 1/σ as the conjugated variable to 〈ρ〉.
Define the dynamical variables q(τ) = 〈ρ(τ)〉 and p(τ) = 1/σ(τ) to describe the SIS model. In addition, consider the following Hamiltonian Plugging these expressions in Eq. (10), one obtains the equations of motion: Thus, at first glance  appears to be a valid candidate to describe the SIS model. Even more, replacing (q, p) by Eq. (9) in Eq. (11) shows that the Hamiltonian is a constant of motion 1/2 . The upper index in  ∞ is a reminder that calculations take place in the absence of finite size corrections.
However, taking finite size corrections into account changes drastically the notion of  as a constant of motion. In fact, as Fig. 2 depicts, numerical simulations for finite populations reveal  changes continuously along time until equilibrium sets in, akin to a non-conservative system. The precise meaning of  in the epidemiological context is still murky, at best. A detailed analysis of correlations between changes in  and the spreading pattern of real outbreaks is mandatory to understand the action-reaction analogy. In the meantime, it is instructive to study  for τ 1 and τ 1 (see Fig. 2). For τ 1, where incidentally fluctuations vary the most (see Fig. 1), a remarkable feature appears via the relation ~τ λ − with λ = 1/2. In particular, the exponent λ seems insensitive to changes in the epidemiological parameter γ. This parameter-free behavior is not observed for the remaining statistics, 〈ρ(τ)〉 and σ(τ). Power-laws are crucial to identify scaling relations and the emergence of universal features, and they are usually related to the symmetry of the problem rather than microscopic details. Here, evidence of universal behavior is captured by the data collapse ρ / 0 2  (Fig. 2a inset). From these observations, we can infer fluctuations play a larger role in the early disease spreading, being largely independent of exact values of epidemiological parameters.
An effective decay τ τ − e / eff describes the general behavior of  away from the outbreak onset. The relaxation time τ eff depends on N and the ratio γ/α, and it can be estimated from data by fitting  to an exponential function plus a constant. Alternatively, it can be evaluated as eff 0 From a formal point of view, the evaluation of τ eff requires the solutions of Eqs (5a) and (5b) in Eq. (11), followed by an integration. Surely, the procedure is arguably more demanding than estimating R 0 . However, as others have reasoned before, R 0 provides a naive estimation on secondary infections because the growth rate of the outbreak changes continuously along time 56 . In contrast, τ eff mimics a constant of motion.

Lagrangian and Canonical Transformations
Another insight from τ eff links the temporal integral of  with the mechanical action S. A formal connection with S is desirable because it brings a large machinery revolving around variational principles and conservation laws. However, the action   is a functional of the Lagrangian . It turns out that  can obtained from  by inspection. From Eqs (15) and (7) 0 where we have used Eq. (6a) and considered only the positive root. Thus,  is proportional to the standard deviation while the action entails the accumulated deviation over the course of the outbreak. To check our result for large populations N 1  , the minimal action recovers Eq. (1) as expected for a noise-free system. In general, the equation of motion reads (τ) to the variance of the system σ 2 (τ): This picture is consistent with addition of a noise function σ ext 2 (τ) to Eq. (6a). The perturbed Lagrangian ′  describes, ultimately, the time evolution of the disease prevalence in environments with noise. Note that this description differs from the usual derivation of Langevin equations, in which the noise function (force) r(τ) couples linearly with q, i.e.,   τ τ ′ = − r q ( ) ( . The covariance matrix can estimated or modeled directly from data, pro- www.nature.com/scientificreports www.nature.com/scientificreports/ moting further understanding on the spreading of co-existing diseases, where facilitation or competition processes are in place. With both Hamiltonian and Lagrangian formalisms secured, canonical transformations become available. These transformations are particularly useful to highlight properties of the dynamical systems and to solve them. They change the old variables (q, p) into new variables (Q, P), while preserving Hamilton's equations. There are a large number of transformations available: it would render impossible to cover all of them here. Instead, we show that at least one canonical transformation exists, and that it promotes the interpretation of the stochastic spreading process as effective mechanical systems. Consider: P 1 (t) = 2p 1/2 q and Q 1 (t) = −p 1/2 . The Poisson bracket {Q 1 , P 1 } q,p = (∂Q 1 /∂q)(∂P 1 /∂p) − (∂Q 1 /∂p)(∂P 1 /∂q) = 1 shows the transformation is canonical. Setting m = 2, the Hamiltonian in terms of the canonical variables (Q 1 , P 1 ) becomes One may interpret − 1 as the Hamiltonian of an effective mechanical problem in one-dimension, in which the particle has mechanical momentum P 1 (τ), with generalized coordinate Q 1 (τ), subjected to a velocity dependent potential.

Conclusion
The description of several real-world problems often includes stochastic fluctuations. The SIS epidemic model includes them due to uncertainties associated with pathogen transmission. For small fluctuation amplitudes, 〈ρ(τ)〉 and σ 2 (τ) are adequate descriptors. Our findings demonstrate 〈ρ(τ)〉 and 1/σ(τ) are conjugated variables, and they satisfy Hamilton's equation. These results link the stochastic SIS epidemic model with a pure dynamical system, which can be solved and manipulated using standard analytical tools. We find the Hamiltonian is a constant of motion for N  1. However, finite size effects break the temporal symmetry of the system:  1/2 τ − follows a power-law around the outbreak onset. A clear explanation for this scaling is still lacking. The relaxation time τ eff portrays the decay of  until equilibrium sets in, meaning that it can also be used to characterize the SIS epidemic. Unlike popular estimates of epidemic growth rate, τ eff remains constant along time and can be extracted from data values of . Finally, our results also suggest a way to incorporate interactions into the SIS model via the Lagrangian function. This finding has intriguing implications for our understanding of facilitation-competition mechanisms between co-occurring diseases since it does not replicate the canonical procedure to obtain Langevin equations.