Improved susceptible–infectious–susceptible epidemic equations based on uncertainties and autocorrelation functions

1. Introduction

Communicable diseases are health disorders caused by pathogens transmitted from infected individuals to susceptible ones [1]. In general, the transmission process occurs with variable success rate, subjected to stochastic uncertainties during the infectious period of the host. These uncertainties comprehend aspects related to biological transmission mechanisms and availability of adequate contact between hosts and susceptible individuals. For large and well-connected populations, stochastic factors are discarded in favour of deterministic differential equations, also known as compartmental or mean-field equations [2–4]. Recent advances in network theory [5] provided a far more clear picture of interactions among elements of the population, improving predictions for heterogeneous social structures. Generalizations for compartmental equations have been able to reproduce pandemics and prove analytical results, taking into account complex network topologies, highlighting the role of central hubs in general disease spreading dynamics [6–11].

By contrast, the stochastic nature of disease transmission cannot be omitted for a number of scenarios. It becomes more pronounced for small populations, where the characteristics of each agent forming the population are relevant variables to the spreading process [12]. Incidentally, this is often the case in emerging diseases [13]. Because the population cannot be treated as homogeneous, average values are no longer adequate, impacting the accuracy of compartmental equations. Stochastic models deal with the issue by proposing simpler rules to express the disease transmission, taking the relevant stochastic factors into account. More importantly, the stochastic analysis expands the machinery used to study the problem beyond population averages. It includes tools such as correlations [12] and autocorrelation functions, which extract inner details of the stochastic dynamics and subsequently provide insights to solve them. For instance, in the standard Brownian motion, the autocorrelation function of the position displays a delta-like behaviour due to white noise, i.e. ⟨x(t)x(t′)⟩ ∝ δ(t − t′). This means that the position of a particle x(t) is uncorrelated to its position x(t′) at time t′, except when t = t′. This feature leads to the well-known linear growth of the spatial variance with time [14]. In disease spreading, autocorrelation functions have also been used to study time series of epidemiological data and assess the impact of spatial influences on stochastic fluctuations [15–20].

Here, we derive exact differential equations for both the instantaneous average density of infected agents, ⟨ρ(t)⟩, and its corresponding variance, σ²(t), in the susceptible–infectious–susceptible (SIS) epidemic model with N agents. We find that uncertainties play an important role in small populations or small prevalence of the disease, impacting estimates of epidemiological parameters from data. Numerical and analytical evidence allow us to formulate two closure relations for ⟨ρ³(t)⟩, and derive systems of differential equations for ⟨ρ(t)⟩ and σ²(t). The selection of the appropriate closure relation depends solely on the nature of the fluctuations present in the system. This issue has been examined in detail before [21–24]. It turns out that the nature of the fluctuations can be assessed from the normalized autocorrelation function D_ρρ(t), including scenarios with finite population sizes. Non-Gaussian fluctuations develop whenever the absorbing state (disease eradication) influences the outcome of disease spreading [21]. We exploit the relationship between D_ρρ(t) and ⟨ρ³(t)⟩ to craft a closure relation in this case. For non-Gaussian fluctuations, a different closure relation emerges as a consequence of vanishing skewness coefficient κ₃(t). The resulting differential equations for Gaussian fluctuations have been reported before [25,26], and also derived in a more general formulation for population dynamics based on Langevin equations [27]. We combine the system of equations into a single nonlinear second-order differential equation, and discuss an analytical solution. The new equations provide significant improvements over the traditional compartmental equation, as they account for stochastic effects, while being far more amenable to analytical studies than the master equation of the disease spreading process.

This paper is organized as follows. Section 2 opens our discussion with compartmental equations of the SIS model, with emphasis on general aspects of parameter estimation. Section 3 reviews the spreading process under the agent-based approach. Improved differential equations for the SIS model are derived. Analytical and numerical properties of D_ρρ(t) are investigated in §4, leading to dynamics for non-Gaussian fluctuations. Dynamics for Gaussian fluctuations are addressed in §5. In each case, the systems of differential equations are combined producing two distinct second-order differential equations for ⟨ρ(t)⟩. In §6, we present our closing arguments, remarks and potential applications.

