Stochastic fluctuations of the transmission rate in the susceptible-infected-susceptible epidemic model

We analyze the dynamics of the susceptible-infected-susceptible epidemic model when the transmission rate displays Gaussian white noise fluctuations around its mean value. We obtain analytic expressions for the final size distribution of infectives and the mean infected number of individuals. The model displays a variety of noise-induced transitions as a function of the basic reproductive rate and the noise intensity. We derive a threshold criterion for epidemic invasion in the presence of external noise and determine the mean epidemic duration and the mean epidemic extinction time.


I. INTRODUCTION
Epidemics are a set of complex processes based on the transmission of a disease between the individuals of a population [1].Mathematical models have become important tools in analyzing the spread and control of infectious diseases and are used in comparing, planning, implementing, evaluating, and optimizing various detection, prevention, therapy, and control programs [1][2][3].They provide flexible and useful tools to study aspects such as passive immunity, gradual loss of immunity acquired either by vaccine or disease, stages of infection, vertical transmission, disease vectors, macroparasitic loads, age structure, social and sexual mixing groups, spatial spread, vaccination, quarantine, and chemotherapy.Specific models have been formulated for diseases such as measles, rubella, chickenpox, whooping cough, smallpox, malaria, gonorrhea, and HIV/AIDS, to name only a few.The choice of which epidemiological classes (compartments) to include in a model depends on the characteristics of the particular disease being modeled and the purpose of the model.Early models are largely based on deterministic equations.They represent a mean-field description that considers the behaviour of averaged quantities and assumes homogeneous mixing of the population.Their justification lies in the assumption that the stochastic fluctuations (demographic noise) around the averaged quantities tend to zero as the system size grows.A mean-field description fails in systems that display near extinction events, where a class of the population becomes so small that the stochastic effects become important and can no longer be neglected.
Demographic noise or internal fluctuations are due to the discrete nature of the constituents of the system (i.e., the susceptible and infected individuals).The effects of such fluctuations on the dynamics of epidemics have been widely explored previously in the literature of the field (see, e.g, Refs.[3][4][5][6][7][8]).Much less attention has been paid, however, to the effects of external fluctuations, which account for other sources of stochasticity such as environmental fluctuations due to changes in meteorological conditions, variability within and between individuals, and other factors.External noise can be taken into account by replacing the characteristic parameters of the deterministic mean-field equations, such as birth and death rates [9] or the transmission rate [10], for example, by random quantities.In this case, the external noise appears multiplicatively in the equations and is able to modify dramatically the mean dynamical behavior of the population [11].
We study the role of external noise in the transmission rate of a susceptible-infected-susceptible (SIS) model.The origin of the external fluctuations lies, for example, in the varying susceptibility to infection of an individual due to circadian rhythms and daily activities, varying susceptibilities between individuals, and variability in the number of contacts between infected and susceptible individuals.These random variations occur on a time scale of hours or less, which is short compared to the time scale of the epidemic dynamics.Therefore the fluctuations in the transmission rate can be modeled by a white noise process [11][12][13].
The paper is organized as follows.In Sec.II we introduce the SIS model with external noise.We analyze the initial dynamics of the infectious disease in Sec.III and define the basic reproductive rate in the presence of random fluctuations.Section IV deals with the behavior of the fluctuations in the number of infected individuals near the boundaries of the state space and establishes the criterion for extinction.In Sec.V we determine analytically the probability distribution of the final size of the epidemic and three of its important characteristics, namely the most likely number of infectives, the mean number of infectives, and the variance of the number of infectives.We determine the mean first passage for the epidemic to evolve from a small initial fraction of infectives to the mean final size as well as the mean epidemic extinction time in Sec.VI.We present concluding remarks in Sec.VII.
The SIS model is the simplest description of the population dynamics for bacterial and other diseases, such as many sexually transmitted diseases [1], where following treatment an infectious individual recovers from the disease but is again susceptible to infection.The system of equations reads ) where S is the number of susceptible and I the number of infected individuals (i.e., individuals able to transmit the disease).The transmission rate is denoted by β and the recovery rate by γ .Adding (2.1a) and (2.1b), we find that the total population number N = S + I is constant, and (2.1b) can be written as We include external fluctuations in the model by allowing the transmission rate β to fluctuate around the mean value β 0 such that denote averaging over the realizations of the noise.)To begin exploring the effects of a fluctuating transmission rate on the dynamics of infectious diseases, we consider the case that is most tractable analytically.We assume, as discussed above, that the transmission rate fluctuates rapidly compared to the evolution of the disease, so that the noise term ζ (t) can be modeled by Gaussian white noise, where D is the noise intensity, a measure of the amplitude of fluctuations in the infection rate with respect to the mean value β 0 .The frequency spectrum of a real noise with a short correlation time is nearly flat [12].
where the noise enters multiplicatively.
Internal fluctuations can also be described by a Langevin equation with multiplicative noise (see, e.g., Refs.[12,14,15]).The functional form of the noise term is different, and more importantly, the noise intensity scales as N −1 .Note that the noise intensity D for the case of external fluctuations scales as N 0 .The effects of internal fluctuations and external noise can be differentiated by their scaling behavior.If the population, that is, N , is large enough, then external noise dominates the dynamics of the epidemic.
Throughout this paper we will interpret (2.5) as a Stratonovich stochastic differential equation, where W (t) is the standard Wiener process and • denotes the Stratonovich integral [11,16].We are interested in analyzing the conditions under which the Langevin equation (2.5), or the stochastic differential equation (2.6), predicts an epidemic invasion and an endemic state and in finding the final size of the epidemics, that is, the stationary probability density P st (I ).
It is also useful to determine three important characteristics of the final size distribution, namely, the mean number of infected individuals I ∞ , the variance of the number of infectives, and the most likely number of infectives in the endemic state.An equation similar to (2.6), though interpreted as an Ito equation, was considered in Ref. [17].That study focuses mainly on the existence and uniqueness of the solution and does not address most of the topics we deal with here.In particular no analytic expression for the final size distribution is derived.Further, since Gaussian white noise is an idealization of real noise with a short correlation time, we prefer the Stratonovich interpretation.As shown by Wong and Zakai and others [12,18,19], if one considers a differential equation with a noise term with some nonzero correlation time, then in the limit that the noise becomes δ correlated (i.e., has a flat spectrum) the differential equation becomes a Stratonovich stochastic differential equation.Consequently, Langevin equations that describe system with real noise should be interpreted as a Stratonovich equation.While Ito and Stratonovich equations often, but not always, predict the same qualitative noiseinduced behavior, the thresholds where transitions occur are generally different [11].

