# Outbreak-size distributions under fluctuating rates

Jason Hindes1, Luis Mier-y-Teran-Romero 2, Ira B. Schwartz1 and Michael Assaf
3 1U.S. Naval Research Laboratory, Washington, DC 20375, USA  2Leidos, Reston,
VA 20190, USA  3Racah Institute of Physics, Hebrew University of Jerusalem,
Jerusalem 91904, Israel

###### Abstract

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, 15, 4]. 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, 16, 12, 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, 15, 4, 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 [17, 1, 21, 22, 23, 24, 19, 25]. For
instance, recently we developed a theoretical approach that allows for
calculating the distribution of outbreak sizes in well-mixed 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, best-fit 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[34, 35, 36, 37, 38, 1, 31, 32], which in the
case of human epidemics correlate with more general social contact rates[39].
For instance, techniques for extracting time-dependent parameters have been
applied to the recent COVID-19 pandemic as well[40, 41, 31, 33], 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 timescales111If re-
infection happens on short timescales, then models such as SIRS may be more
appropriate[1]. 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 $\beta SI\\!/\\!N$ (for a well-mixed population), where
$\beta$ is the infectious contact rate[1, 2, 4]. Similarly, the probability
per unit time that the number of infected decreases by one is $\gamma I$,
where $\gamma$ 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:

$\dot{x}_{s}=-\beta x_{s}x_{i},\quad\dot{x}_{i}=\beta x_{s}x_{i}-\gamma
x_{i},\quad\dot{x}_{r}=\gamma x_{i},$ (1)

where $x_{s}+x_{i}+x_{r}=1$. Starting from a small initial infection density,
$x_{i}(t\\!=\\!0)\\!\ll\\!1$, the final fraction of susceptibles in Eqs.(1)
$\overline{x_{s}^{*}}\equiv x_{0}$ satisfies $x_{0}=e^{-R_{0}(1-x_{0})}$,
where $R_{0}=\beta/\gamma$ 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 $\overline{x_{r}^{*}}=1-x_{0}$,

$\overline{x_{r}^{*}}=1-e^{-R_{0}\overline{x_{r}^{*}}}.$ (2)

Note that if $R_{0}\leq 1$ in Eq.(2) then $\overline{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 $\tau$, while the noise variances
are $\sigma_{\beta}^{2}$ and $\sigma_{\gamma}^{2}$ for the infection and
recovery rates, respectively. Thus, we write:
$\beta(t)=\beta_{0}(1+\xi_{\beta}(t))$ and
$\gamma(t)=\gamma_{0}(1+\xi_{\gamma}(t))$, and augment Eq.(1) into the
stochastic system

$\displaystyle\dot{x}_{s}$ $\displaystyle=$
$\displaystyle-\beta_{0}(1+\xi_{\beta})x_{s}x_{i},$ $\displaystyle\dot{x}_{i}$
$\displaystyle=$
$\displaystyle\beta_{0}(1+\xi_{\beta})x_{s}x_{i}-\gamma_{0}(1+\xi_{\gamma})x_{i},$
$\displaystyle\dot{\xi}_{\beta}$ $\displaystyle=$
$\displaystyle-\frac{\xi_{\beta}}{\tau}+\sqrt{\frac{2}{\tau}}\Big{(}\frac{\sigma_{\beta}}{\beta_{0}}\Big{)}\eta_{\beta}(t),$
$\displaystyle\dot{\xi}_{\gamma}$ $\displaystyle=$
$\displaystyle-\frac{\xi_{\gamma}}{\tau}+\sqrt{\frac{2}{\tau}}\Big{(}\frac{\sigma_{\gamma}}{\gamma_{0}}\Big{)}\eta_{\gamma}(t).$
(3)

Here, $\eta_{\beta}$ and $\eta_{\gamma}$ are Gaussian white noises, while
$\xi_{\beta}$ and $\xi_{\gamma}$ are OU processes. Note that by construction,
$\beta(t)$ and $\gamma(t)$ are assumed to be wide-sense stationary Gaussian
processes with $\left<\beta\right>\\!=\\!\beta_{0}$,
$\left<\gamma\right>\\!=\\!\gamma_{0}$,
$\left<(\beta(t)-\beta_{0})(\beta(t+\Delta)-\beta_{0})\right>\\!=\\!\sigma_{\beta}^{2}e^{-\Delta/\tau}$,
and
$\left<(\gamma(t)-\gamma_{0})(\gamma(t+\Delta)-\gamma_{0})\right>\\!=\\!\sigma_{\gamma}^{2}e^{-\Delta/\tau}$,
where $\left<\cdot\right>$ denotes the expectation operator. In general, one
can simulate the system of equations (II) and find the final outbreak-size
distribution for fluctuating rates with any magnitude and correlation time.

Figure 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).