2. Compartmental equations

Let ρ(t) be the density of infected agents in a population of size N in the SIS model. In the compartmental approach, the population is assumed to be large, homogeneous and highly interconnected. As a result, agents can be regarded as statistically equivalent, and ρ(t) becomes a good descriptor of the system. The assumptions impose that the system must be, on average, invariant under permutations. The simplest way to satisfy permutation invariance assumes agents connected to each other. Incidentally, this population structure shares the same characteristics as the complete graph [4,28].

The other relevant assumption concerns the transmission mechanism. Because the population is taken as homogeneous, the adequate interaction between infected and susceptible agents occurs with probability proportional to (1 − ρ)ρ. This assumption constitutes the basis for the random mixing hypothesis [3]. At the same time, recovery events are proportional to the infected density ρ. Following this notation, the SIS compartmental equation for ρ(t) reads

where α and γ are the transmission and recovery rate, respectively. Explicit generalizations are available for several different networks [6–8], including complex networks. These special network structures highlight the role of super-spreaders in real-world spreading processes [3].

For data fitting and parameter estimation purposes, it is convenient to consider the relative variation of ρ(t) over time. Rearranging equation (2.1) and defining the steady-state density ρ_eq = 1 − γ/α, we obtain

From epidemiological data, equation (2.2) provides a simple way to extract α and γ by a linear fit. For α ≥ γ, dividing equation (2.2) by ρ(t) and plugging the solution $\rho (t)={\rho }_{\mathrm{eq}}/(1-C_1{\text{e}}^{-{\rho }_{\mathrm{eq}}\alpha t})$, with C₁ = 1 − (ρ_eq/ρ(0)), leads to a simple exponential decay

where the decay rate depends only on epidemiological parameters. Again, equation (2.3) can be used to extract α and γ using a linear fit in logscale.

It should be clear by now that equation (2.2) is an important tool to extract epidemiological parameters. What would be the implications for epidemiological studies if equation (2.2) had additional terms or corrections? Figure 1 displays the values of α⁻¹ (d/dt) ln ρ using equation (2.2) with data from numerical simulations (see Data accessibility and [29] for further details). In this controlled computational experiment, predictions for α⁻¹(d/dt) ln ρ deviate from ρ_eq − ρ. Even more, figure 1 shows that early estimates of epidemiological parameters, typical during the onset of epidemics, underestimate the transmission rate. Since agents are equivalent to each other in this setting, the only remaining source of error is due to the discrete nature of transmission and recovery events. Therefore, inherent stochastic events in spreading processes affect predictions whenever uncertainties cannot be neglected. Thus, it seems reasonable to examine more closely this limitation of the compartmental equations, which are far more familiar to epidemiology practitioners [30].

In what follows, we calculate corrections for equation (2.2) using a stochastic agent-based approach to better grasp the emergence of uncertainties in the SIS model. The approach allows for a direct comparison with numerical simulations, and it does not require the coupling of dynamical equations with an external noise source to mimic fluctuations (Langevin formulation). We also note that the averaging procedure employed in the numerical simulations is equivalent to the ensemble averaging. In an ensemble average, the averages are estimated from a large set of independent realizations of a given stochastic process that share the same initial conditions. One way to create an approximate ensemble from real epidemiological data consists in partitioning the system into smaller subsets that are weakly interacting with each other. This is the basis of ensemble formation in physics [31] and the core hypothesis of metapopulation models [32]. Even so, this is only a coarse representation of the idealized ensemble. The problem is somewhat reduced in numerical simulations as the number of realizations can be increased in exchange for computing time. In short, the advantage of ensemble averaging is that it makes it possible to find equations that describe the general behaviour of stochastic variables—for instance, the diffusion equation for the random walker.

3. Stochastic formalism