III. INITIAL BEHAVIOR
To find the conditions under which the system evolves into an endemic state, we linearize (2.5) around the infection-free steady state (S = N,I = 0) and obtain Integration of (3.1) yields where Averaging (3.2) we obtain For any Gaussian random variable z with zero mean exp(z and (3.4) turns into Equation (3.6) implies that the intrinsic growth rate μ is given by The basic reproductive rate R 0 in the deterministic SIS model is defined by where R 0 = β 0 N/γ .If we assume by analogy that this definition also holds in the mean in the presence of noise, then In Figure 1 we plot the intrinsic growth rate from (3.7) and compare it to numerical simulations.Note that even in cases where R 0 < 1 (i.e., no invasion occurs in the deterministic case) the epidemic can grow initially, μ > 0, because of the noise, leading to a transient epidemic advance.Equation (3.7) implies that if R 0 < 1, noise will induce initial epidemic growth if In Fig. 2 we plot the number of infected individuals as a function of time by integrating the Langevin equation (2.5) numerically.Each color represents a different stochastic realization.As the noise intensity D increases, the fluctuations around the mean number of infected individuals also increase.The strength of the fluctuations can be quantified by the dependence of the standard deviation of the number of infected individuals on the noise intensity D. This analysis will be carried out Sec. V.

IV. BEHAVIOR NEAR THE BOUNDARIES I = 0 AND I = N
As the noise intensity increases, the noise causes excursions of the infectives population that approach closer and closer to the lower boundary I = 0 (i.e., epidemic extinction) (see Fig. 2).Excursions do not appear to approach the upper boundary I = N as closely.It is therefore useful to investigate the nature of these boundaries.
Consider the general Stratonovich stochastic differential equation or the equivalent Ito stochastic differential equation, where f (X) = f (X) + Dg (X)g(X).In order for the model to be acceptable, the deterministic dynamics must not drive the population of infectives out of the interval [0,N ] [i.e., the kinetic term f (I ) must be such that f (0) 0 and f (N ) 0].Both conditions are fulfilled for (4.4), f (0) = 0 and f (N ) = −γ N < 0. It is equally desirable that the external noise does not drive the population of infectives out of the interval [0,N ] (i.e., we require I = 0 and I = N to be intrinsic boundaries of the stochastic differential equation A boundary b is called natural, if it is attained only with probability zero even if time goes to infinity.The boundary is natural if diverges [11].Here b * is near b, and From (4.4) we find that where c is a constant that depends on b * and ) The integral diverges if and only if , we find that the boundary I = 0 is natural (i.e., unattainable at all times) if R 0 > 1.This means that if R 0 > 1 the epidemic never completely goes to extinction, regardless of the noise intensity.This condition coincides with the criterion for invasion in the deterministic model.
If α 2 + 1 > 0 (i.e., R 0 < 1) we need to consider the quantity L 2 .The boundary b is attracting if L 1 (b) is finite and diverges [11].The behavior near an attracting boundary is as follows: If the random process X(t) starts at t = 0 at x 0 ∈ (b,b * ), then it either leaves this interval in a finite time via the right point, or it never leaves this interval and then X(t) → b as t → ∞.We find which always diverges.The boundary I = 0 is an attracting boundary for R 0 < 1.Since f (0) = 0 and g(0), the boundary I = 0 is also a stationary point of (2.5), or (2.6).The fact that it is attracting for R 0 < 1 implies that in that case the stationary probability will be concentrated entirely on zero, If the deterministic basic reproductive rate is less than one, the final size of the epidemic is zero and the epidemic becomes extinct, even in the presence of Gaussian white noise fluctuations of the transmission rate.Though the epidemic will initially grow in the presence of noise if R * < R 0 < 1, the dynamics are unable to establish an endemic state under these conditions.
To determine the nature of the boundary I = N , we need to evaluate Carrying out the variable transformation ξ = N − I , we obtain The integrand goes to infinity exponentially as ξ → 0, and L 1 (N) = ∞ for all values of the parameters.The upper boundary is always a natural boundary, as expected since f (N ) < 0 and g(N ) = 0.

V. FINAL SIZE OF THE EPIDEMIC
To determine the final size of the epidemic we employ the Fokker-Planck equation, an evolution equation for P (I,t), the probability that the number of infected individuals is I at time t.The following Fokker-Planck equation corresponds to the Langevin equation ( 2 The final size distribution of the epidemic is given by the stationary solution P st (I ) of (5.1).If I = 0 and I = N are natural boundaries, then where C is a constant to be determined from the normalization condition N 0 P st (I )dI = 1.Introducing the new variable x = I/N, (5.2) can be rewritten as ( where χ = 1/Dγ R 2 0 and (•) is the Gamma function.The normalization constant Z is well defined (finite) if and only if 1 + α 1 > 0, which is equivalent to the condition that R 0 > 1.As shown in Sec.IV, that is exactly the requirement for I = 0 to be a natural boundary.As also discussed in that section, for R 0 < 1 the final size distribution is concentrated entirely on zero and an endemic state is impossible.In Fig. 3 we show the results given by (5.3) and (5.4) for the stationary distribution for different values of R 0 and provide a comparison with results from Langevin dynamics simulations using (2.5), which shows an excellent agreement.
While the stationary probability distribution provides the most complete characterization of the epidemic's final size, various numerical characteristics of P st (I ), such as the mean and the variance and the local extrema, are often useful in  applications.The mean value of the final size is given by and from (5.2) we obtain where W a,b (z) is the Whittaker function.In Fig. 4 we plot the mean final size of the epidemic I ∞ /N as D increases (right panel) and as R 0 increases (left panel).As shown above, the mean final size is 0 for R 0 < 1 and increases monotonically with R 0 if R 0 > 1.The mean final size decreases monotonically as the noise intensity increases.As we showed in Sec.III (see Fig. 1) the intrinsic growth rate is positive in the presence of noise for R * < R 0 < 1.The epidemic grows initially, but approaches extinction as time goes to infinity.The noise can cause a transient epidemic in situations where R * < R 0 < 1, but it cannot sustain an endemic state.However, although the presence of external noise does not modify the epidemic invasion threshold, its intensity affects the mean epidemic size.
Figure 2 shows that the fluctuations around the mean value of the number of infected individuals increase with increasing noise intensity.This dependence can be quantified by determining the variance of the final size distribution.We calculate the variance using (5.3) and (5.4) and find (χ ) (χ ) . (5.7) If the noise intensity is small, D 1, the fluctuations depend linearly on the noise intensity, (5.8) In Figure 5 we compare the variance σ 2 obtained by numerical simulations, the exact result (5.7), and the approximate value (5.8).The fluctuations around the mean value increase indeed linearly with the noise intensity in the small noise regime.
FIG. 5. Variance of the number of infected individuals.Solid curves represent the exact analytic result given by (5.7), the broken curves correspond to the approximation (5.8), and the symbols represent the results of numerical simulations.The circles correspond to γ = 0.5 and β 0 = 1, and the triangles to γ = 1 and β 0 = 2.
The local extrema of the stationary probability distribution for the stochastic differential equation (4.1) are the roots of the following equation [11]: (5.9) From (4.4), (4.5), and (5.9) we obtain Consequently either I = 0 or We carry out the variable substitution x = I/N in (5.11), divide both sides by Dβ 2 0 N 2 , and obtain (5.12) The roots of this equation are The roots are real if the discriminant and x m,± are real if R 0 > R d .To represent local extrema of P st (x), the roots x m,± must also be in the interval (0,1).If where then x m,− is negative and x m,+ ∈ (0,1).If R 0 > max{R p,+ ,R d }, then both x m,− and x m,+ are in (0,1).Note that the thresholds R p,± are real and distinct for Dγ < 0.25.They coalesce and become complex conjugates at Dγ = 0.25.As the noise intensity D goes to zero, x m,+ goes to the deterministic steady state, the endemic state, (R 0 − 1)/R 0 , and x m,− goes to minus infinity.If the basic reproductive rate R 0 becomes large (i.e., goes to infinity) then x m,+ → 1 and x m,− → 1/2.
The roots x m,± correspond to the local extrema.To gain a complete picture of the maxima and minima of the final size distribution, we need to examine also the behavior of P st (x) near the boundaries.Equation (5.3) implies that P st (x) → 0 exponentially as x → 1. Near the boundary 0, P st (x) ∼ x α 1 and P st (x) → 0 as x → 0 if α 1 > 0. Otherwise, the final size distribution diverges integrably near zero for R 0 > 1.In other words, x = 0 changes from a maximum to a minimum at which coincides with the condition (5.16) for the roots to pass through zero.We draw the regions of the possible shapes of the final size distribution P st (x) in Fig. 6.In region A, defined by R < max{R p,− ,R d }, the roots x m,± are complex conjugate or negative.The final size distribution P st (x) diverges integrably at x = 0 and decreases monotonically as x increases in this region.In other words, the most likely final size of the epidemic is zero in this region.The thresholds R p,− and R d coalesce at Dγ = 0.22.For Dγ < 0.22, a noise-induced transition occurs at R 0 = R p,− .In region B, x m,+ lies between 0 and 1, while x m,− is negative.The final size distribution is zero at x = 0 and has a maximum at x m,+ .It is only for R 0 > R p,− that the SIS model with a fluctuating transmission rate leads to a most likely nonzero final size.The external noise makes it harder for the epidemic to take hold in the population.At R 0 = R p,+ a re-entrant noise-induced transition occurs.In region C, the final size distribution diverges again integrably at x = 0, and both x m,+ and x m,− lie between 0 and 1.The former corresponds to a maximum and the latter to a minimum.
For 0.22 < Dγ < 0.25, we have a succession of noiseinduced transitions as R increases.At R = R d the final size distribution changes from shape A to shape C. The final size distribution continues to diverge integrably at x = 0, but a local minimum, x m,− and a local maximum, x m,+ appear, giving rise to bistability.At R = R p,− , the bistability and the divergence at zero disappear.The final size distribution P st (x) changes to shape B; it is monomodal with a local maximum at x m,+ .At R = R p,+ a re-entrant noise-induced transition back to shape C, bistability, occurs.For Dγ > 0.25, only one noise-induced transition occurs at R = R d , where the final size distribution changes from shape A to shape C.
Region C corresponds to epidemics with a basic reproductive rate R 0 bounded away from 1. Counterintuitively, the model predicts bistability for larger R 0 .The most likely final sizes are x = 0 and x = x m,+ .The epidemic either dies out or establishes itself in the population.

VI. MEAN FIRST PASSAGE TIMES
To gain further insight into the dynamics of the SIS with a fluctuating transmission rate, we have determined the mean first passage time between two given sizes of the epidemic.In terms of the variable x = I/N, the stochastic differential equation (2.6) Since the boundaries 0 and 1 are natural boundaries for R 0 > 1, the mean first passage time (MFPT) from y to z, M 1 (y → z), with 0 y z < 1, is given by the expression [20] M and the MFPT from y to z with 0 < z y 1 is given by [20] where φ(x) is given by (4.8) with x = I/N.Equation (6.2) allows us to calculate the mean epidemic duration time, that is, the mean time for the number of infected individuals to evolve from the initial value to the final size.Figure 7 shows the dependence of the MFPT on the parameters of the model.Figures 7(a), 7(b), and 7(c) display the mean epidemic duration computed from an initial size of y = 10 −3 to the mean final size z = I ∞ /N given by (5.6), in terms of the recovery rate γ , the noise intensity D, and the mean transmission rate β 0 .The mean epidemic duration increases with the recovery rate, if the other parameters D and β 0 are fixed [see Fig. 7(a)].If γ increases, then R 0 decreases and the epidemics evolves slower to the final state by virtue of (3.7).Therefore, we expect the mean epidemic duration to decrease with γ .The mean epidemic duration decreases with the noise intensity if one keeps γ and β fixed.By virtue of (3.7), the In panel (d) we plot the mean extinction time versus β 0 .The symbols correspond to the results of numerical simulations, and the solid curves are obtained from the expressions (6.2) and (6.3) for the mean epidemic duration and mean epidemic extinction time, respectively.epidemic evolves faster as D increases, and in consequence the mean epidemic duration is expected to decrease with D. In Fig. 7(c) we show how the mean epidemic duration decreases, as the mean transmission rate increases, with D and γ fixed.If β 0 increases, then R 0 increases too, and the mean epidemic duration is expected to decrease.
The MFPT given by (6.3) provides the mean epidemic extinction time, that is, the mean time for the final epidemic size I ∞ /N to evolve to 10 −3 .In Fig. 7(d) we show the mean epidemic extinction time as function of β 0 fixing D and γ .In this case the dependence is not monotonic.There exists a value of β 0 for which the mean extinction time reaches its maximum value.This is due to opposite effects caused by an increase in β 0 .As we have seen in Fig. 7(c), if β 0 increases the epidemics evolves faster and the mean epidemic duration decreases.At the same time the final size increases with β 0 (i.e., with R 0 ) (see Fig. 4), which implies that the epidemic duration increases with β 0 .We find excellent agreement between theoretical predictions and results from numerical simulations.

VII. CONCLUSION
Rohani et al. have studied the effects of environmental stochasticity, causing the transmission rate to be a normally distributed random variable, in the dynamics of childhood diseases [10].Fluctuations in the transmission rate have also been considered to account for seasonal fluctuations due to rural-urban migration or seasonal migrations [21].The effects of rapid fluctuations in the transmission rate on the dynamics of the SIS model are then expected to play a significant role in dynamics of epidemics, and we have analyzed the role of fluctuations both analytically and numerically.We have defined the stochastic basic reproductive rate based on an analysis of the initial behavior of the infectious disease.That analysis also revealed that external noise can induce initial epidemic growth, even if the deterministic basic reproductive rate is less than 1.We have derived the criterion for extinction and persistence of the epidemic by determining the nature of the boundary I = 0 of the state space.An analytic expression for the final size distribution of infectives is obtained by solving the stationary Fokker-Planck equation and is compared with numerical simulations.
In the presence of Gaussian white noise, there is an endemic state for R 0 > 1.The value of the final size depends on the noise intensity as well as on the transmission and recovery rate.The epidemic becomes extinct, that is, the final size of the epidemic is zero, regardless of the noise intensity if R 0 < 1.If R * < R 0 < 1, the noise induces an initial epidemic growth, but the final outcome is still extinction.We find that various noise-induced transitions occur in the stochastic SIS model.If the noise intensity is sufficiently small, a re-entrant transition occurs, where I = 0 goes from being a maximum to a minimum back to a maximum of the final size distribution as R 0 increases from 1.If the noise intensity surpasses a certain threshold, Dγ = 0.25, then I = 0 is a maximum of the final size distribution for all R 0 larger than 1.As R 0 crosses a threshold, R 0 = R d , the final size distribution develops an inflection point that splits into a minimum and a maximum, which migrate towards 0.5 and 1, respectively, as R 0 becomes large.The final size distribution displays noise-induced bistability.Finally, we have computed the mean epidemic duration, the mean time the system takes to evolve from the initial condition of x(0) = 10 −3 to reach the mean final size, and the mean epidemic extinction time, the mean time the system takes to evolve from the mean final size to x = 10 −3 , and their behavior in terms of the parameters of the model.Our findings and methodology could be useful to study the evolution of specific epidemics where environmental stochasticity can play an important role, as in the case of malaria or pertussis.
We have focused on situations where external noise dominates the dynamics of the epidemic and internal fluctuations are negligible.This is the case in sufficiently large populations.An open problem is the interplay between external and internal fluctuations in populations of small to moderate size, where both types of noise are expected to play a role.Interactions between the two noises could possibly give rise to unexpected phenomena; this is a subject that deserves further exploration in future works.

FIG. 1 .
FIG. 1. (Color online) Intrinsic growth rate as a function of R 0 for different values of the noise intensity D. Symbols represent the results of numerical simulations and solid curves correspond to (3.7).

FIG. 3 .
FIG. 3. (Color online) Stationary probability distribution (i.e., final size distribution, of the number of infected individuals for different values of R 0 and γ = 1).The symbols represent the results of the numerical simulations, and the solid curves correspond to the exact analytic results from (5.3) and (5.4).

FIG. 4 .
FIG. 4. (Color online) Mean final size of the epidemic as function of the basic reproductive rate (left panel) and as function of the noise intensity (right panel).Symbols correspond to the results of numerical simulations and solid curves to the exact results (5.6).(γ = 1).

FIG. 6 .
FIG. 6. (Color online) Diagram for the different shapes of the final size distribution.

FIG. 7 .
FIG.7.Mean first passage time.In panels (a), (b), and (c) we plot the mean epidemic duration versus γ , D, and β 0 , respectively.In panel (d) we plot the mean extinction time versus β 0 .The symbols correspond to the results of numerical simulations, and the solid curves are obtained from the expressions (6.2) and (6.3) for the mean epidemic duration and mean epidemic extinction time, respectively.