Outbreak-size distributions under fluctuating rates

We study the effect of noisy infection (contact) and recovery rates on the distribution of outbreak sizes in the stochastic SIR model. The rates are modeled as Ornstein-Uhlenbeck processes with finite correlation time and variance, which we illustrate using outbreak data from the RSV 2019-2020 season in the US. In the limit of large populations, we find analytical solutions for the outbreak-size distribution in the long-correlated (adiabatic) and short-correlated (white) noise regimes, and demonstrate that the distribution can be highly skewed with significant probabilities for large fluctuations away from mean-field theory. Furthermore, we assess the relative contribution of demographic and reaction-rate noise on the outbreak-size variance, and show that demographic noise becomes irrelevant in the presence of slowly varying reaction-rate noise but persists for large system sizes if the noise is fast. Finally, we show that the crossover to the white-noise regime typically occurs for correlation times that are on the same order as the characteristic recovery time in the model.


I. INTRODUCTION
Epidemic models are useful for understanding the spreading dynamics of general contagious processes and effectively describe a wide variety of phenomena from spreading diseases, to rumors, fads, panics, computer viruses, and even election dynamics [1][2][3][4][5][6][7][8].In addition, epidemic models are practically useful, since epidemiologists rely on models to quantify risks of local epidemic outbreaks of various sizes and formulate optimal control strategies for many diseases including Sars-Cov-2, Ebola, and Dengue [9][10][11][12][13][14]. Within a given population, outbreak dynamics are typically described in terms of compartmental models [1,4,15].For example, starting from some small initial infection, over time, individuals in a population make stochastic transitions between some number of discrete disease states (susceptible, exposed, infectious, etc.) based on prescribed probabilities for a particular population and disease [9,12,[16][17][18][19][20].In the limit of large populations and non-fluctuating parameters the stochastic dynamics approach deterministic (mean-field) differential equations for the expected fraction of a population in each state [1,4,15,21].
Yet, for real finite populations with time-fluctuating parameters outbreak dynamics have a wide range of outcomes for each initial condition, which are not predicted by mean-field models [1,17,19,[21][22][23][24][25].For instance, recently we developed a theoretical approach that allows for calculating the distribution of outbreak sizes in wellmixed populations under demographic noise.This approach provided a closed-form expression for the asymptotic outbreak distribution in the Susceptible-Infected-Recovered (SIR) model and more general SIR model extensions with fixed population sizes (N ) and static infection/recovery rates [13].However, many data analyses have shown that for a multitude of diseases, bestfit epidemic model parameters can fluctuate significantly in time [26][27][28][29][30][31][32][33].For instance, by measuring the relative changes in reported disease incidence and hospitalization, one can compute an effective infectious contact rate between individuals in a population over time.Doing so one often finds fluctuating and/or periodic rates in general [1,31,32,[34][35][36][37][38], which in the case of human epidemics correlate with more general social contact rates [39].For instance, techniques for extracting timedependent parameters have been applied to the recent COVID-19 pandemic as well [31,33,40,41], in order to account for fluctuations in contact rates, rendering the usual SIR class of forecasting models time dependent.In addition, here we give an another example based on 2019-2020 hospitalization data of the respiratory syncytial virus (RSV) season in the US [42], and find the data effectively parameterized in terms of two general metrics for quantifying temporal variations about a mean: the infection rate's standard deviation and correlation time.
Despite the theoretical importance of understanding noise effects in canonical non-equilibrium epidemic models, as well as the practical importance for quantifying uncertainty in real epidemics, a general analytical approach for describing small and large fluctuations in outbreak dynamics due to parameter fluctuations is still lacking.
Here we develop such an approach within the context of the SIR model with noisy reaction rates with finite variances and correlation times.We motivate our use of these standard noise characteristics by extracting infectious contact rate fluctuations in RSV outbreak data from the U.S. in 2019-2020 using a Bayesian model inference.In terms of general model analysis, we focus on the outbreak-size distribution and quantify the probabilities for outbreaks that differ from the mean-field predictions.In particular, we calculate the distribution in the limit of adiabatic and white noises, and demonstrate several important properties including: the skewness of the outbreak distribution toward unusually small outbreaks, and the existence of optimal values of the basic reproductive number that maximize the outbreak variance.We also study the cross-over of the outbreak distribution with finite population size and noise-correlation time and analyze when the limiting theories of demographic, adiabatic and white reaction-rate noise apply.

II. SIR MODEL WITH REACTION-RATE NOISE
We are interested in outbreak dynamics in which the epidemic wave is fast compared to both demographic and re-infection time scales; the latter denotes the possibility for individuals to be infected multiple timescales [43].The canonical epidemic model for this regime is the SIR model [1,2], in which individuals are either susceptible (capable of getting infected), infected, or recovered (or removed/deceased), and can make transitions between these states through two basic processes: infection and recovery.Denoting the total number of susceptibles S, infected I, and recovered R in a population of fixed size N , the probability per unit time that the number of susceptibles decreases by one and the number of infected increases by one is βSI/N (for a well-mixed population), where β is the infectious contact rate [1,2,4].Similarly, the probability per unit time that the number of infected decreases by one is γI, where γ is the recovery rate [1,2,4].As a result, the deterministic rate equations in the limit of large N describing the mean densities of susceptibles, infected and recovered, x s = S/N , x i = I/N and x r = R/N , respectively read: where Starting from a small initial infection density, x i (t = 0) ≪ 1, the final fraction of susceptibles in Eqs.(1) x * s ≡ x 0 satisfies x 0 = e −R0 (1−x0) , where R 0 = β/γ is the basic reproduction number [1,2].Hence, in the mean-field theory the total fraction of the population infected over the whole epidemic wave is x Note that if R 0 ≤ 1 in Eq.( 2) then x * r = 0, giving us the usual condition R 0 = 1 as the epidemic threshold.
As noted in Sec.I, in many cases the parameters for the SIR model are time fluctuating.As a simple model, we allow the infection and recovery rates to be generated by independent Ornstein-Uhlenbeck (OU) processes with some correlation times and variances.For simplicity, we assume the correlation times are identical for both rates and equal τ , while the noise variances are σ 2 β and σ 2 γ for the infection and recovery rates, respectively.Thus, we write: β(t) = β 0 (1 + ξ β (t)) and γ(t) = γ 0 (1 + ξ γ (t)), and augment Eq.(1) into the stochastic system Here, η β and η γ are Gaussian white noises, while ξ β and ξ γ are OU processes.Note that by construction, β(t) and γ(t) are assumed to be widesense stationary Gaussian processes with ⟨β⟩ = β 0 , β e −∆/τ , and ⟨(γ(t) − γ 0 )(γ(t + ∆) − γ 0 )⟩ = σ 2 γ e −∆/τ , where ⟨•⟩ denotes the expectation operator.In general, one can simulate the system of equations ( 3) and find the final outbreak-size distribution for fluctuating rates with any magnitude and correlation time.

A. RSV model fit
Finite correlation time and variance are general physical metrics that quantify temporal fluctuations around a mean -the sort of temporal variation observed in many epidemic data analyses [1, 26-38, 40, 41].We can further motivate our study of the SIR model with temporally fluctuating reaction rates by extracting such noise characteristics from data on the 2019-2020 respiratory syncytial virus (RSV) season in the U.S.
We perform a parameter inference from RSV hospitalization data assuming a discretized version of Eqs.(3) with daily time steps and a fixed recovery rate γ = γ 0 = 1/7 days −1 [44].We use the well-known platform Stan via the R package rstan [45,46] to do the numerical Bayesian inference by tying the dynamical model to the number of recorded daily hospitalizations, as obtained from the CDC [42].The parameters for the inference are: β 0 , the inverse correlation time α, σ β , the hospitalization rate, and initial conditions for the SIR [47]; output examples are shown in Fig. 1.In panel (a) we plot the daily hospitalization numbers and compare to the median prediction of the model (in red) along with its credible intervals.A similar plot is shown in panel (b) for the daily infectious contact rate, which drives the predictions for (a).Further details are given in App.A [48].
Our inference uncovers significant temporal fluctuations in the most-likely RSV infectious contact rate.A summary of the output that is relevant for our analysis includes: R0 = 1.37 in [1.32, 1.44], α = 0.11 days −1 in [0.045, 0.20] days −1 , and σβ = 0.026 days −1 in [0.014, 0.040] days −1 , where ˆdenotes the median within a quartile range specified by the square brackets.From these we observe a fairly tight value of the inferred timeaveraged R 0 , but with substantial temporal fluctuations between 10 − 20%.On the other hand, the correlation time estimate α −1 is quite broad ranging from 5−20 days.Nevertheless, note that such time scales are significantly smaller than seasonal effects, which are expected to occur on time scales on the order of a year.The noise-inference, therefore, quantifies temporal fluctuations distinct from seasonality [49].We can situate the inferred noise characteristics of the RSV season within the results of our analytical theory, see SecV.We begin by analyzing outbreak statistics driven by the fluctuations in Eqs.(3).

III. LIMIT OF ADIABATIC NOISE
In order to gain insight on the outbreak distributions generated from the general Eqs.( 3) and temporal fluctuations of the sort we inferred from RSV data, we first consider limiting regimes.We start with the limit of adiabatic noise, τ ≫ 1.Here, the underlying assumption is that, during the epidemic wave, the rates do not change appreciably, and hence the problem reduces to that of quenched noise on system (1).For simplicity and illustration of the adiabatic limit, here we deal with the case where only β varies and γ is constant, such that σ γ = 0.In order to simplify the equations, we take γ 0 = 1, which merely specifies the time units and results in R 0 = β.
To find the distribution of the final outbreak size P (x * r ), we have to compute the following integral: The conditional probability P (x * r |β) is a dirac delta function around the mean-field value of the outbreak at β.
, where the mean-field final outbreak fraction x * r satisfies Eq.( 2).Taking a Gaussian distribution for P (β) with mean β 0 and standard deviation σ β , Eq. (4) becomes We point out that in order for the SIR model to make physical sense β ≥ 0. Therefore, when plugging in an unrestricted Gaussian in Eq.( 5), σ β cannot be too large [50].
Otherwise other distributions, e.g. that vanish at β = 0 can be used instead; yet this does not change the results qualitatively.We also note that since x * r (β) vanishes for β ≤ 1, the lower boundary in the integral in Eq. ( 5) can be taken to be 1 without affecting the distribution.
Changing variables from β to x * r , and using the fact , we can explicitly perform the integration by plugging instead of β, − log(1−x * r )/x * r , which is the solution of x * r = x * r (β).6) is shown in solid black, while the small-noise limit Eq.( 7) is shown with a dashed line.For both panels γ0 = 1.
As a result, Eq. ( 5) reduces to From Eq.( 6) we can derive, e.g., the typical fluctuations around the mean-field given by σ a -the standard deviation associated with the adiabatic outbreak PDF.In particular, in the limit of small σ β , the variance becomes Note that for adiabatic noise we can repeat our calculation of the outbreak-size distribution and variance for any quenched distribution of β (or γ), and not just Gaussians, which may better reflect the data in a given application.
Next, we can plot the probability distribution function (PDF) for adiabatic infection-rate noise and explore its qualitative features.An example prediction is shown in Fig. 2 (a) with a solid line for fixed values of β 0 = 2 and σ β = 0.1β 0 .The solution from Eq.( 6) can be compared to stochastic simulations of Eqs.(3) for large τ .Note that the agreement with simulations is quite good.Qualitatively, one of the most important features that we observe in the PDFs is the high degree of skewness toward small outbreaks.We can get a quantitative measure of this skewness by examining the exponent of Eq.( 6), called the action (for reasons explained in Sec.IV), for two limiting values of the outbreak fraction: x * r = 0 and x * r = 1, i.e., small and large outbreaks.Indeed, the PDF [Eq.(6)] can be described effectively as P ∼ exp[−S/σ 2 β ], where S = (ln(1 − x * r )/x * r + R 0 ) 2 /2.When x * r → 0 the action remains finite, i.e, S → (R 0 − 1) 2 /2.On the other hand, when x * r → 1, S → ∞.Hence, minimally small outbreaks occur with finite probability for finite R 0 , while maximally large outbreaks can never occur when the reaction-rate noise is finite, which is why the outbreak distribution's tails are skewed toward small outbreaks.
In addition to the PDFs, we can examine the variance of the outbreak PDF for adiabatic noise as a function of R 0 .Examples can be seen in Fig. 2(b), where we plot simulated outbreak variances for three large values of τ with γ = 1.Here another interesting qualitative feature emerges: the existence of a maximum in the outbreak variance for some value of R 0 .On the one hand, as σ β → 0, the maximum approaches R 0 = 1.On the other hand as σ β increases, the maximum variance occurs for an R 0 that is an increasing function of σ β .For example, in Fig. 2 (b) we observed a maximum near R 0 = 1.1.However, the saddle-point equation for the maximum variance in the adiabatic limit cannot be solved analytically.
In general, we observe good agreement with the predicted variance of Eq.( 6) (solid line) and the small-noise limit Eq.(7) (dashed line), including the existence of a maximum, which the former captures.Yet, as R 0 → 1, eventually all the simulation results have discrepancy with both adiabatic predictions.The reason is, as we approach the epidemic threshold, the SIR dynamics slow down, meaning that even a large τ may not be "slow" with respect to the underlying process.

IV. LIMIT OF WHITE NOISE
So far we have assumed that the dynamics of the noise is slow compared to the dynamics of Eqs.(1), but what happens if it is fast?In this latter limit, τ → 0, instead of Eqs.(3) we can write which we denote as the white reaction-rate limit.Here, ζ β and ζ γ are white Gaussian noises.In order to coincide with Eqs.(3) as τ → 0, one must demand that To analyze the outbreak-size PDF given Eq. ( 8), we follow the approach in [13], and construct the equivalent Fokker-Planck equation for the probability to observe densities x s and x i at time t (assuming Itô calculus): To simplify notation, henceforth, we will assume that σ 2 2 = f σ 2 1 , with f > 0 for simplicity, and again rescale time t → γ 0 t, so that β 0 is replaced by the basic reproduction number, R 0 = β 0 /γ 0 .Next, we employ the WKB approximation P (x s , x i ) ∼ exp −S(x s , x i )/σ 2  1 , which is the expected scaling-form for solutions to Eq.( 10) in the limit of small noise and large deviations [51][52][53][54], and which we observe in simulations of Eqs.(3). Figure 3 shows several examples of the expected scaling with noise-variance for different values of the final outbreak size.Indeed, the logarithm of the probability tends to straight lines as 1/σ 2 1 is varied, with slopes that change with the outbreak size.Using this insight, we substitute the exponential ansatz into Eq.(10) and arrive at a Hamilton-Jacobi equation, ∂S/∂t+H = 0, in the leading order in σ 1 ≪ 1, with In this formalism, H is called the Hamiltonian, S is the action, while p s = ∂ xs S and p i = ∂ xi S are the conjugate momenta, just as in analytical mechanics [51,53,54].
In order to compute probabilities for different outbreak sizes, we need to know the action S, which means calculating the integrals S = p s dx s + p i dx i − Hdt given the dynamics of Hamilton's equations: ẋs = ∂ ps H, ẋi = ∂ pi H, ṗs = −∂ xs H, and ṗi = −∂ xi H.We can simplify the action computation by noting that, first, since we are interested in outbreaks that emerge from initially small levels of infection x i (t = 0) → 0 the "energy" of outbreaks is zero, H = 0.As the Hamiltonian (11) has no explicit time dependence, it is a constant of motion, namely H(t) = 0. Second, we can rewrite the Hamiltonian as i , using ẋs and ẋi .Third, by substituting the zero-energy condition into ṗi , we obtain that ṗi As a result, the Hamiltonian (11) simplifies to: H(t) = p s ẋs + (d/dt)(x i p i ).Integrating both sides of this equation with respect to time over the full course of an outbreak, yields 0 = p s dx s + x i (t → ∞)p i (t → ∞) − x i (t = 0)p i (t = 0).As the fraction of the population infected goes to zero both at short and long times (assuming no reinfection), we derive the useful fact that p s dx s = 0.As a consequence, the action associated with an outbreak in the white-noise limit is simply A.
Phase-space trajectories for outbreaks In order to compute the outbreak probabilities we need to solve Hamilton's equations and substitute the resulting trajectories into Eq.(12).To do so, we must understand the phase-space structure of outbreak paths.First we recall that in the mean-field system Eq.( 1), the outbreak dynamics follow a heteroclinic trajectory, which starts at t = 0 at a fixed point (x s = 1, x i = 0) and ends at the final state (x s = x 0 , x i = 0) as t → ∞.In our Hamiltonian system this corresponds to a special trajectory with p s = 0 and p i = 0, or in phase space (x s , x i , p s , p i ) = (1, 0, 0, 0) for t = 0.However, in general there are an infinite number of related x i = 0 initial conditions with non-zero momenta, which one can find by solving ẋs = 0, ẋi = 0, ṗs = 0, and ṗi = 0, given x s = 1 and x i = 0.It is straightforward to show that the general fixed-point initial conditions are where δ ≡ p i (t = 0) is a free parameter.As pointed out in [13] for the case of demographic noise, if we propagate each of the possible initial conditions forward in time, they tend to unique final outbreak values; namely, one x * s (x * r ) for each δ.Examples are shown in Fig. 4 (b), where we plot the outbreak sizes as a function of δ for three different values of R 0 .A simple algorithm for generating the outbreak distribution numerically for a fixed value of R 0 is to: (1) pick a δ, (2) propagate forward with Hamilton's equation given Eq.( 11) (assuming some small perturbation from the chosen fixed point), (3) compute the integral in Eq.( 12) from the resulting trajectory, and (4) repeat for another value of δ.Each δ results in a unique x * s and S(x * s ).The slopes of the lines in Fig. 3, were computed in just this way, and correspond to numerical solutions for the outbreak paths and associated S(x * s ) for the chosen values of x * r = 1−x * s .Similarly, by sweeping over values of δ we can compute the full white-noise distributions for any x * r .Examples are plotted in Fig. 2(a) (the narrow PDF prediction) for infection-rate fluctuations, and in Fig. 4(a) for different combinations of infection and recovery noise.For all Figs.2(a), 3, and 4(a) the white-noise WKB theory and simulations agree well, which demonstrates the accuracy of our general approach.In fact, by combining the method presented with the results of [13], we have a complete algorithmic solution for generating the outbreak PDF of the SIR model with general and multicomponent white noise, which we return to in Sec.V.

B. Outbreak variance
In the general case of noise in both β and γ, it seems that Hamilton's equations cannot be solved analytically in a simple manner -apart from constructing a powerseries expansion in the initial-condition parameter δ.The primary reason for this, in contrast to [13], is that for the reaction-rate noise discussed in this work there is no conservation of momentum.Therefore, we proceed to first calculate the variance of the outbreak-size distribution, which is related to the lowest order contribution to Eq (12) in δ.A complete solution for the case of recoveryonly fluctuations, at all orders, is given in App.B.
For the variance calculation, we attempt to find the action in the vicinity of the mean-field final outbreak fraction 1 − x 0 .First, let us assume p i ≪ 1, to be verified a-posteriori.Equating H = 0, yields: p s = (1 − 1/(R 0 x s ))p i , i.e., p s ≪ 1 as well.Second, we show that p i (t) remains small during the entire epidemic duration as long as the initial momentum δ is small.Writing down the Hamilton's equation for ṗi = −∂ xi H, and using Eq. ( 11) we have: ṗi i .The solution of this differential equation is where p i (0) = δ is the initial condition, and we have used the fact that in the leading order in δ ≪ 1, x r = x i dt = ln(x s )/R 0 .This is legitimate as the action will have a δ 2 dependence, see below, so we can substitute in O(δ 2 ) terms their mean-field O(δ 0 ) approximation.
To compute δ, we can use Hamilton's equations for ẋs , and ẋr = − ẋs − ẋi , and compute ẋr / ẋs = dx r /dx s .This yields a differential equation, which can be solved with initial conditions x r (t = 0) = 0 and x s (t = 0) = 1, assuming that during the epidemic duration, p i (t) is almost constant within O(δ).Using Eq. ( 14) and that when the outbreak ends x * r = 1 − x * s , and assuming x * s − x 0 ∼ O(δ) (to be confirmed a-posteriori), we find This confirms a-posteriori that δ ≪ 1, under the assumption that the final susceptible fraction is close to its meanfield counterpart, i.e. x * s − x 0 ≪ 1.Finally, to compute the integral in Eq. ( 12), it is more convenient to change variables to x s , see Eq. ( 14).Thus, we write: Here, the Jacobian can be found using the Hamilton's equations: where again we have used mean-field results for the O(δ) terms, namely x i = 1 − x s + ln(x s )/R 0 .Putting it all together, and using Eqs.( 14) and ( 15), we can perform the integration in Eq. ( 12) over x s from 1 to x * s , which yields the action, in the leading order in x * s − x 0 ∼ δ: Indeed, having obtained a δ 2 dependence of the action corroborated our assumptions a-posteriori.Here, v = v(R 0 ) is the (rescaled) variance of the outbreak-size distribution.Namely, remembering that we have sought the outbreak-size PDF as P (x * s ) ∼ exp −S(x * s )/σ 2 1 , the variance of the outbreak-size distribution in the limit of white reaction-rate noise, σ 2 w , is We can test our predictions for the outbreak variance in the white-noise regime by performing stochastic simulations of Eqs.(8) with different values of R 0 and different combinations of noise.Results are shown in Fig. 4(b).First, one can see an interesting behavior where the variance receives a maximum at R 0 ≃ 1.33, similar to the adiabatic regime shown in Fig. 2 (b).Here, however, the maximum variance occurs for an R 0 that is independent of the noise amplitude and noise combination, unlike adiabatic noise.The reason for the maximum appearing around 1.33 is that for this value of R 0 the mean outbreak fraction is approximately obtained at x * r ≃ 0.5 which maximizes the variance possibility.In addition, we note that the predicted outbreak variance in the white-noise regime only depends on R 0 and the total variance of the reaction-rate noise, σ 2 1 + σ 2 2 .For example, in Fig. 4(b) we show the predicted scaling collapse to a single function of R 0 of the simulated outbreak variance resulting from different combinations of noise.In general, we would expect infection-rate and recovery-rate noise to produce additive variance (since the two noise sources are independent), but the fact that their prefactor dependence on R 0 is identical is interesting.On the other hand, one can check that this symmetry between infection and recovery noise disappears for higher-order statistics, e.g., by repeating the above calculation to O(δ 3 ) for the third central moment.

V. CROSSOVER WITH CORRELATION TIME AND SYSTEM SIZE
Now that we have analyzed the outbreak-distribution in limiting cases (including App.B), we next address when the various limiting regimes apply.In particular, we examine the cross-over behavior of the stochastic SIR model as a function of the reaction-rate noise correlation time and population size; the latter has been assumed infinite so far.We use as our metric the variance of the outbreak-size distribution since it is the lowest order statistic not captured by mean-field theory.
First, we remain in the N → ∞ limit, and try to understand how small (large) τ has to be in order to produce effectively white (adiabatic) outbreak statistics.To do so, we plot in Fig. 5(a) the variance of the outbreak-size PDF found from simulating Eqs.(3) versus the (inverse) correlation time τ for three values of R 0 and fixed σ β .Note that the outbreak variance is normalized by the adiabatic limit, Eq.( 7), so that each simulation series approaches unity for small τ −1 .In addition to the adiabatic limit, for comparison we plot predictions for white-noise, Eq.( 16), with lines.In the latter case, the τ -dependence comes from the definition of the white-noise variance, Eqs.(9).
Figure 5(a) has several important features.First, we point out that the outbreak variance has a maximum in the adiabatic limit, meaning that for fixed infection-rate noise variance, the SIR model dynamics is most sensitive to slow noise.This effect is observed in other population models as well [55,56].For the SIR model, the primary reason is that even relatively small fluctuations in β can bring an epidemic closer to the R 0 = 1 threshold.If the noise is slowly varying, in particular, the effect is felt over the full time-course of the epidemic wave, resulting in potentially much smaller outbreak than the mean-field.As mentioned in Sec.III, this produces highly skewed PDFs with significant probabilities for small outbreaks, and hence large outbreak variance.In contrast, in Fig. 5(a) we can see that the white-noise predictions (the lines) are accurate for quite large values of τ .In fact, for each value of R 0 = 2, 1.5, and 1.2 (from top to bottom), we can see that the white noise prediction remains valid for correlation times on the order of the recovery time, τ ∼ γ −1 0 = 1.In general, the crossover point in τ between white-and adiabatic-noise regimes has some R 0 dependence: the smaller R 0 , the larger τ can be for the white-noise results to be valid, since effectively as the epidemic gets closer to threshold the dynamics slows down, making even slowly-varying noise potentially fast.An estimate for the crossover time, τ c , can be found by solving σ 2 a = σ 2 w , or where we have used Eqs.( 7) and ( 16), valid for small noise.Evidently, τ c depends only on R 0 and not, e.g., on the noise variance for small noise.The crossover time is plotted in the inset of Fig. 5(a), which for the typical model parameters of R 0 − 1 ∼ O(1) remains near the recovery time scale (or unity in our chosen units).Now that we have an estimate for cross-over times, we can situate the inferred RSV contact-rate fluctuations and determine what regime they fall into.By plugging in the median and quartile inferred parameter values given in Sec.II A into Eq.(18), and using Eqs.( 7) and ( 16), we find that the ratio of the noise correlation time to the cross-over time, τ /τ c , falls between 0.1−0.2.As the τ estimates are substantially smaller than the cross-over times, we expect the outbreak-size statistics to be well approximated by the white-noise theory.Hence, our analytical results can be used to make quantitative estimates for future RSV outbreak size probabilities, assuming parameters remain relatively similar to the 2019-2020 epidemic.
The second crossover that we consider is that of finite system size.Namely, how large does a population have to be before demographic noise becomes irrelevant compared to reaction-rate noise?For this exploration we perform a discrete time stochastic simulation (with small time steps) of the discrete state reactions defined for the SIR model in Sec.II (above Eqs.1) while the reaction-rates fluctuate according to the OU processes in Eqs.(3).In Fig. 5(b) we plot the outbreak-size variance as a function of N −1 for three values of the correlation time τ = 10, 1, and 0.1 (from top to bottom).The curves are the expected total white-noise variance, σ 2 w,tot , which is a sum of reaction-rate and demographic noise The demographic piece was calculated in [13].Note that here we have assumed that the total variance is the sum of the variances from the independent noise sources [57], and that this result holds for σ 2 1 , σ 2 2 , N −1 ≪ 1.
In Fig. 5(b) we can see that for large system sizes the variance becomes flat with respect to changes in N and approaches approximately the expected white-noise limit, Eq.( 16) -especially for the two smaller values of τ where the white-noise approximation is more appropriate.On the other hand, the crossover can occur for quite large system sizes, e.g., N ∼ 10 5 for τ = 0.1 and N ∼ 10 4 for τ = 1, meaning that demographic noise tends to persist if the reaction-rate noise is fast, but disappears quickly with N if the noise is slow; notice that the top series with τ = 10 has almost no N -dependence.

VI. CONCLUSIONS
Temporal fluctuations in the parameters that control contagion dynamics are inevitable, and have been shown in many epidemic data analyses.Motivated by this, we analyzed the effects of fluctuating infection and recovery rates on the outbreak-size distribution in the canonical SIR model.The SIR reaction rates were modeled with Ornstein Uhlenbeck noise, allowing us to extract the outbreak statistics as a function of the noise standard deviations and correlation times.Our simple choice was demonstrated by performing a model inference of the 2019-2020 RSV season in the US, where we observed significant temporal fluctuations in infectious contact rates.
In terms of analytical results, we found solutions for the outbreak-size distribution in the adiabatic and white noise regimes, and showed that the distributions can be highly skewed with significant probabilities for large fluctuations away from mean-field predictions.Interestingly, we discovered that the outbreak variance is generally maximized for a value of the basic reproductive number that depends on the correlation time of the noise, which in the white-noise limit is independent of where noise resides (infection or recovery).In addition, we compared the typical fluctuations emerging from demographic and reaction-rate noise and determined the population sizes, correlation times, and reproductive numbers that places noisy SIR systems in the various limiting regimes.Altogether then, our work illustrated a rich interplay between noise and outbreak dynamics -depending sensitively on fundamental noise characteristics and population size.
Currently the theory presented pertains to well-mixed populations in which individuals come into contact with a contagion with homogeneous rates.In actuality the contact rates within a population can be highly heterogeneous and/or spatially distributed, and therefore, an important extension of our work will include the generalization of the outbreak-size distribution to complex network topology.Another common assumption that we used, which is only an approximation, was the implicit exponential waiting times for both the infection and recovery processes.Future work should include generalization to gamma and other more realistic waiting-time distributions.Finally, our work has relied substantially on small-noise assumptions allowing us to focus on the dominant, exponential contribution to the outbreak-size distribution.Corrections to this approach, which would include next-order contributions for larger noise amplitudes, are an important avenue for future analysis. VII.

25 FIG. 1 .
FIG. 1. RSV model inference.(a) Weekly RSV hospitalizations (black dots) and 2-week rolling average (black line) from the 2019-2020 season in the United States[42].Results from the Bayesian inference model are overlaid with the data (median: red line, shaded bands represent the inter-quartile range and the 95% credible intervals).(b) Inferred infectious contact rate obeying a time-discretized version of the OU process (median: black line, shaded bands represent the inter-quartile range and the 95% credible intervals).

FIG. 5 .
FIG. 5. Variance of the final outbreak size versus the inverse of the correlation time (left) and N (right).Panel (a): The variance (rescaled by the predicted adiabatic-limit) as function of τ −1 for σ β = 0.04β0 and β0 = 2.0, 1.5 and 1.2 (top to bottom).The lines are the white-noise predictions.The inset shows the crossover time, Eq.(18) versus R0.Panel (b): The variance versus N −1 for τ = 10, 1 and 0.1 (top to bottom) with β0 = 2.0 and σ β = 0.1β0.The lines are the white-noise predictions, which are the sum of the variances from reactionrate and demographic noise.For both panels γ0 = 1.