In the agent-based approach [33,34], the population consists of N distinguishable agents connected to each other according to a predefined adjacency matrix A (N × N). Each agent (k = 0, 1, · · ·, N − 1) may assume one of two possible health states n_k in the SIS model, either susceptible (n_k = 0) or infected (n_k = 1). Following [35,36], there are 2^N available configurations in the canonical basis |μ⟩, with μ = 0, 1, · · ·, 2^N−1. Configurations are readily extracted from the binary construction μ = n₀2⁰ + n₁ 2¹ + · · · + n_N−12^N−1. As an example, for N = 4, the configuration |0⟩ = |0 0 0 0⟩ represents the infected-free configuration, whereas all agents are infected in |15⟩ = |1 1 1 1⟩.

Here, we treat the disease spreading process as a Markov process. The corresponding master equation reads

where $|P(t)⟩=\sum _{\mu =0}^{2^N-1}{P}_{\mu }(t)|\mu ⟩$ is the probability vector, with P_μ(t) being the instantaneous probability to find the system in the configuration |μ⟩; and $\hat{H}$ is the generator of time translations, given by the following expression:

Operators are assigned the hat symbol to distinguish them from scalars. The operators ${\hat{n}}_k$ extract the health state of the k-th agent, ${\hat{n}}_k|n_0\cdots n_k\cdots ⟩=n_k|n_0\cdots n_k\cdots ⟩$, while ${\hat{\sigma }}_k^{\pm }$ are the usual spin-1/2 ladder operators, i.e. ${\hat{\sigma }}_k^{+}|n_0\cdots 0_k\cdots ⟩=|n_0\cdots 1_k\cdots ⟩$ and ${\hat{\sigma }}_k^{-}|n_0\cdots 1_k\cdots ⟩=|n_0\cdots 0_k\cdots ⟩$, respectively. The main advantage of using equations (3.1) and (3.2) lies in their applicability for arbitrary networks, without further assumptions on the probability distribution.

As an example, consider a system in which the adjacency matrix describes a linear chain with periodic boundary conditions, i.e. A_kℓ = δ_k,ℓ±1 and A_0,N−1 = A_N−1,0 = 1. Clearly, the connections between agents are invariant under translations, so that agents are statistically equivalent. However, the number of connections is now reduced to two instead of N − 1 and violates the random mixing hypothesis. Figure 2 exhibits numerical simulations for the linear chain and predictions using the compartmental equation with effective transmission rate (2/N)α. It should come as no surprise that the predictions become increasingly worse for vanishing γ, since agents can only infect their nearest neighbours, thus introducing correlations. By contrast, the agreement between simulated data and the prediction provided by equations (3.1) and (3.2) is far more accurate, reinforcing their validity for general networks.

Despite the known effects of network structures on the dynamics of epidemics [11,28,37,38], there are instances in which the uniqueness of agents can be a minor concern. In these cases, uncertainties stem from the stochastic nature of disease spreading processes. They produce additional corrections to the dynamical equations, with enhanced effects for a finite population of size N. We set aside the complexities associated with network structures by adopting the complete graph. The complete graph replicates the random-mixing hypothesis because each agent interacts with the remaining N − 1 agents, A_ij = 1 − δ_ij. The choice also allows an adequate comparison with the compartmental equations.

Equations (3.1) and (3.2) can be used to evaluate statistics relevant to the epidemic model. Among them, the average density of infected agents,