###  II.1  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[26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37, 38, 1, 31, 32, 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.(II) with daily time steps and a fixed recovery
rate $\gamma\\!=\\!\gamma_{0}\\!=\\!1/7\;\text{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: $\beta_{0}$, the inverse correlation time $\alpha$,
$\sigma_{\beta}$, the hospitalization rate, and initial conditions for the
SIR222The Bayesian framework estimates the initial state for the SIR model
that is used for prediction in a single year. The state includes the fraction
of the population recovered from previous years. Since, the model does not
include birth, death, re-infection, etc. it is not appropriate for multiple-
year predictions[49], but is reasonable for a single season.; 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: $\hat{R}_{0}\\!=\\!1.37$ in $[1.32,1.44]$,
$\hat{\alpha}\\!=\\!0.11\;\text{days}^{-1}$ in
$[0.045,0.20]\;\text{days}^{-1}$, and
$\hat{\sigma}_{\beta}\\!=\\!0.026\;\text{days}^{-1}$ in
$[0.014,0.040]\;\text{days}^{-1}$, where $\hat{\;\;}$ denotes the median
within a quartile range specified by the square brackets. From these we
observe a fairly tight value of the inferred time-averaged $R_{0}$, but with
substantial temporal fluctuations between $10-20\%$. On the other hand, the
correlation time estimate $\alpha^{-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.(II).

##  III  LIMIT OF ADIABATIC NOISE

In order to gain insight on the outbreak distributions generated from the
general Eqs.(II) and temporal fluctuations of the sort we inferred from RSV
data, we first consider limiting regimes. We start with the limit of adiabatic
noise, $\tau\gg 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
$\beta$ varies and $\gamma$ is constant, such that $\sigma_{\gamma}\\!=\\!0$.
In order to simplify the equations, we take $\gamma_{0}\\!=\\!1$, which merely
specifies the time units and results in $R_{0}\\!=\\!\beta$.

To find the distribution of the final outbreak size $P(x_{r}^{*})$, we have to
compute the following integral:

$P(x_{r}^{*})=\int_{-\infty}^{\infty}P(x_{r}^{*}|\beta)P(\beta)d\beta.$ (4)

The conditional probability $P(x_{r}^{*}|\beta)$ is a dirac delta function
around the mean-field value of the outbreak at $\beta$. Namely,
$P(x_{r}^{*}|\beta)=\delta(x_{r}^{*}-\overline{x_{r}^{*}})$, where the mean-
field final outbreak fraction $\overline{x_{r}^{*}}$ satisfies Eq.(2). Taking
a Gaussian distribution for $P(\beta)$ with mean $\beta_{0}$ and standard
deviation $\sigma_{\beta}$, Eq. (4) becomes

$P(x_{r}^{*})=\frac{1}{\sqrt{2\pi\sigma_{\beta}^{2}}}\int_{1}^{\infty}\delta(x_{r}^{*}-\overline{x_{r}^{*}}(\beta))e^{-\frac{(\beta-\beta_{0})^{2}}{2\sigma_{\beta}^{2}}}d\beta.$
(5)

We point out that in order for the SIR model to make physical sense
$\beta\\!\geq\\!0$. Therefore, when plugging in an unrestricted Gaussian in
Eq.(5), $\sigma_{\beta}$ cannot be too large 333we must assume that e.g.,
$\sigma_{\beta}/\beta_{0}\lesssim 1/3$ so that over $99\%$ of the distribution
is in the physical domain.. Otherwise other distributions, e.g. that vanish at
$\beta=0$ can be used instead; yet this does not change the results
qualitatively. We also note that since $\overline{x_{r}^{*}}(\beta)$ vanishes
for $\beta\leq 1$, the lower boundary in the integral in Eq. (5) can be taken
to be $1$ without affecting the distribution.

Changing variables from $\beta$ to $x_{r}^{*}$, and using the fact that
$d\beta/dx_{r}^{*}=[x_{r}^{*}+(1-x_{r}^{*})\ln(1-x_{r}^{*})]/[(1-x_{r}^{*})(x_{r}^{*})^{2}]$,
we can explicitly perform the integration by plugging instead of $\beta$,
$-\log(1-x_{r}^{*})/x_{r}^{*}$, which is the solution of
$x_{r}^{*}=\overline{x_{r}^{*}}(\beta)$. As a result, Eq. (5) reduces to

$P(x_{r}^{*})=\frac{d\beta}{dx_{r}^{*}}\exp\\!\left[-\frac{\left((1/x_{r}^{*})\ln(1-x_{r}^{*})+\beta_{0}\right)^{2}}{2\sigma_{\beta}^{2}}\right].$
(6)

From Eq.(6) we can derive, e.g., the typical fluctuations around the mean-
field given by $\sigma_{\text{a}}$ – the standard deviation associated with
the adiabatic outbreak PDF. In particular, in the limit of small
$\sigma_{\beta}$, the variance becomes

$\sigma_{\text{a}}\\!=\\!\sigma_{\beta}\\!\left.\\!\frac{dx_{r}^{*}(\beta)}{d\beta}\\!\right|_{\beta=\beta_{0}}\\!=\frac{\sigma_{\beta}x_{0}(1-x_{0})}{1-R_{0}x_{0}}.$
(7)

Note that for adiabatic noise we can repeat our calculation of the outbreak-
size distribution and variance for any quenched distribution of $\beta$ (or
$\gamma$), and not just Gaussians, which may better reflect the data in a
given application.

Figure 2: Outbreak statistics for adiabatic noise. (a) Simulated PDFs of the
final outbreak fraction $x_{r}^{*}$ from Eqs.(II), in the case of time-
correlated reaction-rate noise for $\tau=10^{0}$-$10^{3}$ (from narrowest to
widest) with $\sigma_{\beta}=0.1\beta_{0}$ and $\beta_{0}\\!=\\!2$. The solid
and dashed lines denote the adiabatic and white-noise predictions,
respectively. (b) Variance of the outbreak PDF (normalized by infection noise)
versus $R_{0}\\!=\\!\beta_{0}/\gamma_{0}$ for $\tau=10^{3}$-$10^{5}$ (bottom
to top) and $\sigma_{\beta}=0.04\beta_{0}$. The variance of the PDF (6) is
shown in solid black, while the small-noise limit Eq.(7) is shown with a
dashed line. For both panels $\gamma_{0}\\!=\\!1$.

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
$\beta_{0}\\!=\\!2$ and $\sigma_{\beta}=0.1\beta_{0}$. The solution from
Eq.(6) can be compared to stochastic simulations of Eqs.(II) for large $\tau$.
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\\!\sim\\!\exp[-S/\sigma_{\beta}^{2}]$, where
$S\\!=\\!(\ln(1-x_{r}^{*})/x_{r}^{*}+R_{0})^{2}/2$. When
$x_{r}^{*}\\!\rightarrow\\!0$ the action remains finite, i.e,
$S\\!\rightarrow\\!(R_{0}-1)^{2}/2$. On the other hand, when
$x_{r}^{*}\\!\rightarrow\\!1$, $S\\!\rightarrow\\!\infty$. 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 $\tau$
with $\gamma\\!=\\!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 $\sigma_{\beta}\to 0$, the maximum approaches
$R_{0}\\!=\\!1$. On the other hand as $\sigma_{\beta}$ increases, the maximum
variance occurs for an $R_{0}$ that is an increasing function of
$\sigma_{\beta}$. 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}\\!\rightarrow\\!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
$\tau$ 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,
$\tau\to 0$, instead of Eqs. (II) we can write

$\displaystyle\dot{x}_{s}$ $\displaystyle=$
$\displaystyle-\beta_{0}(1+\sigma_{1}\zeta_{\beta}(t))x_{s}x_{i},$
$\displaystyle\dot{x}_{i}$ $\displaystyle=$
$\displaystyle\beta_{0}(1+\sigma_{1}\zeta_{\beta}(t))x_{s}x_{i}-\gamma_{0}(1+\sigma_{2}\zeta_{\gamma}(t))x_{i},$
(8)

which we denote as the white reaction-rate limit. Here, $\zeta_{\beta}$ and
$\zeta_{\gamma}$ are white Gaussian noises. In order to coincide with Eqs.
(II) as $\tau\\!\to\\!0$, one must demand that

$\displaystyle\sigma_{1}\\!=\\!\sqrt{2\tau}\sigma_{\beta}/\beta_{0},\quad\sigma_{2}\\!=\\!\sqrt{2\tau}\sigma_{\gamma}/\gamma_{0}.$
(9)

To analyze the outbreak-size PDF given Eq. (IV), 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):

$\displaystyle\dfrac{\partial P}{\partial t}=$
$\displaystyle-\dfrac{\partial}{\partial
x_{s}}\Big{[}\\!\\!-\\!\beta_{0}x_{s}x_{i}P\Big{]}-\dfrac{\partial}{\partial
x_{i}}\Big{[}(\beta_{0}x_{s}x_{i}-\gamma_{0}x_{i})P\Big{]}\;+$ (10)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\Big{(}\dfrac{\partial^{2}}{\partial
x_{s}^{2}}+\dfrac{\partial^{2}}{\partial x_{i}^{2}}-2\dfrac{\partial}{\partial
x_{s}}\dfrac{\partial}{\partial
x_{i}}\Big{)}\Big{[}\frac{1}{2}\beta_{0}^{2}x_{s}^{2}x_{i}^{2}\sigma_{1}^{2}P\Big{]}\;+$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\dfrac{\partial^{2}}{\partial
x_{i}^{2}}\Big{[}\frac{1}{2}\gamma_{0}x_{i}^{2}\sigma_{2}^{2}P\Big{]}.$

To simplify notation, henceforth, we will assume that
$\sigma_{2}^{2}\\!=\\!f\sigma_{1}^{2}$, with $f\\!>\\!0$ for simplicity, and
again rescale time $t\to\gamma_{0}t$, so that $\beta_{0}$ is replaced by the
basic reproduction number, $R_{0}=\beta_{0}/\gamma_{0}$. Next, we employ the
WKB approximation
$P(x_{s},x_{i})\\!\sim\\!\exp\\!\left[-S(x_{s},x_{i})/\sigma_{1}^{2}\right]$,
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.(II). 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/\sigma_{1}^{2}$ 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,
$\partial{S}/\partial{t}\\!+\\!H\\!=\\!0$, in the leading order in
$\sigma_{1}\ll 1$, with

$H=x_{i}[p_{i}(R_{0}x_{s}-1)-R_{0}x_{s}p_{s}]+\frac{1}{2}x_{i}^{2}[R_{0}^{2}x_{s}^{2}(p_{s}-p_{i})^{2}+fp_{i}^{2}].$
(11)