where ${\eta }_{\mu }\equiv \sum _k⟨\mu |{\hat{n}}_k|\mu ⟩$ is the total number of infected agents in the configuration |μ⟩. By virtue of equation (3.3), it is clear that the time derivative of ⟨ρ(t)⟩ depends solely on dP_μ/dt. In turn, equation (3.1) states $\text{d}{P}_{\mu }/\text{d}t=-\sum _{\nu }⟨\mu |\hat{H}|\nu ⟩P_{\nu }=-\sum _{\nu }{H}_{\mu \nu }{P}_{\nu }$, which concerns the calculation of the matrix elements H_μν. Although their explicit evaluation exists, we are actually interested in the summation $\sum _{\mu }{\eta }_{\mu }{H}_{\mu \nu }$. The latter can be easily calculated noting that $\sum _{\mu }{\eta }_{\mu }⟨\mu |\sum _k{\hat{\sigma }}_k^{+}|\nu ⟩=({\eta }_{\nu }+1)(N-{\eta }_{\nu })$ and $\sum _{\mu }{\eta }_{\mu }⟨\mu |\sum _k{\hat{\sigma }}_k^{-}|\nu ⟩=({\eta }_{\nu }-1){\eta }_{\nu }$. More specifically, the non-vanishing matrix elements $⟨\mu |\sum _k{\hat{\sigma }}_k^{+}|\nu ⟩$ connect configurations whose number of infected agents differ by one, and the number of possible configurations is N − η_ν. For example, with N = 4 and |ν⟩ = |0010⟩, the configurations in question are |1010⟩, |0110⟩, |0011⟩. A similar argument can be made for $\sum _k{\hat{\sigma }}_k^{-}$, reducing η_μ by one and η_ν matching configurations. Therefore,

with instantaneous variance σ²(t) = ⟨ρ²⟩ − ⟨ρ⟩² (see appendix A for details). A brief inspection of equation (3.4) shows that the correction −ασ²(t) always slows down the growth rate of ⟨ρ(t)⟩. As a result, it directly affects the estimation of epidemiological parameters as shown in figure 1.

We emphasize that the inherent fluctuations of the disease spreading process is summarized by σ²(t) in equation (3.4). An initial uncertainty evolves during the course of the spreading process, restricted by the fact that agents can only be either susceptible or infected, i.e. there is no half infection nor half cure. In a sense, σ²(t) shares the concept of shot noise in condensed matter physics [39]. Moreover, equation (3.4) recovers equation (2.1) for vanishing σ²(t), a situation that often arises for large populations since the relative uncertainty scales with N^−1/2. For small populations or small values of ⟨ρ(t)⟩, equation (3.4) highlights the influence of noise in the spreading process, even if agents are statistically equivalent.

Noting that σ²(t) depends on time, there must exist an additional differential equation for σ²(t). Indeed, the same rationale behind equation (3.4) can be used to find (d/dt)σ², as detailed in appendix A ( in accordance with [26] or [27]). The equation of motion reads

where Δ₃(t) = ⟨ρ³(t)⟩ − ⟨ρ(t)⟩³. Again, numerical simulations support our findings (see figures 3 and 4). Despite the encouraging results, equation (3.5) creates an explicit dependence on the third statistical moment ⟨ρ³(t)⟩, even if o(1/N) corrections are omitted. Formally, we could calculate the differential equation for Δ₃(t) but then we would have to deal with ⟨ρ⁴(t)⟩ and so on, creating a set of hierarchic equations for the statistical moments of ρ(t).

Let us briefly assume that it is possible to estimate a surrogate dynamic for Δ₃(t) by some ingenious method. In this case, equations (3.4) and (3.5) form a system of differential equations for ⟨ρ(t)⟩ and σ²(t). However, Δ₃(t) also measures the fluctuation strength and it can change radically for different sets of parameters (N,γ/α) as long as N remains finite. The inset in figure 4a exhibits the changes in σ²(t) as one reduces N. Holding N fixed and varying γ/α also triggers this phenomenon. Figure 5 provides a concrete example of two distinct behaviours for fluctuations for fixed N: Gaussian and non-Gaussian fluctuations. An existing relationship between Δ₃(t) and the instantaneous coefficient of skewness, κ₃(t), provides a way to investigate symmetric fluctuations [26]. Likewise, the density autocorrelation function provides insights on Δ₃(t) for non-symmetric fluctuations. Since the nature of these two types of fluctuations is so dissimilar, we shall study them separately.

4. Autocorrelation function

Uncertainties described by Gaussian fluctuations are expected to play a significant role in widespread epidemics. However, the situation changes when a small fraction of the agents is infected. The SIS model used in this paper does not account for external infection sources, such as wild animals or immigration; once the number of infected vanishes the spreading process comes to a halt. This constraint means that the absorbing state |0⟩ prevents the occurrence of symmetric probability distributions around low densities. The effect can be found in large populations but it is enhanced in small populations: fluctuations can eradicate the disease. Thus, we need to look for a statistics other than κ₃(t) to model the dynamics of Δ₃(t) for non-Gaussian fluctuations. The statistics should involve, at most, ρ(t) up to the power two; otherwise, it could reintroduce higher statistical moments. In that regard, two-point autocorrelation functions fulfil these requirements.

Let C_ρρ(t) be the instantaneous autocorrelation function between ρ(t) and ρ(t + δt), lagged by a single time window

Here, averages are evaluated by considering samples from an ensemble instead of the usual Fourier transform, as the ergodic hypothesis is unavailable. For Markov processes,

The evaluation of this expression involves the same rationale used for equation (3.4), as detailed in appendix B. Plugging the result into equation (4.1), we find C_ρρ(t) = σ²(t) + αδt[ρ_eq⟨ρ²⟩ − ⟨ρ³⟩] + o(δt²). The crucial information here is the relationship between ⟨ρ³(t)⟩ and C_ρρ(t): provided C_ρρ(t) can be fitted from epidemiological data, it seems plausible to use it to model ⟨ρ³(t)⟩ and, thus, create a surrogate dynamics for Δ₃(t). Unfortunately, the lack of a simple functional form prevents the fitting of C_ρρ(t) with at most two parameters.

Instead, consider the normalized autocorrelation function

For vanishing σ²(t) and N ≫ 1, D_ρρ(t) ≈ ρ_eq − ⟨ρ(t)⟩ recovers the r.h.s. of equation (2.2). Hence, D_ρρ(t) can be associated with (d/dt)ln⟨ρ⟩ in the same limit.

According to equation (2.3), an exponential decay of D_ρρ(t)/⟨ρ(t)⟩ occurs whenever ⟨ρ(t)⟩ is reasonably described by compartmental equations. As the system evolves, D_ρρ(t)/⟨ρ(t)⟩ experiences a strong divergence (figure 6). Afterwards, D_ρρ(t)/⟨ρ(t)⟩ either converges to a constant value; or engages in a regime of exponential growth (figure 7). The first case signals that ⟨ρ(t)⟩ describes the spreading process with uncertainties summarized by σ²(t). Fluctuations that increase ρ(t) are as likely as those that decrease it. Thus, the probability density function (pdf) associated with of the fluctuations of ⟨ρ(t)⟩ is symmetrical. We call them Gaussian fluctuations for the lack of a better name.

By contrast, an exponential growth of D_ρρ(t)/⟨ρ(t)⟩ exposes the influence of the absorbing state on the evolution of the system. Its impact becomes more noticeable as ⟨ρ(t)⟩ approaches the disease eradication, and for small population sizes. In such cases, the fluctuation pdf becomes asymmetrical, resulting in the degradation of ⟨ρ(t)⟩ in contradiction with equation (2.1). The fluctuations, in this case, are non-Gaussian. Therefore, D_ρρ/⟨ρ⟩ separates fluctuations into two distinct classes: Gaussian and non-Gaussian.

One could argue that a reciprocal timescale ${\tau }^{-1}=(1/2)(\text{d}/\text{d}t)\ln (D_{\rho \rho }^2/{⟨\rho ⟩}^2)$ emerges because the exponential decay becomes the dominant mode of ⟨ρ(t)⟩ ∝ ρ₁e^−t/ξ, after some time instant t, with peak value ρ₁. Table 1 exhibits a few estimates for τ⁻¹ and ξ⁻¹ from which one can infer τ = ξ/2. As a result, D_ρρ(t) ∝ e^t/λ with λ = ξ after non-Gaussian fluctuations are in place.