In this formalism, $H$ is called the Hamiltonian, $S$ is the action, while
$p_{s}\\!=\\!\partial_{x_{s}}S$ and $p_{i}\\!=\\!\partial_{x_{i}}S$ are the
conjugate momenta, just as in analytical mechanics[51, 53, 54].

Figure 3: Scaling of the outbreak-size distribution. (a) The natural log of
simulated PDF values for $x_{r}^{*}\\!=\\!0.785$ (green stars), $0.780$ (blue
circles), $0.775$ (red squares), and $0.770$ (magenta diamonds) from Eqs.(II)
versus the reaction-rate noise variance. The slopes are predictions from the
white-noise theory. Other model parameters are $\beta_{0}\\!=\\!2$,
$\gamma_{0}\\!=\\!1$, $\tau\\!=\\!0.1$, and $f\\!=\\!0$. (b) The final
outbreak size plotted as a function of the infected momentum initial condition
multiplied by $R_{0}\\!-\\!1$ for: $R_{0}\\!=\\!1.1$ (solid black),
$R_{0}\\!=\\!1.r$ (dashed red), and $R_{0}\\!=\\!2$ (dotted blue). Other model
parameters are $\gamma_{0}\\!=\\!1$ and $f\\!=\\!0.65$.

In order to compute probabilities for different outbreak sizes, we need to
know the action $S$, which means calculating the integrals $S\\!=\\!\int
p_{s}dx_{s}+\int p_{i}dx_{i}-\int Hdt$ given the dynamics of Hamilton’s
equations: $\dot{x}_{s}\\!=\\!\partial_{p_{s}}H$,
$\dot{x}_{i}\\!=\\!\partial_{p_{i}}H$, $\dot{p}_{s}\\!=\\!-\partial_{x_{s}}H$,
and $\dot{p}_{i}\\!=\\!-\partial_{x_{i}}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)\\!\rightarrow\\!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
$H\\!=\\!p_{s}\dot{x}_{s}+p_{i}\dot{x}_{i}\\!-\\!(1/2)R_{0}^{2}x_{s}^{2}x_{i}^{2}\big{[}p_{s}-p_{i}\big{]}^{2}\\!-\\!(1/2)fx_{i}^{2}p_{i}^{2}$,
using $\dot{x}_{s}$ and $\dot{x}_{i}$. Third, by substituting the zero-energy
condition into $\dot{p}_{i}$, we obtain that
$\dot{p}_{i}=-(1/2)R_{0}^{2}x_{s}^{2}x_{i}\big{[}p_{s}-p_{i}\big{]}^{2}-(1/2)fx_{i}p_{i}^{2}$.
As a result, the Hamiltonian (11) simplifies to:
$H(t)=p_{s}\dot{x}_{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=\int\\!p_{s}dx_{s}+x_{i}(t\rightarrow\infty)p_{i}(t\rightarrow\infty)-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
$\int\\!p_{s}dx_{s}\\!=\\!0$. As a consequence, the action associated with an
outbreak in the white-noise limit is simply

$S=\int\\!p_{i}dx_{i}.$ (12)

###  IV.1  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\\!\rightarrow\\!\infty$. 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
$\dot{x}_{s}\\!=\\!0$, $\dot{x}_{i}\\!=\\!0$, $\dot{p}_{s}\\!=\\!0$, and
$\dot{p}_{i}\\!=\\!0$, given $x_{s}\\!=\\!1$ and $x_{i}\\!=\\!0$. It is
straightforward to show that the general fixed-point initial conditions are

$(x_{s},x_{i},p_{s},p_{i})_{{}_{\scaleto{t=0}{5pt}}}\;=\;(1,0,\delta(R_{0}-1)/R_{0},\delta),$
(13)

where $\delta\equiv 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 $\delta$.
Examples are shown in Fig.4 (b), where we plot the outbreak sizes as a
function of $\delta$ 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 $\delta$, (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 $\delta$. Each $\delta$ 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 $\delta$ 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
multi-component white noise, which we return to in Sec.V.

Figure 4: Outbreak statistics for white noise. (a) The simulated white-noise
PDFs from Eqs.(IV) for $\sigma_{1}=0.1$ and $f=0.1,\;1,\;2$ and $5$ (from
narrowest to widest) with $\beta_{0}\\!=\\!2$. White-noise predictions for
each combination are shown with curves overlaying the simulation results. (b)
Variance of the outbreak PDF (normalized by noise variance) for white
reaction-rate noise versus $R_{0}\\!=\\!\beta_{0}/\gamma_{0}$. The noise
combinations are: blue circles ($\sigma_{1}\\!=\\!0.1$, $f\\!=\\!0.1$), red
squares ($\sigma_{1}\\!=\\!0.1$, $f\\!=\\!1$), green diamonds
($\sigma_{1}\\!=\\!0.1$, $f\\!=\\!5$), magenta triangles
($\sigma_{1}\\!=\\!0.2$, $f\\!=\\!0.1$), and black hexagrams
($\sigma_{1}\\!=\\!0.04$, $f\\!=\\!1$). The black curve shows the white-noise
predictions. In both panels $\gamma_{0}\\!=\\!1$.

###  IV.2  Outbreak variance

In the general case of noise in both $\beta$ and $\gamma$, it seems that
Hamilton’s equations cannot be solved analytically in a simple manner – apart
from constructing a power-series expansion in the initial-condition parameter
$\delta$. 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 $\delta$. A complete solution for the case of recovery-only
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}\ll 1$, to be verified a-posteriori. Equating $H=0$, yields:
$p_{s}=(1-1/(R_{0}x_{s}))p_{i}$, i.e., $p_{s}\ll 1$ as well. Second, we show
that $p_{i}(t)$ remains small during the entire epidemic duration as long as
the initial momentum $\delta$ is small. Writing down the Hamilton’s equation
for $\dot{p}_{i}=-\partial_{x_{i}}H$, and using Eq. (11) we have:
$\dot{p}_{i}=-(1/2)(1+f)x_{i}(t)p_{i}^{2}$. The solution of this differential
equation is

$p_{i}(t)\simeq\delta\left\\{1+[(1+f)/(2R_{0})]\ln(x_{s})\delta\right\\},$
(14)

where $p_{i}(0)=\delta$ is the initial condition, and we have used the fact
that in the leading order in $\delta\ll 1$, $x_{r}=\int
x_{i}dt=\ln(x_{s})/R_{0}$. This is legitimate as the action will have a
$\delta^{2}$ dependence, see below, so we can substitute in ${\cal
O}(\delta^{2})$ terms their mean-field ${\cal O}(\delta^{0})$ approximation.

To compute $\delta$, we can use Hamilton’s equations for $\dot{x}_{s}$, and
$\dot{x}_{r}=-\dot{x}_{s}-\dot{x}_{i}$, and compute
$\dot{x}_{r}/\dot{x}_{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 ${\cal O}(\delta)$. Using Eq. (14) and that when the outbreak ends
$x_{r}^{*}=1-x_{s}^{*}$, and assuming $x_{s}^{*}-x_{0}\sim{\cal O}(\delta)$
(to be confirmed a-posteriori), we find

$\delta\simeq\frac{2(1-R_{0}x_{0})}{(1+f)(1-x_{0})x_{0}(2-R_{0}(1+x_{0}))}(x_{s}^{*}-x_{0}).$
(15)

This confirms a-posteriori that $\delta\ll 1$, under the assumption that the
final susceptible fraction is close to its mean-field counterpart, i.e.
$x_{s}^{*}-x_{0}\ll 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: $\int p_{i}dx_{i}=\int
p_{i}(dx_{i}/dx_{s})dx_{s}$. Here, the Jacobian can be found using the
Hamilton’s equations:
$dx_{i}/dx_{s}\\!=\\!1/(R_{0}x_{s})\\!-\\!1\\!-\\![(1\\!+\\!f)/(R_{0}x_{s})(1\\!-\\!x_{s}\\!+\\!\ln
x_{s}/R_{0})]\delta$, where again we have used mean-field results for the
${\cal O}(\delta)$ 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}\sim\delta$:

$\displaystyle S$ $\displaystyle=$
$\displaystyle\frac{(x_{s}^{*}-x_{0})^{2}}{2v^{2}}+\mathcal{O}\big{(}(x_{s}^{*}-x_{0})^{3}\big{)},$
$\displaystyle v^{2}$ $\displaystyle=$
$\displaystyle\frac{(1+f)(x_{0}-1)R_{0}x_{0}^{2}[2-R_{0}(1+x_{0})]}{2(1-R_{0}x_{0})^{2}}.$
(16)

Indeed, having obtained a $\delta^{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}^{*})\\!\sim\\!\exp\\!\left[-S(x_{s}^{*})/\sigma_{1}^{2}\right]$, the
variance of the outbreak-size distribution in the limit of white reaction-rate
noise, $\sigma_{w}^{2}$, is