By virtue of equation (4.3), we now exploit the relationship between Δ₃(t) and D_ρρ(t) to propose an equation for the expected dynamics of non-Gaussian fluctuations in equation (3.5):

where s(t) = [(2 − ⟨ρ⟩ − ρ_eq)⟨ρ⟩/2] − σ²(t). Equation (4.4) agrees well with simulated data for the entire time interval considered (figure 8) However, the same agreement is not observed for the approximate formula D_ρρ(t) = −D₁e^t/ξ for the entire time interval. In fact, away from the non-Gaussian regime where the fit is accurate, most of the data fall off the proposed curve. Figure 9 explains the reason: D_ρρ[ρ] = a⟨ρ⟩⁻¹ + b is a straight line that intercepts the origin only after the time interval enclosed by the rectangle. The width of the segment shrinks with increasing values of N. Prior to this interval, the curve D_ρρ[ρ] slightly deviates from a straight line, with non-vanishing intercept b. Although an estimate of D_ρρ[ρ] could be useful in this regime, a far more accurate calculation can be obtained by different means, as we show in the next section.

For practical purposes, one can either monitor how $D_{\rho \rho }^2(t)/{⟨\rho (t)⟩}^2$ evolves along time or D_ρρ as a function of ⟨ρ⟩⁻¹. Both methods capture the transition between fluctuation regimes. In the non-Gaussian regime, one should use equations (3.4) and (4.4). Approximations for D_ρρ(t) should be used with care as indicated in figure 8.

For the sake of completeness, we derive the second-order differential equation for ⟨ρ(t)⟩. Taking the time derivative of equation (3.4) and using equation (4.4), one arrives at the desired expression

Results show an excellent agreement with simulated data, regardless of fluctuation type (figure 10). Furthermore, one can employ the approximation D_ρρ ≈ −D₁/⟨ρ(t)⟩, with D₁ ≥ 0 for fixed N and epidemiological parameters as well. Under this assumption, agreement is observed only in the non-Gaussian regime, as expected. It is instructive to study equation (4.5) when o(1/N) corrections are neglected

The characteristic equation provides a coarse estimate for

This expression allows one to quickly grasp the dependence between ξ_est and the parameter D₁. However, there are several issues with ξ_est. The most important one deals with the hypothesis that o(1/N) terms contribute less than other terms in equation (4.5). In fact, they are similar in magnitude and should not be discarded. A far more reliable estimate can be obtained assuming σ²(t) can be written as a power series, i.e. ${\sigma }^2(t)\approx \sum _{m=1}{\sigma }_m^2{\text{e}}^{-m(t/\xi )}$. Collecting only terms proportional to e^−t/ξ, one deduces D₁ in equation (4.7) should be replaced by $D_1+(1/N)[({\sigma }_1^2/{\rho }_1)-(2-{\rho }_{\mathrm{eq}})/2]$. For instance, numerical data suggest ξ_est = 0.201 (N = 20 and γ/α = 0.9), with D₁ = 0.065, ρ₁ = 0.25, σ₁ = 0.063.

5. Gaussian fluctuations

For large population sizes N ≫ 1, stochastic effects are well represented by Gaussian fluctuations and dominated by finite second moments. Noting that the skewness coefficient κ₃ = (Δ₃ − 3⟨ρ⟩σ²)/σ³ vanishes for Gaussian distributions, we conclude ${\mathrm{\Delta }}_3^{\mathrm{gauss}}\approx 3⟨\rho (t)⟩{\sigma }^2(t)$. Indeed, figure 5 shows the ansatz is not too far-fetched since ${\mathrm{\Delta }}_3(t)-{\mathrm{\Delta }}_3^{\mathrm{gauss}}(t)\sim o({\sigma }^2/N)$ for ratios γ/α = 0.1 and 0.5 for N = 50.

Ignoring o(1/N) corrections in equation (3.5), the following differential equations are obtained:

and

Both equations have been derived previously (see [26,27]). As long as σ²(0) > 0, uncertainties play a role in the SIS epidemic model; σ²(0) = 0 implies σ²(t) = 0 and warrants the validity of equation (2.2). Thus, the instantaneous factor σ²(t)/⟨ρ(t)⟩ in equation (5.1a) improves compartmental predictions if σ²(t) ≠ 0. Figure 11 portrays the corresponding direction field.

Despite the insights provided by equations (5.1a) and (5.1b), some issues still remains. The most relevant one deals with the evaluation of σ²(0) from real epidemiological data. In essence, σ²(0) encapsulates the measurement of ignorance about the system at t = 0. In practice, one would rely on clever measurements—possibly, with bias—to estimate σ²(0). Alternatively, the issue can be avoided entirely by combining the system of differential equations for ρ(t) and σ²(t) into a single differential equation

Recalling that setting σ²(t) = 0 is equivalent to using compartmental equations, one can borrow inspiration from projective transformations and rational functions to search for solutions of equation (5.2) starting from equation (2.1). More specifically, assume

with $x(t)={\text{e}}^{-\alpha {\rho }_{\mathrm{eq}}t}$. The coefficients a_k and b_k are obtained for fixed integer m, in consonance with the three critical points discussed earlier, and with the solution of compartmental equations (case σ² = 0). A suitable candidate is

Plugging the expression above in equation (5.2) and solving the coefficients, one obtains two solutions in addition to the trivial solutions. The analytical solution that encircles all values (ρ, σ²) below the separatrix, for α > γ, reads

with constants a₁ and b₂ determined by initial conditions ρ(0) and (dρ/dt)_t=0. Note that σ²(t) can be computed from equation (5.1a), ${\sigma }_1^2(t)={⟨\rho (t)⟩}_1^2(a_1^2-b_2){\text{e}}^{-2\alpha {\rho }_{\mathrm{eq}}}$. The constraint ${\sigma }_1^2(t)\ge 0$ implies $a_1^2-b_2\ge 0$, while the solution of the compartmental equation is obtained by setting $b_2=a_1^2$.

The remaining solution is

with b₁ = [ρ_eq/ρ(0)] − 2. It corresponds to the case ${\sigma }_2^2(t)={\rho }_2^2(t)$ and includes the third critical point $({\rho }_{\mathrm{eq}}/2,{\rho }_{\mathrm{eq}}^2/4)$, along the separatrix. The role of the separatrix can be understood in terms of the signal-to-noise ratio s(t) = ⟨ρ(t)⟩²/σ²(t). Below the separatrix, s(t) > 1 and the average ⟨ρ(t)⟩ becomes more relevant than σ²(t), leading to the equilibrium density ρ_eq. At the separatrix, s(t) = 1 and it indicates that both signal and noise are present in equal measures. Indeed, at the critical point $({\rho }_{\mathrm{eq}}/2,{\rho }_{\mathrm{eq}}^2/4)$ one would expect ⟨ρ(t ≫ 1)⟩ to fluctuate around ρ_eq/2, confined between the other critical points. For large ρ_eq, it also means large deviations. By contrast, noise becomes predominant for s(t) < 1, leading to non-biological dynamics as illustrated in figure 11. Perhaps one can argue s(t) < 1 implies some of the samples used to calculate ⟨ρ⟩ acquire negative values. In this case, the Gaussian description becomes inadequate to portray the biological system. To reinforce this conclusion, one can consider the limiting case with vanishing ⟨ρ(t)⟩ but finite σ²(t): equation (5.1b) is approximated by dσ²/dt ≈ 2ρ_eqσ² so that ${\sigma }^2(t)={\sigma }^2(0){\text{e}}^{2{\rho }_{\mathrm{eq}}t}$. Therefore, the negative parcel in equation (5.1a) grows exponentially along time, producing negative solutions.