$\displaystyle\sigma_{w}^{2}=\sigma_{1}^{2}v^{2}.$ (17)

We can test our predictions for the outbreak variance in the white-noise
regime by performing stochastic simulations of Eqs.(IV) 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}\simeq 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}^{*}\simeq 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,
$\sigma_{1}^{2}+\sigma_{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 $\mathcal{O}(\delta^{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\\!\rightarrow\\!\infty$ limit, and try to
understand how small (large) $\tau$ 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.(II) versus the
(inverse) correlation time $\tau$ for three values of $R_{0}$ and fixed
$\sigma_{\beta}$. Note that the outbreak variance is normalized by the
adiabatic limit, Eq.(7), so that each simulation series approaches unity for
small $\tau^{-1}$. In addition to the adiabatic limit, for comparison we plot
predictions for white-noise, Eq.(IV.2), with lines. In the latter case, the
$\tau$-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 $\beta$ 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
$\tau$. 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,
$\tau\\!\sim\\!\gamma_{0}^{-1}\\!=\\!1$.

In general, the crossover point in $\tau$ between white- and adiabatic-noise
regimes has some $R_{0}$ dependence: the smaller $R_{0}$, the larger $\tau$
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,
$\tau_{c}$, can be found by solving $\sigma_{a}^{2}\\!=\\!\sigma_{w}^{2}$, or

$\displaystyle\tau_{c}=\dfrac{R_{0}^{2}}{2v^{2}}\dfrac{\sigma_{a}^{2}}{\sigma_{\beta}^{2}}=\frac{R_{0}(x_{0}-1)}{(1+f)(2-R_{0}(1+x_{0})},$
(18)

where we have used Eqs. (7) and (IV.2), valid for small noise. Evidently,
$\tau_{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\sim\mathcal{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.1 into Eq.(18), and using Eqs. (7) and (IV.2), we find that the ratio
of the noise correlation time to the cross-over time, $\tau/\tau_{c}$, falls
between $0.1\\!-\\!0.2$. As the $\tau$ 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.

Figure 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 $\tau^{-1}$ for
$\sigma_{\beta}=0.04\beta_{0}$ and $\beta_{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 $R_{0}$. Panel (b): The variance versus
$N^{-1}$ for $\tau=10$, $1$ and $0.1$ (top to bottom) with $\beta_{0}=2.0$ and
$\sigma_{\beta}=0.1\beta_{0}$. The lines are the white-noise predictions,
which are the sum of the variances from reaction-rate and demographic noise.
For both panels $\gamma_{0}\\!=\\!1$.

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.(II). In
Fig.5(b) we plot the outbreak-size variance as a function of $N^{-1}$ for
three values of the correlation time $\tau=10,\;1,$ and $0.1$ (from top to
bottom). The curves are the expected total white-noise variance,
$\sigma^{2}_{w,tot}$, which is a sum of reaction-rate and demographic noise

$\displaystyle\sigma^{2}_{w,tot}=\dfrac{(\sigma_{1}^{2}+\sigma_{2}^{2})(x_{0}-1)R_{0}x_{0}^{2}[2-R_{0}(1+x_{0})]}{2(1-R_{0}x_{0})^{2}}\;\;+$
$\displaystyle\dfrac{x_{0}(1-x_{0})(R_{0}^{2}x_{0}+1)}{N(R_{0}x_{0}-1)^{2}}.$
(19)

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
$\sigma_{1}^{2},\sigma_{2}^{2},N^{-1}\\!\ll\\!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.(IV.2) – especially for the two smaller values of $\tau$
where the white-noise approximation is more appropriate. On the other hand,
the crossover can occur for quite large system sizes, e.g.,
$N\\!\sim\\!10^{5}$ for $\tau\\!=\\!0.1$ and $N\\!\sim\\!10^{4}$ for
$\tau\\!=\\!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 $\tau\\!=\\!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  ACKNOWLEDGMENTS

JH and IBS were supported by the U.S. Naval Research Laboratory funding
(N0001419WX00055), and the Office of Naval Research (N0001419WX01166) and
(N0001419WX01322). MA was supported by the Israel Science Foundation Grant No.
531/20.

##  VIII  APPENDICES

##  Appendix A  RSV data analysis

We implement a time-discrete version of an SIR model to carry out Bayesian
parameter inference on the hospitalization data for the 2019-2020 season of
respiratory syncytial virus (RSV) in the US[42]. At day $d$, the daily model
updates follow the dynamics:

$\displaystyle i_{d}$ $\displaystyle=\beta_{d}\cdot S_{d-1}\cdot I_{d-1},$
(20a) $\displaystyle S_{d}$ $\displaystyle=S_{d-1}-i_{d-1},$ (20b)
$\displaystyle I_{d}$ $\displaystyle=I_{d-1}+i_{d}-\gamma\cdot I_{d-1}$ (20c)
$\displaystyle R_{d}$ $\displaystyle=R_{d-1}+\gamma\cdot I_{d-1}$ (20d)

where $(S_{d},I_{d},R_{d})$ are the susceptible, infected and recovered
fractions ($S_{d}+I_{d}+R_{d}=1$), respectively, and $i_{d}$ is the fractional
incidence. The parameters $\beta_{d}$ and $\gamma$ refer to the contact and
recovery rates, respectively.

In our inference model, the dynamics of $\beta_{d}$ are given by a daily-
discretized version of the Ornstein-Uhlenbeck process:

$\displaystyle\beta_{1}$ $\displaystyle\\!=\\!$
$\displaystyle\mathcal{N}(\beta_{0},\sigma_{\beta}),$ (21)
$\displaystyle\beta_{d}$ $\displaystyle\\!=\\!$
$\displaystyle\beta_{d-1}\\!-\\!\alpha(\beta_{d-1}\\!-\\!\beta_{0})\\!+\\!\mathcal{N}(0,2\alpha\sigma_{\beta}),\;\;\;d\\!>\\!1\\!,$
(22)

where $\mathcal{N}(\mu,\sigma)$ is a normal random variable with mean $\mu$
and variance $\sigma^{2}$.

The parameter inference process is done by tying the daily-discretized
statistical model to the number of daily hospitalizations contained in the
data $H_{d}$ by the Poisson observation process

$\displaystyle H_{d}=\text{Pois}(\eta\cdot N\cdot i_{d})$ (23)

where $\eta$ is the hospitalization rate and $N$ denotes the total population.
With this model we do Bayesian parameter inference using the platform Stan via
the R package rstan[45, 46]; R code for the model inference is available upon
request.

##  Appendix B  Action for white-noise fluctuations in recovery

In addition to small fluctuations in the outbreaks, we can gain further
analytical insight into the outbreak distribution for white noise by looking
at other limiting-case scenarios. One such scenario is when the recovery-noise
dominates over infection noise, namely $f\\!\gg\\!1$ while
$\sigma_{\gamma}\\!\ll\\!1$. In this limit, the Hamiltonian reduces to

$H=x_{i}[p_{i}(R_{0}x_{s}-1)-R_{0}x_{s}p_{s}]+\frac{1}{2}fx_{i}^{2}p_{i}^{2}.$
(24)

As in Sec.IV.2, given the simpler Hamiltonian (24), the $H\\!=\\!0$ condition
can be combined with Hamilton’s equations $\dot{x}_{s}$ and $\dot{p}_{i}$, to
give $p_{i}$ as an explicit function of $x_{s}$ and the initial condition for
$\delta\equiv p_{i}(t\\!=\\!0)$:

$p_{i}(x_{s})=\frac{2\delta}{2-f\delta\ln(x_{s})/R_{0}}.$ (25)

Similarly, if we substitute Eq. (25) and the zero-energy condition into
$\dot{x}_{i}$ and divide by $\dot{x}_{s}$, we get the following differential
equation for $x_{i}$ as a function of $x_{s}$:

$\frac{dx_{i}}{dx_{s}}=-1+\frac{1}{R_{0}x_{s}}-\frac{2fx_{i}\delta}{R_{0}x_{s}\Big{(}2-f\delta\ln(x_{s})/R_{0}\Big{)}}.$
(26)

Next, we can solve for $x_{i}(x_{s})$ by separating the fraction of the
population infected into a product of two functions that depend on $x_{s}$,
i.e., $x_{i}(x_{s})\\!=\\!u(x_{s})v(x_{s})$. By substituting the product form
into Eq.(26), setting $u\;dv/dx_{s}\\!=\\!-1+1/(R_{0}x_{s})$, and conditioning
$u(x_{s}\\!=\\!1)\\!=\\!1$ and $v(x_{s}\\!=\\!1)\\!=\\!0$, we find:

$u(x_{s})\\!=\\!\frac{1}{4}\Bigg{(}\\!\\!\frac{f\delta}{R_{0}}\ln(x_{s})-2\\!\Bigg{)}^{\\!\\!2}\\!\\!,\;\;v(x_{s})\\!=\\!\int_{1}^{x_{s}}\\!\\!\\!\dfrac{1-R_{0}x_{s}^{\prime}}{R_{0}x_{s}^{\prime}u(x_{s}^{\prime})}\;dx_{s^{\prime}},$
(27)

where the integral in Eq.(27) can be expressed in terms of incomplete
exponential integrals, though the formula is cumbersome. As
$x_{i}(t\\!\rightarrow\\!\infty)\\!\rightarrow\\!0$, Eq. (27) gives us a
condition for the final outbreak $1-x_{s}^{*}$. Namely, given $R_{0}$ and
$\delta$, we can solve $v(x_{s}^{*})=0$ for the unique value of $x_{s}^{*}$.

The final step for calculating the action associated with a given outbreak in
the limit of recovery-only fluctuations, is to differentiate
$x_{i}(x_{s})\\!=\\!u(x_{s})v(x_{s})$ and substitute Eqs.(25) and (27) into
Eq.(12). The result is the following limiting-case action

$S=\bigintss_{1}^{x_{s}^{*}}\\!\\!\dfrac{2\delta\Big{(}-1+\dfrac{1}{R_{0}x_{s}}\Big{)}}{\Big{(}f\delta\ln(x_{s})/{R_{0}}-2\Big{)}}\;dx_{s},$
(28)

where as mentioned, the boundary condition for the integral (the final
outbreak size) can be determined numerically from the condition
$v(x_{s}^{*})=0$.

## References

  * Keeling and Rohani [2008] M. J. Keeling and P. Rohani,  _Modeling Infectious Diseases in Humans and Animals_ (Princeton University Press, 2008).
  * Andersson and Britton [2000] H. Andersson and T. Britton,  _Stochastic epidemic models and their statistical analysis_ (Springer-Verlag, New York, 2000).
  * Pastor-Satorras  _et al._ [2015] R. Pastor-Satorras, C. Castellano, P. Van; Mieghem, and A. Vespignani, Rev. Mod. Phys. 87, 925 (2015).
  * Rodrigues [2016] H. S. Rodrigues, “Application of sir epidemiological model: new trends,” (2016), arXiv:1611.02565 [physics.soc-ph] .
  * Billings  _et al._ [2002] L. Billings, W. M. Spears, and I. B. Schwartz, Physics Letters A 297, 261 (2002).
  * Dorogovtsev  _et al._ [2008] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Rev. Mod. Phys. 80, 1275 (2008).
  * Hindes and Assaf [2019] J. Hindes and M. Assaf, Physical review letters 123, 068301 (2019).
  * Volkening  _et al._ [2020] A. Volkening, D. F. Linder, M. A. Porter, and G. A. Rempala, SIAM Review 62, 837 (2020), https://doi.org/10.1137/19M1306658 .
  * Reiner and et al. [2021] R. C. Reiner and et al., Nature Medicine 27, 94 (2021).
  * Ghaemi  _et al._ [2022] M. Ghaemi, M. Shojafar, Z. Zabihinpour, and Y. Asgari, PLOS ONE 17, e0265065 (2022).
  * Rodriguez-Barraquer  _et al._ [2014] I. Rodriguez-Barraquer, L. M. y Teran-Romero, I. B. Schwartz, D. S. Burke, and D. A. Cummings, Vaccine 32, 514 (2014).
  * Catching  _et al._ [2021] A. Catching, S. Capponi, M. T. Yeh, S. Bianco, and R. Andino, medRxiv (2021), 10.1101/2020.08.12.20173047.
  * Hindes  _et al._ [2022] J. Hindes, M. Assaf, and I. B. Schwartz, Physical Review Letters 128, 078301 (2022).
  * Hindes  _et al._ [2021] J. Hindes, S. Bianco, and I. B. Schwartz, PLOS ONE 16, 1 (2021).
  * Brauer [2017] F. Brauer, Infectious Disease Modelling 2, 113 (2017).
  * Ray and et al. [2020] E. L. Ray and et al., medRxiv (2020), 10.1101/2020.08.19.20177493.
  * Ganyani  _et al._ [2021] T. Ganyani, C. Faes, and N. Hens, Annual Review of Statistics and Its Application 8, 69 (2021), https://doi.org/10.1146/annurev-statistics-061120-034438 .
  * Batistela  _et al._ [2021] C. M. Batistela, D. P. F. Correa, Á. M. Bueno, and J. C. Piqueira, Chaos, solitons, and fractals 142, 110388 (2021).
  * Miller [2019] J. C. Miller, “Distribution of outbreak sizes for sir disease in finite populations,” (2019), arXiv:1907.05138 [q-bio.PE] .
  * Aron and Schwartz [1984] J. L. Aron and I. B. Schwartz, Journal of Theoretical Biology 110, 665 (1984).
  * Allen [2017] L. J. S. Allen, Infectious Disease Modelling 2, 128 (2017).
  * House  _et al._ [2013] T. House, J. V. Ross, and D. Sirl, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 469, 20120436 (2013).
  * Ball [1986] F. Ball, Advances in Applied Probability 18, 289 (1986).
  * Ball and Clancy [1993] F. Ball and D. Clancy, Advances in Applied Probability 25, 721 (1993).
  * Dykman  _et al._ [1994] M. I. Dykman, E. Mori, J. Ross, and P. M. Hunt, J. Chem. Phys. 100, 5735 (1994).
  * Hong and Li [2020] H. G. Hong and Y. Li, PLOS ONE 15, 1 (2020).
  * Kahn  _et al._ [2020] R. Kahn, C. M. Peak, J. Fernández-Gracia, A. Hill, A. Jambai, L. Ganda, M. C. Castro, and C. O. Buckee, Proceedings of the National Academy of Science 117, 5067 (2020).
  * Vossler  _et al._ [2022] H. Vossler, P. Akilimali, Y. Pan, W. R. KhudaBukhsh, E. Kenah, and G. A. Rempała, Scientific Reports 12, 5534 (2022).
  * Dharmaratne  _et al._ [2020] S. Dharmaratne, S. Sudaraka, I. Abeyagunawardena, K. Manchanayake, M. Kothalawala, and W. Gunathunga, Virology Journal 17, 144 (2020).
  * Bsat  _et al._ [2022] R. Bsat, H. Chemaitelly, P. Coyle, P. Tang, M. R. Hasan, Z. Al Kanaani, E. Al Kuwari, A. A. Butt, A. Jeremijenko, A. H. Kaleeckal, A. N. Latif, R. M. Shaik, G. K. Nasrallah, F. M. Benslimane, H. A. Al Khatib, H. M. Yassine, M. G. Al Kuwari, H. E. Al Romaihi, M. H. Al-Thani, A. Al Khal, R. Bertollini, L. J. Abu-Raddad, and H. H. Ayoub, J Glob Health 12, 05004 (2022).
  * Spannaus  _et al._ [2022] A. Spannaus, T. Papamarkou, S. Erwin, and J. B. Christian, Scientific Reports 12, 10761 (2022).
  * Dureau  _et al._ [2013] J. Dureau, K. Kalogeropoulos, and M. Baguelin, Biostatistics 14, 541 (2013), https://academic.oup.com/biostatistics/article-pdf/14/3/541/17739168/kxs052.pdf .
  * Kucharski  _et al._ [2020] A. J. Kucharski, T. W. Russell, C. Diamond, Y. Liu, J. Edmunds, S. Funk, R. M. Eggo, F. Sun, M. Jit, J. D. Munday, N. Davies, A. Gimma, K. van Zandvoort, H. Gibbs, J. Hellewell, C. I. Jarvis, S. Clifford, B. J. Quilty, N. I. Bosse, S. Abbott, P. Klepac, and S. Flasche,  _The Lancet Infectious Diseases_ , The Lancet Infectious Diseases 20, 553 (2020).
  * London and Yorke [1973] W. London and J. Yorke, American Journal of Epidemiology 98, 453 (1973).
  * Yorke  _et al._ [1979] J. Yorke, N. Nathanson, G. Pianigiani, and J. Martin, American Journal of Epidemiology 109, 103 (1979).
  * Fine and Clarkson [1980] P. Fine and J. Clarkson, International Journal of Epidemiology 11, 5 (1980).
  * Grenfell  _et al._ [2002] B. Grenfell, O. Bjørnstad, and B. Finkenstädt, Ecological Monographs 72, 185 (2002).
  * Schwartz [1992] I. B. Schwartz, Journal of Mathematical Biology (1992).
  * Drolet  _et al._ [2022] M. Drolet, A. Godbout, M. Mondor, G. Béraud, L. Drolet-Roy, P. Lemieux-Mellouki, A. Bureau, É. Demers, M.-C. Boily, C. Sauvageau, G. De Serres, N. Hens, P. Beutels, B. Dervaux, and M. Brisson, BMC Public Health 22, 1032 (2022).
  * Chen  _et al._ [2020] Y.-C. Chen, P.-E. Lu, C.-S. Chang, and T.-H. Liu, IEEE Transactions on Network Science and Engineering 7, 3279 (2020).
  * Setianto and Hidayat [2023] S. Setianto and D. Hidayat, Scientific Reports 13 (2023), 10.1038/s41598-023-31714-5.
  * [42] “Respiratory syncytial virus hospitalization surveillance network,” .
  * Note [1] If re-infection happens on short timescales, then models such as SIRS may be more appropriate[1].
  * Rose  _et al._ [2018] E. Rose, A. Wheatley, G. Langley, S. Gerber, and A. Haynes, MMWR Morb Mortal Wkly Rep 67, 71 (2018).
  * Stan Development Team [2023] Stan Development Team, “RStan: the R interface to Stan,”  (2023), r package version 2.21.8.
  * Hawryluk  _et al._ [2020] I. Hawryluk, T. A. Mellan, H. Hoeltgebaum, S. Mishra, R. P. Schnekenberg, C. Whittaker, H. Zhu, A. Gandy, C. A. Donnelly, S. Flaxman, and S. Bhatt, J R Soc Interface 17, 20200596 (2020).
  * Note [2] The Bayesian framework estimates the initial state for the SIR model that is used for prediction in a single year. The state includes the fraction of the population recovered from previous years. Since, the model does not include birth, death, re-infection, etc. it is not appropriate for multiple-year predictions[49], but is reasonable for a single season.
  * [48] Code available upon request.
  * White  _et al._ [2007] L. J. White, J. N. Mandl, M. G. M. Gomes, A. T. Bodley-Tickell, P. A. Cane, P. Perez-Brena, J. C. Aguilar, M. M. Siqueira, S. A. Portes, S. M. Straliotto, M. Waris, D. J. Nokes, and G. F. Medley, Math Biosci 209, 222 (2007).
  * Note [3] We must assume that e.g., $\sigma_{\beta}/\beta_{0}\lesssim 1/3$ so that over $99\%$ of the distribution is in the physical domain.
  * Dykman [1990] M. Dykman, Physical Review A 42, 2020 (1990).
  * Freidlin  _et al._ [2012] M. Freidlin, J. Szücs, and A. Wentzell,  _Random Perturbations of Dynamical Systems_, Grundlehren der mathematischen Wissenschaften (Springer, 2012).
  * Assaf and Meerson [2017] M. Assaf and B. Meerson, Journal of Physics A: Mathematical and Theoretical 50, 263001 (2017).
  * Forgoston and Moore [2018] E. Forgoston and R. O. Moore, SIAM Review 60, 969 (2018), https://doi.org/10.1137/17M1142028 .
  * Assaf  _et al._ [2013] M. Assaf, M. Mobilia, and E. Roberts, Phys. Rev. Lett. 111, 238101 (2013).
  * Taitelbaum  _et al._ [2020] A. Taitelbaum, R. West, M. Assaf, and M. Mobilia, Phys. Rev. Lett. 125, 048105 (2020).
  * Thattai and van Oudenaarden [2001] M. Thattai and A. van Oudenaarden, Proc Natl Acad Sci U S A 98, 8614 (2001).


# Degrees in link graphs of regular graphs

Itai Benjamini and John Haslegrave

###### Abstract

We analyse an extremal question on the degrees of the link graphs of a finite
regular graph, that is, the subgraphs induced by non-trivial spheres. We show
that if $G$ is $d$-regular and connected but not complete then some link graph
of $G$ has minimum degree at most $\lfloor 2d/3\rfloor-1$, and if $G$ is
sufficiently large in terms of $d$ then some link graph has minimum degree at
most $\lfloor d/2\rfloor-1$; both bounds are best possible. We also give the
corresponding best-possible result for the corresponding problem where
subgraphs induced by balls, rather than spheres, are considered.

We motivate these questions by posing a conjecture concerning expansion of
link graphs in large bounded-degree graphs, together with a heuristic
justification thereof.

##  1 Link graphs with large degrees

For a graph $G$, radius $r>0$ and vertex $v$, the sphere of radius $r$ about
$v$, which we denote $\mathcal{S}_{r}(v)$, is the set of vertices at distance
exactly $r$ from $v$, i.e. $\mathcal{S}_{r}(v):=\\{w\in V(G):d(v,w)=r\\}$. If
$\mathcal{S}_{r}(v)\neq\varnothing$, the $r$-link graph of $v$, which we
denote $\mathcal{L}_{r}(v)$, is the induced subgraph $G[\mathcal{S}_{r}(v)]$.
The graph $\mathcal{L}_{1}(e)$, where $G$ is a Cayley graph of a finitely-
generated group, was used in [8] to give a sufficient condition for property
(T).

Suppose that $G$ is a $d$-regular graph other than $K_{d+1}$. Can we choose
$G$ in such a way as to ensure that each of the link graphs has large minimum
degree? In other words, we wish to find the maximum value of
$\min_{v,r}\delta(\mathcal{L}_{r}(v))$ over all non-complete $d$-regular
graphs.

######  Theorem 1.

If $G$ is a non-complete connected $d$-regular graph then

$\min_{v,r}\delta(\mathcal{L}_{r}(v))\leq\lfloor 2d/3\rfloor-1,$

and this bound can be attained for any $d\equiv 2\pmod{3}$.

###### Proof.

We verify the second statement first. Fix $d=3k-1$, so that $\lfloor
2d/3\rfloor-1=2k-2$, and let $G$ be the graph obtained by blowing up each
vertex of $C_{5}$ to a clique of order $k$. Then $G$ is $(3k-1)$-regular and
vertex-transitive with diameter $2$. For each $v$, the link graph
$\mathcal{L}_{1}(v)$ consists of three cliques of order $k$, $k-1$ and $k$,
with all vertices between the clique of order $k-1$ and the other two cliques,
so has minimum degree $2k-2$. Furthermore, $\mathcal{L}_{2}(v)\cong K_{2k}$
and has minimum degree $2k-1$.

Next we prove the upper bound. Let $G$ be any non-complete $d$-regular graph,
and note that this implies $G$ has no universal vertex. Suppose that every
link graph of $G$ has minimum degree at least $m$. For any vertex $v$, any
vertex in $\mathcal{S}_{1}(v)$ has at least $m$ neighbours within
$\mathcal{S}_{1}(v)$ (and is a neighbour of $v$), so has at most $d-m-1$
neighbours in $\mathcal{S}_{2}(v)$. Choose $x\in\mathcal{S}_{2}(v)$ and
$w\in\mathcal{S}_{1}(v)\cap\mathcal{S}_{1}(x)$. Now we must have
$\deg_{\mathcal{L}_{1}(x)}(w)\geq m$. Since $w$ has at most $d-m-1$ neighbours
in $\mathcal{S}_{2}(v)$, one of which is $x$, it has at most $d-m-2$
neighbours in $\mathcal{S}_{2}(v)\cap\mathcal{S}_{1}(x)$. Consequently, $w$
has at least $m-(d-m-2)=2m+2-d$ neighbours in
$\mathcal{S}_{1}(v)\cap\mathcal{S}_{1}(x)$, and it follows that $x$ has at
least $2m+3-d$ neighbours in $\mathcal{S}_{1}(v)$. However, considering
$\deg_{\mathcal{L}_{2}(v)}(x)$, it also has at least $m$ neighbours in
$\mathcal{S}_{2}(v)$. Since $\deg_{G}(x)=d$, we must have $3m+3-d\leq d$, i.e.
$m\leq 2d/3-1$ and $m$ is an integer. ∎

Note that the example constructed to show that the bound of Theorem 1 is tight
is small relative to $d$, and has diameter $2$. We can do better if we exclude
such examples.

######  Theorem 2.

For each $r\geq 3$, if $G$ is a connected $d$-regular graph with diameter at
least $r$ then

$\min_{v,r}\delta(\mathcal{L}_{r}(v))\leq\lfloor rd/(2r-1)\rfloor-1.$

###### Proof.

Suppose $G$ is such a graph and every link graph has minimum degree at least
$m$. Fix a vertex $v$ such that $\mathcal{S}_{r}(v)\neq\varnothing$. We claim
that for each $2\leq j\leq r$ every vertex in $\mathcal{S}_{j}(v)$ has at
least $(2j-2)m-(j-1)d+2j-1$ neighbours in $\mathcal{S}_{j-1}(v)$. The result
will follow from the case $j=r$ of the claim, since any vertex in
$\mathcal{S}_{r}(v)$ also has at least $m$ neighbours in $\mathcal{S}_{r}(v)$,
and so $(2j-1)m-(j-1)d+2j-1\leq d$, giving $m\leq jd/(2j-1)-1$.

We prove the claim by induction. The case $j=2$ was shown in the proof of
Theorem 1. For $j>2$, let $x$ be a vertex in $\mathcal{S}_{j}(v)$ and $w$ be a
vertex in $\mathcal{S}_{j-1}(v)\cap\mathcal{S}_{1}(x)$. Since $w$ has at least
$(2j-4)m-(j-2)d+2j-3$ neighbours in $\mathcal{S}_{j-2}(v)$ by the induction
hypothesis, and at least $m$ neighbours in $\mathcal{S}_{j-1}(v)$, it has at
most $(j-1)d-(2j-3)m-2j+2$ neighbours in $\mathcal{S}_{j}(v)\setminus\\{x\\}$.
Since $w$ has at least $m$ neighbours in $\mathcal{S}_{1}(x)$, it must have at
least $(2j-2)m-(j-1)d+2j-2$ neighbours in
$\mathcal{S}_{1}(x)\cap\mathcal{S}_{j-1}(v)$, so
$\lvert\mathcal{S}_{1}(x)\cap\mathcal{S}_{j-1}(v)\rvert\geq(2j-2)m-(j-1)d+2j-1$,
as required. ∎

From Theorem 2 it follows that if $\operatorname{diam}(G)$ is sufficiently
large in terms of $d$ (in particular, is at least $(d+1)/2$), and consequently
whenever $\lvert G\rvert$ is sufficiently large in terms of $d$, then

$\min_{v,r}\delta(\mathcal{L}_{r}(v))\leq\lfloor d/2\rfloor-1.$

In fact this bound is best possible.

######  Proposition 3.

For each even $d$ there exist connected $d$-regular graphs of arbitrarily
large diameter satisfying $\min_{v,r}\delta(\mathcal{L}_{r}(v))=d/2-1$.

###### Proof.

Set $k=d/2$ and choose any $n$ satisfying $n-1\equiv a\pmod{2k}$ for some
$a\in\\{k,\ldots,2k\\}$. Consider the $k$th power of the $n$-cycle,
$C_{n}^{k}$. This is $2k$-regular, and has diameter
$s:=\lfloor(n-1)/2k\rfloor$, which may be made arbitrarily large. Every link
graph of radius strictly smaller than $s$ consists of two cliques of order $k$
with possibly some edges in between, so has minimum degree at least $k-1$. A
sphere of radius $s$ is a set of $b:=n-1-2k(s-1)$ consecutive vertices; since
this value is congruent to $n-1$ modulo $2k$, positive and at most $2k$, we
have $b=a\geq k$. Consequently $\mathcal{L}_{s}(v)$ will have minimum degree
$k-1$. ∎

In this example, the diameter grows linearly with the order of the graph. It
is natural to ask about large graphs which are well-connected, having diameter
growing logarithmically with their order. In this case it is still possible
for all link degrees to be linear in $d$. For example, start from a large
cubic expander and blow up each vertex to a clique of order $(d+1)/4$, where
$d$ is fixed. The graph obtained is $d$-regular, and retains good expansion
and logarithmic diameter. Each link graph consists of some cliques, possibly
with edges between them, and has minimum degree at least $(d-3)/4$. However,
we do not know whether the constant $1/4$ in this example can be improved.

Additionally, in the example of Proposition 3, almost all links are
$(d/2-1)$-regular. We might ask whether this is a necessary feature; in
particular, do there exist $d$-regular graphs of large diameter for which all
links have average degree greater than $d/2$? Note that in general there exist
examples where all links have average degree significantly higher than the
minimum degree over all links. For example, take the Cartesian product of a
triangle and a long odd cycle, and blow up all vertices to cliques of order
$(d+1)/5$. The average degree of each link is at least $(3d-7)/10$, but most
links have minimum degree $(d-4)/5$. However, we do not know of a similar
example where the minimum degree is close to $d/2$.

##  2 Induced subgraphs on balls

In this section we consider a natural extension, replacing link graphs with
the subgraphs induced by balls, i.e. we consider induced subgraphs of the form
$\overline{\mathcal{L}}_{r}(v):=G[\mathcal{B}_{r}(v)]$ where
$\mathcal{B}_{r}(v):=\\{w\in V(G):d(v,w)\leq r\\}$. In this case we can easily
have all minimum degrees of such graphs close to $d$ by taking $G$ to have
only slightly more than $d+1$ vertices; for example, if $d$ is even then
taking $G$ to be $K_{d+2}$ minus a perfect matching ensures
$\min_{v,r}\delta(\overline{\mathcal{L}}_{r}(v))=d-2$.

However, again such graphs have diameter $2$, and if we require greater
diameter we obtain non-trivial (and in fact tight) bounds.

######  Theorem 4.

If $G$ is a connected $d$-regular graph with diameter at least $3$ then

$\min_{v,r}\delta(\overline{\mathcal{L}}_{r}(v))\leq\lfloor(2d-1)/3\rfloor,$

and this bound can be attained for any $d\equiv 2\pmod{3}$ and arbitrarily
large diameter.

###### Proof.

Fix $d=3k-1$, so that $\lfloor 2d/3\rfloor-1=2k-2$, choose $n\geq 6$
arbitrarily, and let $G$ be the graph obtained by blowing up each vertex of
$C_{n}$ to a clique of order $k$. Then $G$ is $(3k-1)$-regular and vertex-
transitive with diameter $\lfloor(n-1)/2\rfloor\geq 3$. For each $v$ and $r$,
the graph $\overline{\mathcal{L}}_{r}(v)$ consists of the cliques
corresponding to $\min\\{2r-1,n\\}$ consecutive vertices of the cycle, and so
has minimum degree at least $2k-1=(2d-1)/3$.

Suppose $G$ is $d$-regular with diameter at least $3$ and satisfies
$\min_{v,r}\delta(\overline{\mathcal{L}}_{r}(v))\leq\lfloor 2d/3\rfloor-1$.
Choose vertices $v,y$ with $d_{G}(v,y)=3$, and let $vwxy$ be a shortest path
between them. Since $w$ has at least $m$ neighbours in $\mathcal{B}_{1}(v)$,
it has at most $d-m$ neighbours in $\mathcal{S}_{2}(v)$. Since $w$ also has at
least $m$ neighbours in $\mathcal{B}_{1}(x)$, it has at least $2m-d$
neighbours in $\mathcal{B}_{1}(x)\setminus\mathcal{S}_{2}(v)$, which must all
be in $\mathcal{S}_{1}(v)$. It follows that $x$ has at least $2m-d+1$
neighbours in $\mathcal{S}_{1}(v)$. However, $x$ also has at least $m$
neighbours in $\mathcal{B}_{1}(y)$, which is disjoint from
$\mathcal{S}_{1}(v)$. Thus $3m-d+1\leq\deg(x)\leq d$, and so $m\leq(2d-1)/3$.
∎

##  3 A conjecture

The results in Section 1 are in spirit indicating that spheres are not “too
connected”. Below is a conjecture regarding $r$-links.

A finite graph $G$ is called an $h$-expander for

$h=\inf_{S\subset V_{G}:\,\,0<\lvert S\rvert\leq\lvert
G\rvert/2}\frac{\lvert\partial S\rvert}{\lvert S\rvert},$

where $V_{G}$ are the vertices of $G$ and $\partial S$ is the outer vertex
boundary of $S$.

An expander family is a sequence of graphs such that, for some $h>0$, all
graphs in the sequence have expansion at least $h$.

######  Conjecture 1.

There is no sequence of bounded degree finite graphs, with size growing to
infinity, such that all links in all the graphs form an expander family.

Note that size growing to infinity is equivalent to diameter growing to
infinity. In the extremal examples discussed in Section 1 (other than the
example with diameter $2$), the link graphs do not expand, since they are
disconnected (and for the example with logarithmic diameter, typically have
many components).

In what follows, we provide some heuristic support for Conjecture 1.

######  Lemma 5.

Assume $G$ is an infinite graph, and all the $r$-links in the graph are
uniform expanders. Then $G$ is non-amenable.

###### Proof.

Fix a vertex, first note that the size of spheres around the vertex grows to
infinity with $r$, as otherwise if we have infinitely bounded cutsets
separating the vertex from infinity. Take a vertex which is in the middle
between two remote such cuts, say at distance $R_{0}$ from the cuts. The
sphere centred at that vertex of radius $R_{0}$ will be disconnected, and so
not an expander.

Let $S$ be a finite set of vertices in $G$. Pick a vertex $v$ in $G$ far
enough from $S$, such that the $r$-links around $v$ that intersect $S$ has
size larger than $2^{\lvert S\rvert}$. By uniform expansion of the $r$-links,
in each $r$-link that intersects $S$, the intersection of $S$ with the
$r$-link has proportional boundary. and as the radius of the spheres varied
they give a disjoint cover of $S$. Thus every finite set of vertices in $G$
has proportional boundary and thus $G$ is non-amenable.∎

The heuristic for Conjecture 1 follows the strategy of Salez [7]. By Lemma 5,
the BS limit of graphs in which all links are expanders with expansion bounded
away from zero is a.s. non-amenable. A simple random walk on a non-amenable
graph has positive linear speed. The BS limit is a unimodular random graph
(see [1, 5] for definitions) in which all the $r$-links of the root are an
expander family. In [4] it was proved that if you add edges so that the
vertices of each of the levels of a binary tree are uniform expanders, then
the resulting graph has no non-constant bounded harmonic functions. We
_believe_ , (but cannot prove) that the same holds for unimodular random
graphs. Since it is equivalent to zero speed [3], this will lead to a
contradiction.

Our result in Section 2 equally relates to a similar question. In [2] it was
asked: Is there is a sequence of finite bounded-degrees graphs growing in size
to infinity, so that all the induced balls in all the graphs in the sequence
form an expander family? For a related results on heat kernel supports see
[6].

## Acknowledgements

JH was supported by the UK Research and Innovation Future Leaders Fellowship
MR/S016325/1.

## References

  * [1] David Aldous and Russell Lyons,  Processes on unimodular random networks.  Electron. J. Probab. 12 (2007), no. 54, 1454–1508.
  * [2] Itai Benjamini,  Expanders are not hyperbolic.  Israel J. Math. 108 (1998), 33–36.
  * [3] Itai Benjamini and Nicolas Curien,  Ergodic theory on stationary random graphs.  Electron. J. Probab. 17 (2012), no. 93, 20 pp.
  * [4] Itai Benjamini and Gady Kozma,  Nonamenable Liouville Graphs.  Preprint, arxiv:1010.3365 (2010).
  * [5] Itai Benjamini and Oded Schramm,  Recurrence of distributional limits of finite planar graphs.  Electron. J. Probab. 6 (2001), no. 23, 13 pp.
  * [6] Mikolaj Fraczyk and Wouter van Limbeek,  Heat kernels are not uniform expanders.  Preprint, arXiv:1905.13584 (2019).
  * [7] Justin Salez,  Sparse expanders have negative curvature.  Preprint, arXiv:2101.08242 (2021).
  * [8] Andrzej Żuk,  Property (T) and Kazhdan constants for discrete groups.  Geom. Funct. Anal. 13 (2003), no. 3, 643–670.