The main point of equation (5.2) relies on its compatibility with day-to-day epidemiological data, usually built upon the number of infected patients within a fixed time window. Furthermore, mathematical properties of equations (5.1–5.5) lie well beyond the scope of this paper and merit a proper discussion elsewhere.

6. Conclusion

We have investigated the effects of uncertainties in the SIS epidemic model, finding new differential equations for the average density of infected agents, ρ(t), and its corresponding variance, σ²(t). Our findings reconcile the simplicity of canonical compartmental equations with the accuracy of agent-based simulations, creating suitable tools for practitioners of epidemiology and related fields. At the core of this research, we have demonstrated that uncertainty cannot be neglected in the SIS epidemic model whenever the discreteness of the population is important, even when the population comprises statistically equivalent agents. Uncertainties are inherent aspects of stochastic spreading processes, and their time evolutions are key elements to describe how the number of infected agents vary along time. Concerning their nature specifically, numerical simulations in fully connected networks reveal that uncertainties can be organized into two broad classes, namely, Gaussian and non-Gaussian fluctuations. Gaussian fluctuations, also known as symmetric fluctuations, dominate the spreading process whenever ⟨ρ(t)⟩ and σ²(t) are sufficient to describe the outbreak. This scenario implies the skewness coefficient vanishes for large N, producing a simplified system of differential equations for ⟨ρ(t)⟩ and σ²(t). Alternatively, the differential equations can be combined into a second-order differential equation for ⟨ρ(t)⟩, avoiding problems due to poor estimates of initial values of σ²(t) from raw data. Non-Gaussian fluctuations are far more complex to assess, as they emerge as a consequence of large recovery rates in small population. More specifically, the stochastic process tends to perceive the influence of the absorbing state ρ = 0, creating asymmetric fluctuations. As a result, the skewness coefficient does not converge to a simple mathematical form. Instead, differential equations for ⟨ρ(t)⟩ and σ²(t) are written using the normalized autocorrelation function D_ρρ(t). This function is relevant for the spreading process because it can be interpreted as the likelihood of adequate contact between a given infected agent with susceptible ones, for vanishing variance and large population sizes. For non-Gaussian fluctuations, our numerical simulations show that D_ρρ(t) is proportional ⟨ρ(t)⟩. Therefore, the spreading process reduces, again, to a closed system of differential equations for ⟨ρ(t)⟩ and σ²(t) (see equation (4.4)). Finally, we stress that this research evaluates the impact of uncertainties only for homogeneous populations, i.e. connections between agents are described according to the complete graph. An intriguing question is left open concerning the role of uncertainties in disease spreading processes in other networks.

Numerical simulations are performed using the direct Monte Carlo method [40]. It shares the same origins as the Gillespie algorithm [41], differing only on time step selection. In the Gillespie algorithm, the time step is a random variable distributed according to an exponential pdf; and the system always suffers a single modification once the time step is selected, given by the off-diagonal elements of the transition matrix. In the direct Monte Carlo method, the time step is fixed and equal in value to the average time step of the Gillespie algorithm, both consistent with the Poisson hypothesis. Furthermore, after a single time step has elapsed the system has a chance to remain in the same configuration (diagonal elements of the transition matrix), in addition to the off-diagonal transitions. Direct Monte Carlo simulations tend to be slower than Gillespie but allow for a simple evaluation of statistics at discrete time, including their derivatives and autocorrelation functions, without additional processing algorithms or interpolation methods.

Supplementary Material

Ethics

This work does not present research with ethical considerations.

Data accessibility

Data and numerical codes are deposited at https://doi.org/10.17605/OSF.IO/WFCEP.

Authors' contributions

A.S.M., G.C.C. and G.M.N. designed the research; G.M.N. performed the research and wrote the paper; G.C.C. and A.S.M. edited the paper. All authors reviewed the manuscript.

Competing interests

The authors declare no competing interest.

Funding

This work is supported by Capes 88887.136416/2017-00 and CNPq 307948/2014-5.