††thanks: JZ and FPS contributed equally††thanks: JZ and FPS contributed
equally

# How contact patterns destabilize and modulate epidemic outbreaks

Johannes Zierenberg  F. Paul Spitzner  Max Planck Institute for Dynamics and
Self-Organization, 37077 Göttingen, Germany  Viola Priesemann  Max Planck
Institute for Dynamics and Self-Organization, 37077 Göttingen, Germany
Institute for the Dynamics of Complex Systems, University of Göttingen, 37077
Göttingen, Germany Martin Weigel  Institut für Physik, Technische Universität
Chemnitz, 09107 Chemnitz, Germany  Michael Wilczek  Max Planck Institute for
Dynamics and Self-Organization, 37077 Göttingen, Germany  Institute for the
Dynamics of Complex Systems, University of Göttingen, 37077 Göttingen, Germany
Theoretical Physics I, University of Bayreuth, 95440 Bayreuth, Germany

###### Abstract

The spread of a contagious disease clearly depends on when infected
individuals come in contact with others. Yet, it remains unclear how the
timing of contacts, in particular the time of day and day of week, interacts
with the latent and infectious stages of the disease. Here, we use real-world
physical proximity data to study this interaction, and find that, compared to
randomized controls, human contact patterns destabilize epidemic outbreaks and
non-trivially modulate the basic reproduction number. By exploring simple
generative models constrained by those data, we are able to attribute both of
these intriguing observations to distinct aspects of the temporal statistics
of contact patterns. We find the destabilization of outbreaks to be caused by
a high probability of extreme events, i.e., super-spreading and zero-
spreading, which we are able to reproduce in our models by including temporal
clustering of contacts. Furthermore, contact patterns can either increase or
decrease the basic reproduction number depending on the latent period. We find
this modulation to be caused by the alignment of the infectious period with
reoccurring daily and weekly variations in the conditional rate of secondary
contacts, which we reproduce by including a cyclostationary contact rate in
our models. Thus, by identifying and reproducing relevant non-Markovian
statistics of human contacts, our work opens the possibility to systematically
study the non-equilibrium physics of realistic disease spread, as well as of
non-Markovian spreading processes in general.

epidemic spread $|$ human contact behavior $|$ non-Markovian dynamics $|$
latent period

##  I Introduction

As contagious diseases are passed on through contacts, the number of secondary
infections depends crucially on the contact patterns of infectious
individuals. These contact patterns encode relevant information such as the
number of interaction partners and contact timing. However, the majority of
prevailing models for disease spread prioritize simpler descriptions that
neglect these aspects. Although studies that do explicitly model contact
patterns get more involved, they also generally find that the effects of
contact patterns are crucial for the understanding of disease spread: On the
one hand, _interaction partners_ can be modeled by a static complex network
[1], whose structure can affect epidemic spread through the occurrence of hubs
[2, 3], multiscale link communities [4] and influential spreaders [5]. On the
other hand, real-world _interaction times_ generally follow a non-Markovian
process (in contrast to commonly assumed memoryless processes), which can
influence epidemics through the occurrence of bursts [6, 7] and daily and
weekly variations in human interaction [8, 9].

Aiming for a complete description of time-varying interactions leads to
_temporal networks_ with time-varying links [10, 11]. However, constructing
detailed temporal networks from real-world contacts requires extensive amounts
of recorded data, which can be gleaned from field studies [12, 13, 14] but
that are notoriously limited in system size and time duration. The limitations
of real-world data can be partly remedied by generating surrogate data [15,
16], but it is often unclear to which extent they represent the real system.
It is thus necessary to identify those features of real-world human contact
behavior that impact epidemic outbreaks, complementing past insights about the
relevance of inter-individual variability from studies of static networks [1].

In order to reveal temporal statistics of human contact patterns that affect
non-equilibrium epidemic outbreaks, we here combine an analysis of real-world
contact data with a novel approach of extracting disease-spread statistics.
Instead of characterizing full epidemic outbreaks on a large temporal network
(likely under-determined by the data at hand), we investigate the statistics
of potential secondary infections (the expectation value of which is
proportional to the basic reproduction number $R_{0}$). Our analysis reveals
two mechanisms through which non-Poissonian statistics of contact patterns
affect epidemic outbreaks: (i) temporal correlations of contacts in the form
of clusters destabilize outbreaks due to an increased probability of zero-
spreading and super-spreading events, and (ii) temporal variation of the
contact rate modulates the pace of epidemic spread due to the interference of
contact patterns with disease progression. As a proof of principle for a new
class of generative models that incorporate human contact patterns, we
showcase a novel, non-Markovian process that faithfully reproduces the
temporal statistics of real-world contact data and their effect on disease
spread.

##  II Results

Figure 1:  Real-world contacts represented as encounter trains. a: A contact
between two individuals is defined as successive co-location times. Focusing
on contagious disease transmission, we filter contacts to be closer than 2
meters and longer than 15 minutes. b: The start times of remaining contacts
per individual form their encounter train. A raster plot of such encounter
trains show a clear temporal structure of human contact patterns. c,d:
Randomization preserves the number of encounters per train and destroys
temporal structure.

To understand how the temporal structure of human contact patterns determines
epidemic outbreaks, we consider the statistics of potential secondary
infections. We focus on real-world contact patterns from physical proximity
data among university students collected in the Copenhagen Networks Study [17]
(Fig. 1), complementing our results with additional data from SocioPatterns
[18]. We remove data that are presumably irrelevant for disease spread and
hence only include contacts that are closer than
$2\text{\,}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{s}$ and
longer than
$15\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{u}\mathrm{t}\mathrm{e}\mathrm{s}$,
in line with common guidelines used by regulators in the context of the
COVID-19 pandemic (cf. SI for different distances and durations). From the
remaining contact events, we extract for each individual a list of their
contact starting times (encounters). These encounter trains are a point-
process-like representation that captures the non-Poissonian statistics of
real-world contact patterns.

To demonstrate the relevance of the non-Poissonian statistics of contact
patterns, we compare to randomized encounter trains with Poisson-distributed
encounter times. In this surrogate data, we draw for each encounter a new time
uniformly within the duration of the experiment (at the temporal resolution).
We thereby preserve the number of encounters per train (Fig. 1c) — and the
inter-individual variability — but destroy any temporal memory. This becomes
directly apparent in the distribution of daily encounters (Fig. 1d), which
becomes a superposition of Poisson distributions after randomizing, and we
show that this affects the statistics of potential secondary infections.

To estimate the statistics of potential secondary infections, we assume a
contagious disease that can be transmitted only at infectious encounters. If
$\tau$ is the time elapsed since an infection, the disease can be either in a
latent phase $\tau\in[0,T_{\mathrm{lat}})$, during which individuals are
infected but not infectious, or in an infectious phase
$\tau\in[T_{\mathrm{lat}},T_{\mathrm{lat}}+T_{\mathrm{inf}})$, during which
the disease can be transmitted with probability $p_{\mathrm{inf}}$ at an
encounter. Thus, a measure for potential secondary infections is the number of
encounters within the infectious period, which we call potentially infectious
encounters $n_{\mathrm{inf}}$ (Fig. 2a). We estimate the statistics of
$n_{\mathrm{inf}}$ by sampling hypothetical infections that may start at any
encounter of the data set. For the case of fixed $T_{\mathrm{lat}}$ and
$T_{\mathrm{inf}}$, this can be solved by a simple enumeration.

### Human contact patterns destabilize epidemic outbreaks

Figure 2:  Real-world contact patterns increase the probability of rare
spreading events, which increase the extinction probability of an epidemic
outbreak. a: Potentially infectious encounters, $n_{\mathrm{inf}}$, are
estimated from encounter trains for a disease model with latent period
$T_{\mathrm{lat}}$ and infectious period $T_{\mathrm{inf}}$. b:
$P(n_{\mathrm{inf}})$ depends on $T_{\mathrm{lat}}$ for human encounter trains
(red, blue) but not for randomized encounter trains (yellow lines), which
underestimate both zero-spreading as well as super-spreading events. c: Data-
driven branching model where $n_{\mathrm{inf}}$ are drawn from the empirical
distribution $P(n_{\mathrm{inf}})$ and get infected with constant probability
$p_{\mathrm{inf}}$. d: The extinction probability of an epidemic outbreak is
increased for human contact patterns compared to randomized, which can be
explained by the increased probability of small $n_{\mathrm{inf}}$ in (b).

We start by comparing the distributions $P(n_{\mathrm{inf}})$ obtained for two
examples of disease progression with the same $T_{\mathrm{inf}}$ but different
$T_{\mathrm{lat}}$ (Fig. 2b). We find that human contact patterns lead to high
probabilities of very small $n_{\mathrm{inf}}$ (zero-spreading events) and
very large $n_{\mathrm{inf}}$ (super-spreading events) compared to randomized.
An increased probability of low offspring number is known to destabilize
epidemic outbreaks: Especially during the first generations of an outbreak,
case numbers are still low and fluctuations could cause zero offspring,
stopping the outbreak despite a reproductive number $R_{0}>1$ [19].

One way to quantify the destabilization of an outbreak is through the
extinction probability, i.e., the fraction of outbreaks that asymptotically
end up in the absorbing state of zero infections. To quantify how the
extinction probability is affected by the empirical distribution
$P(n_{\mathrm{inf}})$, we construct a data-driven branching process (Fig. 2c).
Here, each infection first generates $n_{\mathrm{inf}}$ encounters according
to $P(n_{\mathrm{inf}})$, each of which is infected independently with
probability $p_{\mathrm{inf}}$ resulting in binomial-distributed secondary
infections. By summing over all possible $n_{\mathrm{inf}}$, this yields the
offspring distribution

$P(x)=\sum_{n_{\mathrm{inf}}=x}^{\infty}P(n_{\mathrm{inf}})\binom{n_{\mathrm{inf}}}{x}p_{\mathrm{inf}}^{x}(1-p_{\mathrm{inf}})^{n_{\mathrm{inf}}-x}.$
(1)

For branching processes, one can calculate the asymptotic extinction
probability with the help of the probability generating function,
$\pi(\theta)=\sum_{x=0}^{\infty}p(x)\theta^{x}$: The extinction probability is
the smallest $\theta^{\ast}$ that solves $\theta=\pi(\theta)$ [19], which we
solve numerically because $P(n_{\mathrm{inf}})$ is an empirical distribution.

For our chosen examples of $T_{\mathrm{lat}}$, we find that human contact
patterns increase the extinction probability of outbreaks compared to
randomized patterns (Fig. 2d). To compare between the examples, we plot the
extinction probability as a function of the basic reproduction number

$R_{0}=p_{\mathrm{inf}}\,\overline{n}_{\mathrm{inf}},$ (2)

where $\overline{n}_{\mathrm{inf}}$ is the average number of potentially
infectious encounters. (Note that this representation accounts for the
surprising differences in $\overline{n}_{\mathrm{inf}}$ between examples, cf.
Fig. 2b, vertical dashed lines).

As expected for branching processes, the extinction probability changes at
$R_{0}=1$: For $R_{0}\leq 1$, $P_{\mathrm{ext}}=1$ such that asymptotically no
outbreak occurs, and for $R_{0}>1$, $P_{\mathrm{ext}}<1$ such that there is a
non-zero chance for an outbreak. In addition, we observe that using
$P(n_{\mathrm{inf}})$ from human encounter trains leads to a higher extinction
probability than using $P(n_{\mathrm{inf}})$ from randomized encounter trains.
Whereas the extinction probability for randomized trains is independent of
$T_{\mathrm{lat}}$, we find for human encounter trains that $T_{\mathrm{lat}}$
determines the extinction probability by changing the probability of low
$n_{\mathrm{inf}}$ (Fig. 2b).

Together, we find that human contact patterns with their non-Poissonian
statistics increase the probability of zero-spreading and super-spreading
events compared to randomized encounter trains, and, thereby, destabilize
epidemic outbreaks.

### Interplay between contact pattern and disease progression modulates basic
reproduction number

Figure 3:  Real-world contact patterns modulate the pace of epidemic spread as
a function of latent period. a: Absolute $\overline{n}_{\mathrm{inf}}$ are a
periodic function of $T_{\mathrm{lat}}$ for human encounter trains but not for
randomized (yellow). b: Relative
$\overline{n}_{\mathrm{inf}}$/$\overline{n}_{\mathrm{inf}}^{\mathrm{rand}}$
reveals periodic modulations on both daily and weekly scale in the full
$(T_{\mathrm{lat}},T_{\mathrm{inf}})$-plane. c: $\Psi(\tau)$ features daily
and weekly modulations for non-Poissonian human encounter trains but is
constant for randomized encounter trains (yellow). Assuming the initial
encounter to be an infection, this explains the modulations of
$\overline{n}_{\mathrm{inf}}$ by combinations of $T_{\mathrm{lat}}$ and
$T_{\mathrm{inf}}$ for which the integral (shaded areas) is dominated by
valleys (blue) or hills (red). d: Continuous-time branching model, where
encounter times are generated from $\Psi(\tau)$ and infected with probability
$p_{\mathrm{inf}}$ for $R_{0}$ average secondary infections. e: Initializing
the model with $I_{0}=100$ random infections in
$[-T_{\mathrm{lat}}-T_{\mathrm{inf}},0)$, the (average) number of infections
grows exponentially. The pace of epidemic spread decreases with
$T_{\mathrm{lat}}$ but is modulated for non-Poissonian $\Psi(\tau)$. f: The
growth rate $\lambda$ for time-independent encounter times in a fixed
$T_{\mathrm{inf}}$ is expected to decrease trivially with $T_{\mathrm{lat}}$.
If not constant, $\Psi(\tau)$ modulates $\lambda$ and causes regimes of
slower-than-random (blue) or faster-than-random (red) pace of epidemic spread.

As noted above, we observe that $\overline{n}_{\mathrm{inf}}$ depends on
$T_{\mathrm{lat}}$ for human encounter trains. This is at first glance
surprising, because for memory-less Poisson processes
$\overline{n}_{\mathrm{inf}}$ is proportional to $T_{\mathrm{inf}}$ but
independent of $T_{\mathrm{lat}}$. Hence, in the following we systematically
vary $T_{\mathrm{lat}}$ and $T_{\mathrm{inf}}$ to study how the interplay
between human contact patterns and disease progression affects
$\overline{n}_{\mathrm{inf}}$, which via Eq. (2) also affects $R_{0}$.

For a fixed
$T_{\mathrm{inf}}=$3\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$, we
observe a systematic dependence of $\overline{n}_{\mathrm{inf}}$ on
$T_{\mathrm{lat}}$ for human encounter trains that is not present for
randomized encounter trains (Fig. 3a). For human encounter trains
$\overline{n}_{\mathrm{inf}}$ can be higher ($T_{\mathrm{lat}}<1$ and
$4<T_{\mathrm{lat}}<8$) or lower ($1<T_{\mathrm{lat}}<4$) as compared to the
randomized case, which suggests a periodic modulation.

To further study the dependence on $T_{\mathrm{inf}}$, we consider
$\overline{n}_{\mathrm{inf}}$ relative to randomized (Fig. 3b). This accounts
for the trivial increase of $\overline{n}_{\mathrm{inf}}$ with increasing
$T_{\mathrm{inf}}$. For small
$T_{\mathrm{inf}}<$1\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}$$, we find daily
modulations as a function of $T_{\mathrm{lat}}$, with regions below-randomized
(blue) and above-randomized (red). This effect diminishes with increasing
$T_{\mathrm{inf}}$ but remains as daily sub-modulation. For larger
$T_{\mathrm{inf}}$, we find extended triangular regions that are below-
randomized (blue) or above-randomized (red). The interfaces between these
regions are located at $T_{\rm lat}+T_{\rm
inf}=$7\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$ and at
$T_{\mathrm{lat}}=$7\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$, which
suggests that these are weekly modulations. Summarizing, we find periodic
modulations of $\overline{n}_{\mathrm{inf}}$ as a function of
$T_{\mathrm{lat}}$, with periods on the order of days for small
$T_{\mathrm{inf}}$ and on the order of weeks (with daily sub-modulations) for
larger $T_{\mathrm{inf}}$.

In order to understand how changes in $T_{\mathrm{lat}}$ result in periodic
modulations of $\overline{n}_{\rm inf}$, we introduce the conditional
encounter rate $\Psi(\tau)$ to describe the average rate of encounters
conditioned on an initial encounter at time $0$ (Fig. 3c). Importantly, when
considering the initial encounter as an infection, $\Psi(\tau)$ measures the
rate of potentially infectious encounters, but neglects fluctuations across
hypothetical infections (cf. Fig. 2b). We find that $\Psi(\tau)$ features a
peak at $0\text{\,}\mathrm{(}$which implies strong clustering [6]) and has
periodic modulations between high and low encounter rates on the scale of days
and weeks, which imply statistically reoccurring contact patterns. Both, the
initial peak as well as periodic modulations are lost for randomized encounter
trains (yellow line).

If we interpret the initial encounter as an infection, then $\Psi(\tau)$ marks
periods of statistically high and low encounter rates after infection, and
thereby implies a time-dependent secondary attack rate. From $\Psi(\tau)$, we
can obtain an estimate of $\overline{n}_{\mathrm{inf}}$ for a particular
disease progression ($T_{\mathrm{lat}}$, $T_{\mathrm{inf}}$) by integrating
within the infectious period (Fig. 3c, shaded areas). For the two previous
examples with fixed
$T_{\mathrm{inf}}=$3\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$, we
reproduce a lower $\overline{n}_{\mathrm{inf}}$ for
$T_{\mathrm{lat}}=$2\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$ (blue
area covering the valley) and a larger $\overline{n}_{\mathrm{inf}}$ for
$T_{\mathrm{lat}}=$6\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$ (red
area covering the 7-day peak). We thus find that the alignment of the
infectious period with regions of low or high $\Psi(\tau)$ determines
$\overline{n}_{\mathrm{inf}}$.

Since $\overline{n}_{\mathrm{inf}}$ is related to the basic reproduction
number $R_{0}$ via Eq. (2), the interplay between contact patterns and disease
progression modulates the pace of epidemic spread. To illustrate this, we
construct a branching process that is continuous in time, where each exposed
individual generates encounters according to $\Psi(\tau)$. During the
infectious period, encounters again have a probability $p_{\mathrm{inf}}$ to
become infected (Fig. 3d). Assuming an outbreak that survived the initial
generations, we prepare the system with $100$ random initial infections in the
interval $[-T_{\mathrm{lat}}-T_{\mathrm{inf}},0]$. The time evolution of daily
new cases shows clear exponential growth (Fig. 3e), where the growth rate
$\lambda$ trivially decreases with the generation time and, thus,
$T_{\mathrm{lat}}$ (Fig. 3f). However, this expected decrease of $\lambda$ for
Poisson contacts (yellow) is modulated by $\Psi(\tau)$ for human contact
patterns, which results in slower-than-random (blue) or faster-than-random
(red) growth.

Summarizing, the dependence of $\overline{n}_{\mathrm{inf}}$ on
$T_{\mathrm{lat}}$ implies that human contact patterns modulate the basic
reproduction number and thereby the growth rate of an epidemic outbreak.

Figure 4:  Specific temporal statistics of surrogate point processes can be
related to specific characteristics of disease spread. top: Tailored renewal
process that captures time-dependent $\rho(t)$ (first column), heavy-tailed
$P(\delta t)$ (second column), and heterogeneous $P(n_{\mathrm{train}})$
(third column), reproduces core characteristics such as modulation of
$\overline{n}_{\mathrm{inf}}$ (fourth column, cf. Fig. 3 and Supplementary
Material for additional figures for $\Psi(\tau)$) and increased probability of
zero-spreading in $P(n_{\mathrm{inf}})$ (fifth column, cf. Fig. 2). center: An
inhomogeneous Poisson process with only time-varying, cyclostationary
$\rho(t)$ and heterogeneous $P(n_{\mathrm{train}})$ reproduces the modulations
in $\overline{n}_{\mathrm{inf}}$, but underestimates the probability of low
$n_{\mathrm{inf}}$ due to a lack of clustering and the corresponding large
inter-cluster times in $P(\delta t)$. bottom: A Weibull renewal process with
contact clustering from $P(\delta t)$ and heterogeneous
$P(n_{\mathrm{train}})$ reproduces the high probability of rare outcomes in
$P(n_{\mathrm{inf}})$, but cannot reproduce the modulations in
$\overline{n}_{\mathrm{inf}}$ due to a lack of cyclic temporal structure. If,
additionally, one relaxes the constraint of heterogeneous
$P(n_{\mathrm{train}})$ and considers statistically identical trains, then
both surrogate processes underestimate the probability of large
$n_{\mathrm{inf}}$ related to super-spreading (see Supplementary Material).

### Destabilization and modulation of epidemic spread can be attributed to
specific temporal statistics of contact patterns

After illustrating that non-Poissonian statistics of human contact patterns
can destabilize and modulate epidemic outbreaks, it seems natural to ask how
they can be included in models of disease spread. In such models it is common
to approximate encounter times between individuals as memory-less Poisson
processes [1]. Assuming independence, these Poisson processes can be merged to
result in encounter trains with Poisson statistics — the same statistics as
our randomized encounter trains. In the following, we construct encounter
trains with non-Poissonian statistics and identify three specific features of
contact patterns that are necessary to reproduce the relevant statistics of
encounters. As a proof of principle, we showcase a novel inhomogeneous
Weibull-renewal process that is constrained by data and reproduces all
features (Fig. 4, top row):

i) Focusing on temporal statistics, the encounter rate $\rho(t)$ averaged
across individuals and weeks is time-dependent but cyclostationary; $\rho(t)$
repeats in a weekly cycle with differences between day and night, and between
weekdays and weekends (Fig. 4, first column). This can already be captured by
an inhomogeneous-Poisson process (Fig. 4, middle row), which reproduces the
periodic modulation of $\overline{n}_{\mathrm{inf}}$ (and $\Psi(\tau)$, see
Supplemental Fig. S3).

ii) The distribution of inter-encounter intervals $P(\delta t)$ has high
probability for small $\delta t$ and a heavy tail of non-vanishing probability
for large $\delta t$ (second column). Because this tail corresponds to long
periods without any encounters, it causes the high probability of
$n_{\mathrm{inf}}\approx 0$ (last column) that is responsible for the
destabilization of epidemic outbreaks. $P(\delta t)$ is dominantly shaped by
the clustering of human contacts and can be well approximated by a Weibull
distribution [7, 20]. Accordingly, a Weibull-renewal processes (last row) well
reproduces $P(\delta t)$ and $P(n_{\mathrm{inf}})$, but it does not have a
time-varying $\rho(t)$ and cannot reproduce the period modulations of
$\Psi(\tau)$ and $\overline{n}_{\mathrm{inf}}$.

iii) Encounter rates vary between individuals (third column). Such variability
can be attributed to intrinsic differences in contact behavior (partly
captured by the degree distribution of the contact network [21]). We find that
this across-individual variability is crucial to reproduce the heavy tailed
$P(n_{\mathrm{inf}})$ (see Supplemental Fig. S3 for generative processes where
individuals share a common rate).

We conclude that epidemic outbreaks are destabilized due to the temporal
clustering of encounters, but are modulated due to a cyclostationary encounter
rate. Clearly, models of disease spread can benefit from a generative process
that reproduces those relevant features of human contact patterns, such as the
tailored Weibull-renewal process showcased here. However, although our process
reproduces all discussed features, it is built heuristically (Methods), and
future work is needed to construct microscopic models that give rise to
cyclostationary rates with clustering in a principled way, while remaining
mathematically tractable.

##  III Discussion

We analyzed real-world human contact patterns and found that their non-
Poissonian timings shape epidemic spread in two important ways. Firstly, they
decrease the probability of an epidemic outbreak from an initial infection
compared to random (Poisson) contact patterns, which we attribute to the
clustering of contacts. While clustering is typically associated with super-
spreading events, it inevitably causes periods of zero contacts that we find
to increase the probability of zero-spreading events. The resulting
destabilization despite super-spreading is consistent with previous results
derived from branching processes with non-Poissonian offspring distributions
due to individual reproduction numbers [22]. In addition to individual
variability, we here identify temporal clustering as a necessary ingredient to
reproduce empirical non-Poissonian offspring distributions. In the regime of
low case numbers, where zero-spreading events can stop an outbreak, human
contact patterns thus destabilize epidemic spread.

Secondly, human contact patterns non-trivially modulate the pace of epidemic
spread depending on the latent period, which we attribute to time-dependent
but cyclostationary encounter rates. A cyclostationary rate leads to periods
of statistically high and low encounter rates conditioned on a potential
infection. How these periods align with the infectious period is affected by
the latent period and determines whether the number of potential secondary
infections, and in turn $R_{0}$, increases or decreases. This modulation of
$R_{0}$ can thus be understood as a resonance following either a constructive
or destructive interference between a periodically changing contact rate and
the disease progression. This resonance is a new mechanism to explain
previously observed slow-down or speed-up of diffusion processes on temporal
networks due to non-Markovian characteristics [23].

In the main manuscript we focus on deterministic disease progression with
fixed periods ($T_{\mathrm{lat}}$, $T_{\mathrm{inf}}$), but we also consider
non-deterministic disease progression with gamma-distributed periods [24, 25,
26] in the Supplementary Material (Fig. S1). We find our main conclusion
verified for non-deterministic disease progression: the probability of zero-
spreading events is reliably higher for human contact pattern compared to
randomized; however, the modulation of $\overline{n}_{\mathrm{inf}}$ with
$T_{\mathrm{lat}}$ is smeared out with increasing variability in the period
durations. Thus, in the unrealistic (but commonly adopted) limit of
exponentially distributed periods, human contact patterns still reduce the
robustness of outbreaks but no longer modulate the pace of epidemic spread.

To reproduce the relevant temporal features of human contact patterns, we
introduced non-Markovian contact dynamics in the form of Weibull-distributed
inter-encounter intervals (clustering) or inhomogeneous encounter rates
(cyclostationarity). Previous studies of non-Markovian disease spread [7]
found that clustering drastically affects the epidemic threshold for
$T_{\mathrm{lat}}=0$, which is caused by the high frequency of small inter-
encounter intervals [27] that, in our context, manifests as a near-zero peak
in the conditional encounter rate (Fig. 3). Although it was shown that some
non-Markovian models can be mapped onto effective Markov models [28, 29], our
results suggest that the non-Markovian and cyclostationary features of human
contact patterns make a similar general mapping elusive. This highlights the
necessity for generative models that are non-Markovian, yet well understood
and simple enough to find broad use in epidemic modeling and beyond.

Our work is a first step towards fulfilling this necessity: we identified
temporal statistics of real-world contact data that affect disease spread, and
faithfully reproduced them with our tailored Weibull renewal process. Thereby,
our work provides an accessible pathway towards including non-Poissonian
statistics into spreading processes, in general, and paves the way to
systematically study their non-equilibrium physics.

## Acknowledgements

We would like to thank Peter Sollich for helpful discussion. J.Z. received
financial support from the Joachim Herz Stiftung. J.Z., F.P.S., V.P., and
M.Wi. acknowledge funding by the Max Planck Society. V.P. acknowledges support
from the SFB 1528 ”Cognition of Interaction”.

## Methods

Extracting contacts from real-world physical proximity data: Consider data
composed of a list of co-locations (physical proximity) described by the tuple
(timestamp, user id A, user id B). We first sort the co-location times into
unique lists for all id pairs (A,B) and (B,A). For each valid A, we then
iterate over all lists (A,B) and merge successive co-location times to
construct pairwise contacts with starting time $s$ and duration $D$. Combining
these contacts yields a list of contacts $\\{(s_{i},D_{i})\\}_{A}$ for each
participant $A$.

From the lists of contacts, we construct a point-process-like representation
for each participant that we call encounter train (see Fig. 1). Throughout the
manuscript, an encounter refers to the starting time of a contact. The
encounter train of participant A is the time-sorted list of all contact
starting times and can formally be written as

$T(t)=\sum_{i}\delta(t-s_{i})$ (3)

The data from the Copenhagen Networks Study [12, 17] is based on Bluetooth
signals between phones of individuals that participated in the study. The
published data is a list of interactions described by the tuple (timestamp, A,
B, RSSI), where user id B can be negative if the interaction is with a device
outside of the study or an empty scan, and RSSI is the received (Bluetooth)
signal strength indicator. The RSSI can be considered as a proxy for
interaction distance, especially since all participants used the same device
[30], with an $\text{RSSI}\approx$-80\text{\,}\mathrm{d}\mathrm{B}\mathrm{m}$$
corresponding to a distance of about
$$2\text{\,}\mathrm{m}$\pm$0.5\text{\,}\mathrm{m}$$. Since the data provides a
maximal RSSI per time window, we consider
$\text{RSSI}<$-80\text{\,}\mathrm{d}\mathrm{B}\mathrm{m}$$ to indicate
interactions to be further apart than $2\text{\,}\mathrm{m}$ throughout the
full time window [30], and exclude them. Consequently, we filter the raw data
to only include those interactions that are within the study (user id B $\geq
0$) and have $\text{RSSI}\gtrsim$-80\text{\,}\mathrm{d}\mathrm{B}\mathrm{m}$$.
The data set covers a duration of
$t_{\mathrm{max}}=$28\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$, with
a temporal resolution of $\Delta
t=$5\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}$$, for 675 encounter trains.

Average time-dependent encounter rate $\rho(t)$: Because encounter trains are
a point-process-like representation, we can define an encounter rate as the
number of encounters in a window of size $\Delta t$. Assuming statistical
independence between weeks and between participants, we determine the average
time-dependent encounter rate $\rho(t)$ by averaging the number of encounters
in a time windows of size $\Delta t=$1\text{\,}\mathrm{h}$$ throughout the
week (i.e. first hour of a Sunday until last hour of a Saturday) across weeks
of the experiment and across participants. Statistical errors are calculated
on the level of participants using delete-$m$ jackknife error analysis.

Inter-encounter interval $\delta t$: To study temporal clustering and contact
bursts, we measure the interval $\delta t$ between consecutive encounter
times. The number of inter-encounter intervals per train ($m_{j}$) differs
because encounter trains have different numbers of encounter times
($n_{j}=m_{j}+1$). Since we are interested in the encounter statistics, each
encounter has the same statistical weight independent of its encounter train
origin. Consequently, the distribution $P(\delta t)$ is simply the
distribution over all observed intervals. Statistical errors are then
evaluated on the level of observed intervals using delete-$m_{j}$ jackknife
error analysis.

Conditional encounter rate $\Psi(\tau)$: To investigate how contact patterns
interact with disease spread, we measure the encounter rate $\Psi(\tau)$ upon
a hypothetical infection at $\tau=0$. To construct $\Psi(\tau)$, we iterate
over all encounters to measure the time-dependent encounter rate with temporal
resolution of the experiment, starting from the encounter time, i.e,
$\tau=t-s_{i}=0$, until $\tau=\tau_{\mathrm{max}}$ (we typically chose
$\tau_{\mathrm{max}}=$10\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$)
or, if $t_{\mathrm{max}}-s_{i}<\tau_{\mathrm{max}}$, until
$\tau=t_{\mathrm{max}}-s_{i}$. We then average over all these time-dependent
encounter rates taking into account their different lengths. Statistical
errors are calculated on the level of encounters using delete-$m_{j}$
jackknife error analysis.

Disease model: We consider a disease that progresses in three discrete states
upon infection: exposed-infectious-recovered. The duration $T_{\mathrm{lat}}$
within the exposed state is called latent period and the duration
$T_{\mathrm{inf}}$ within the infectious state is called infectious period.
For our main results, we consider the simple and intuitive case of a
deterministic disease progression, where these periods are always the same.
This is reminiscent of drawing the periods from delta distributions, which is
quite different to commonly employed approximations that draw periods from
exponential distributions (as expected for Poisson processes that describe
many state transitions, from radioactive decay to chemical reactions). To
control that our results also apply to non-deterministic disease progression,
we repeated our analysis for the more realistic case of gamma-distributed
periods [24, 25, 26] and obtained consistent results (Supplemental Material).

Potentially infectious encounters $n_{\mathrm{inf}}$: To avoid assumptions on
the probability of infection upon encounter, we introduce potentially
infectious encounters as the number of encounters that occur during the
infectious period of a hypothetical disease. For the deterministic disease
progression, we can enumerate the statistics by iterating over all encounters
of the data set. For each encounter $s_{i}$, we check whether the disease
progression still fits into the experimental duration $t_{\mathrm{max}}$
($s_{i}+T_{\rm lat}+T_{\rm inf}\leq t_{\mathrm{max}}$), and if true estimate
$n_{\mathrm{inf}}$ as the number of subsequent encounter $s_{j}$ for which
$T_{\rm lat}<s_{j}-s_{i}<T_{\rm lat}+T_{\rm inf}$. For the non-deterministic
disease progression, we need to sample disease realizations (see Supplemental
Material). Statistical errors are calculated on the level of encounters using
the delete-$m_{j}$ jackknife analysis.

Branching process with empirical distribution: To estimate the survival
probability from the empirical distribution of potentially infectious
encounters, $P(n_{\mathrm{inf}})$, we construct a discrete-time data-driven
branching process (Fig. 2c). In a first step, each infection causes $X\sim
P(n_{\mathrm{inf}})$ potentially infectious encounters. In a second step, each
of these encounters can cause a secondary infection with probability
$p_{\mathrm{inf}}$, such that the number of secondary infections is binomial,
$Y\sim\mathcal{B}(X,p_{\mathrm{inf}})$. From $Z_{t}$ infections in generation
$t$, we thus obtain $Z_{t=1}=\sum_{i=1}^{Z_{t}}Y_{i}$ infections in the next
generation.

Continuous-time branching process with inhomogeneous contacts: To study the
pace of epidemic spread, we construct a continuous-time branching process that
captures the conditional encounter rates but neglects interactions between
infected. Here, each infected individual generates independent encounter
trains starting from the initial infection time as inhomogeneous Poisson
processes with a time-dependent rate given by the conditional encounter rate
(Fig. 3d). Only those encounters that occur during the infectious period cause
secondary infections with a chosen probability $p_{\mathrm{inf}}$. Every
secondary infection then generates a new encounter train and so on. For
simplicity, we restrict our example to deterministic diseases with fixed
latent and infectious periods.

Point process models to approximate human contact patterns: To disentangle the
effect of distinct features of human contact patterns on the statistics of
encounters, we constructed point-process models that captured (i) the
distribution of rates across individuals, (ii), a time-dependent average
encounter rate, and (iii), the distribution of inter-encounter intervals, or a
combination thereof (see Supplementary Information for comparison of
combinations).

To reproduce the inter-individual variability, we match the number of
encounter trains to the data and weight each train with a relative rate, i.e,
$w_{i}=n_{\mathrm{train},i}/\langle n_{\mathrm{train}}\rangle$, where
$\langle\cdot\rangle$ is the average across trains.

To reproduce a time-dependent encounter rate, we employ thinning [31]: Given a
hidden process with rate $\rho_{i}=\max_{t}[\rho(t)]$, we accept events at
time $t$ with probability $p(t)=\rho(t)/\max_{t}[\rho(t)]$. This procedure can
formally only applied for memory-less hidden processes, i.e., Poisson
processes, in which case it results in an inhomogeneous Poisson process. To
further reproduce heterogeneous rates in the inhomogeneous Poisson process, we
fix $p(t)$ but rescale the rates of the hidden processes,
$\rho_{i}=w_{i}\max_{t}[\rho(t)]$.

To reproduce the empirical distribution of inter-encounter intervals, we
construct a Weibull renewal process: Inter-encounter intervals are drawn from
a Weibull distribution with scale parameter $\lambda$ and shape parameter $k$.
The Weibull distribution was parameterized by a fit to the data yielding
$(k,\lambda)=(0.36,3029.69)$. To further reproduce heterogeneous rates in the
Weibull-renewal process, we notice that the mean rate of a Weibull renewal
process is given by $\rho_{i}=\left[\lambda_{i}\Gamma(1+1/k_{i})\right]^{-1}$,
such that we can simply choose $k_{i}=k$ and $\lambda_{i}=\lambda/w_{i}$.

To combine all features in a single model, we construct a tailored renewal
process: Weibull renewal processes with heterogeneous rates and additional
(heuristic) thinning. We start with a set of hidden Weibull renewal processes
with $k_{i}=k$, $\lambda_{i}=\lambda/w_{i}$, and time-dependent acceptance
probability $p(t)$ with time-average $\overline{p(t)}$. The mean rate of each
process is $\rho_{i}=\overline{p(t)}w_{i}/\lambda\Gamma(1+1/k)$. Since we
cannot fit $(k,\lambda)$ of the hidden process, we further constrain the
parameters with the mean rate from data, i.e.,
$\overline{\rho(t)}=\langle\rho_{i}\rangle=\overline{p(t)}/\lambda\Gamma(1+1/k)$,
with $\langle w\rangle=1$ by construction. Since
$\overline{\rho(t)}/\overline{p(t)}=\max_{t}[\rho(t)]$, we thus find
$\lambda=\left[\max_{t}[\rho(t)]\Gamma(1+1/k)\right]^{-1}$, such that $k$
remains the only free parameter. We obtained our best estimate of $k$ by
minimizing the Kullback-Leibler divergence [32] between the distribution tails
($\delta t\gtrsim$0.5\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$) of
model and empirical $P(\delta t)$, finding $k\approx 0.24$.

Jackknife error estimation: To estimate statistical errors of our results, we
use jackknife error estimation while carefully taking into account the size of
the left-out data set. The basic idea of the jackknife method is to estimate
from some data $X=\\{x_{1},\ldots,x_{g}\\}$ the variance of an observable
$\hat{O}=f(X)$ using a systematic resampling approach [33]. Jackknife samples
$O_{j}$ are generated by systematically leaving out data, e.g.,
$\hat{O}_{j}=f(X_{\bar{\jmath}})$ with
$X_{\bar{\jmath}}=\\{x_{1},....,x_{j-1},x_{j+1},...x_{g}\\}$. Importantly,
here each $x_{i}$ can be a block of (differently many) data points. While
typically theses blocks have the same size $m$ (delete-$m$ jacknife), they
could have different sizes $m_{j}$ (delete-$m_{j}$ jackknife), which will be
relevant for some of our cases. From the jackknife samples, one can show that
bias-reduced estimators of the mean and variance are given by [34]

$\displaystyle\hat{O}_{J}$
$\displaystyle=\sum_{j=1}^{g}\frac{1}{h_{j}}\left(h_{j}\hat{O}-(h_{j}-1)\hat{O}_{j}\right),$
$\displaystyle\hat{\sigma}^{2}_{J}$
$\displaystyle=\frac{1}{g}\sum_{j=1}^{g}\frac{1}{h_{j}-1}\left(h_{j}\hat{O}-(h_{j}-1)\hat{O}_{j}-\hat{O}_{J}\right)^{2},$
(4)

where $h_{j}=(\sum_{j=1}^{g}m_{j})/m_{j}$, and $\hat{O}=f(X)$ is the naive
estimate. For blocks of equal size, $m_{j}=m$, we have $h_{j}=g$ and this
simplifies to

$\displaystyle\hat{O}_{J}$
$\displaystyle=g\hat{O}-\frac{g-1}{g}\sum_{j=1}^{g}\hat{O}_{j},$
$\displaystyle\hat{\sigma}^{2}_{J}$
$\displaystyle=\frac{g-1}{g}\sum_{j=1}^{g}\left(\hat{O}_{j}-\frac{1}{g}\sum_{j=1}^{g}\hat{O}_{j}\right)^{2}.$
(5)

In our case, the data $X$ is the set of all encounter trains and in the
resampling step we leave out individual encounter trains. Since trains include
differently many encounters, this can results in removing blocks of different
sizes. In particular, all observables that derive from the number of
encounters, e.g., $\overline{n}_{\text{inf}}$ or $P(n_{\text{inf}})$, require
the delete-$m_{j}$ analysis, (4), to estimate the statistical error. On the
other hand, for observables that depend on time-binned data, e.g., the time-
dependent rate, each encounter train has the same size given by the number of
time bins during the recording such that the delete-$m$ analysis, (5), is
sufficient to estimate the statistical error.

Data availability: Our methods are publicly available [35] and applied to
open-access data-sets [17, 13].

## References

  * Pastor-Satorras  _et al._ [2015] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Epidemic processes in complex networks, Rev. Mod. Phys. 87, 925 (2015).
  * Newman and Park [2003] M. E. J. Newman and J. Park, Why social networks are different from other types of networks, Phys. Rev. E 68, 036122 (2003).
  * Goltsev  _et al._ [2012] A. V. Goltsev, S. N. Dorogovtsev, J. G. Oliveira, and J. F. F. Mendes, Localization and Spreading of Diseases in Complex Networks, Phys. Rev. Lett. 109, 128702 (2012).
  * Ahn  _et al._ [2010] Y.-Y. Ahn, J. P. Bagrow, and S. Lehmann, Link communities reveal multiscale complexity in networks, Nature 466, 761 (2010).
  * Kitsak  _et al._ [2010] M. Kitsak, L. K. Gallos, S. Havlin, F. Liljeros, L. Muchnik, H. E. Stanley, and H. A. Makse, Identification of influential spreaders in complex networks, Nat. Phys. 6, 888 (2010).
  * Barabási [2005] A.-L. Barabási, The origin of bursts and heavy tails in human dynamics, Nature 435, 207 (2005).
  * Van Mieghem and van de Bovenkamp [2013] P. Van Mieghem and R. van de Bovenkamp, Non-Markovian Infection Spread Dramatically Alters the Susceptible-Infected-Susceptible Epidemic Threshold in Networks, Phys. Rev. Lett. 110, 108701 (2013).
  * Du  _et al._ [2018] Z. Du, S. J. Fox, P. Holme, J. Liu, A. P. Galvani, and L. A. Meyers, Periodicity in Movement Patterns Shapes Epidemic Risk in Urban Environments, ArXiv180905203 Phys. (2018), arXiv:1809.05203 [physics] .
  * Towers and Chowell [2012] S. Towers and G. Chowell, Impact of weekday social contact patterns on the modeling of influenza transmission, and determination of the influenza latent period, Journal of Theoretical Biology 312, 87 (2012).
  * Masuda and Holme [2013] N. Masuda and P. Holme, Predicting and controlling infectious disease epidemics using temporal networks, F1000Prime Rep 5 (2013).
  * Holme [2015] P. Holme, Modern temporal network theory: A colloquium, Eur. Phys. J. B 88, 234 (2015).
  * Sekara  _et al._ [2016] V. Sekara, A. Stopczynski, and S. Lehmann, Fundamental structures of dynamic social networks, Proc. Natl. Acad. Sci. 113, 9977 (2016).
  * Génois and Barrat [2018] M. Génois and A. Barrat, Can co-location be used as a proxy for face-to-face contacts?, EPJ Data Sci. 7, 1 (2018).
  * Schlosser  _et al._ [2020] F. Schlosser, B. F. Maier, O. Jack, D. Hinrichs, A. Zachariae, and D. Brockmann, COVID-19 lockdown induces disease-mitigating structural changes in mobility networks, Proc. Natl. Acad. Sci. 117, 32883 (2020).
  * Leitch  _et al._ [2019] J. Leitch, K. A. Alexander, and S. Sengupta, Toward epidemic thresholds on temporal networks: A review and open questions, Appl. Netw. Sci. 4, 105 (2019).
  * Presigny  _et al._ [2021] C. Presigny, P. Holme, and A. Barrat, Building surrogate temporal network data from observed backbones, Phys. Rev. E 103, 052304 (2021).
  * Sapiezynski  _et al._ [2019] P. Sapiezynski, A. Stopczynski, D. D. Lassen, and S. Lehmann, Interaction data from the Copenhagen Networks Study, Sci. Data 6, 1 (2019).
  * [18] SocioPatterns, http://www.sociopatterns.org.
  * Harris [1963] T. E. Harris,  _The Theory of Branching Processes_ (Springer Berlin, 1963).
  * Jiang  _et al._ [2013] Z.-Q. Jiang, W.-J. Xie, M.-X. Li, B. Podobnik, W.-X. Zhou, and H. E. Stanley, Calling patterns in human communication dynamics, Proc. Natl. Acad. Sci. 110, 1600 (2013).
  * Newman [2010] M. E. J. Newman,  _Networks: An Introduction_ (Oxford University Press, Oxford ; New York, 2010).
  * Lloyd-Smith  _et al._ [2005] J. O. Lloyd-Smith, S. J. Schreiber, P. E. Kopp, and W. M. Getz, Superspreading and the effect of individual variation on disease emergence, Nature 438, 355 (2005).
  * Scholtes  _et al._ [2014] I. Scholtes, N. Wider, R. Pfitzner, A. Garas, C. J. Tessone, and F. Schweitzer, Causality-driven slow-down and speed-up of diffusion in non-Markovian temporal networks, Nat. Commun. 5, 5024 (2014).
  * Bailey [1964] N. T. J. Bailey, Some Stochastic Models for Small Epidemics in Large Populations, Appl. Stat. 13, 9 (1964).
  * Anderson and Watson [1980] D. Anderson and R. Watson, On the spread of a disease with gamma distributed latent and infectious periods, Biometrika 67, 191 (1980).
  * Lloyd [2001] A. L. Lloyd, Destabilization of epidemic models with the inclusion of realistic distributions of infectious periods, Proc. R. Soc. Lond. B Biol. Sci. 268, 985 (2001).
  * Masuda and Holme [2020] N. Masuda and P. Holme, Small inter-event times govern epidemic spreading on networks, Phys. Rev. Res. 2, 023163 (2020).
  * Starnini  _et al._ [2017] M. Starnini, J. P. Gleeson, and M. Boguñá, Equivalence between Non-Markovian and Markovian Dynamics in Epidemic Spreading Processes, Phys. Rev. Lett. 118, 128301 (2017).
  * Feng  _et al._ [2019] M. Feng, S.-M. Cai, M. Tang, and Y.-C. Lai, Equivalence and its invalidation between non-Markovian and Markovian spreading dynamics on complex networks, Nat. Commun. 10, 3748 (2019).
  * Sekara and Lehmann [2014] V. Sekara and S. Lehmann, The Strength of Friendship Ties in Proximity Sensor Data, PLOS ONE 9, e100915 (2014).
  * Lewis and Shedler [1979] P. a. W. Lewis and G. S. Shedler, Simulation of nonhomogeneous poisson processes by thinning, Nav. Res. Logist. Q. 26, 403 (1979).
  * Kullback and Leibler [1951] S. Kullback and R. A. Leibler, On Information and Sufficiency, Ann. Math. Stat. 22, 79 (1951).
  * Efron [1982] B. Efron,  _The Jackknife, the Bootstrap, and Other Resampling Plans_ (Society of Industrial and Applied Mathematics, 1982).
  * Busing  _et al._ [1999] F. M. T. A. Busing, E. Meijer, and R. V. D. Leeden, Delete-m Jackknife for Unequal m, Stat. Comput. 9, 3 (1999).
  * noa [2021] Resonance_contact_disease, https://github.com/zierenberg/resonance_contact_disease (2021).
  * Xin  _et al._ [2022] H. Xin, Y. Li, P. Wu, Z. Li, E. H. Y. Lau, Y. Qin, L. Wang, B. J. Cowling, T. K. Tsang, and Z. Li, Estimating the Latent Period of Coronavirus Disease 2019 (COVID-19), Clin. Infect. Dis. 74, 1678 (2022).
  * Génois  _et al._ [2015] M. Génois, C. L. Vestergaard, J. Fournet, A. Panisson, I. Bonmarin, and A. Barrat, Data on face-to-face contacts in an office building suggest a low-cost vaccination strategy based on community linkers, Netw. Sci. 3, 326 (2015).

## Supplementary Information

In the Supplementary Information, we provide additional controls to verify the
robustness of our results in the main manuscript. In particular, we provide a
reanalysis of our main data set where we only consider participants for which
we can ensure a continuous participation in the experiment (Sec. III.4), as
well as an analysis of an additional, independent data set (Sec. III.5). Both
analyses fully support our main results. In addition, we provide complementary
analyses of our main data set. These include (i) a more detailed analysis of
shorter latent and infectious periods compared to the example considered in
the main manuscript that allows us to reach the low-dispersion regime (Sec.
III.2), and (ii) an analysis that verifies that our conclusions remain valid
for infections from outside of the study group (Sec. III.6).

###  III.1 Non-deterministic disease progression

Figure S1: Increased variability of disease-stage periods does not affect
conclusions on robustness of outbreaks but weaken modulation of spreading
pace. a, b: To include person-to-person variability, we draw both latent and
infectious periods from gamma distributions characterized by the dispersion
parameter $k$, which interpolates between exponential ($k=1$) and delta
($k\to\infty$) distributions. For small $k$, $T_{\rm lat}$ and $T_{\rm inf}$
differ in duration from realization to realization — across individuals, the
probability to be infectious at a given time is smeared out. As $k$ increases,
the periods vary less around their expected value and the disease progression
eventually becomes deterministic. c: Probability distributions of potentially
infectious encounters $\varepsilon_{\rm inf}$ for $k\to\infty$ (top) and
$k=10$ (bottom). When randomizing trains, the probability of zero-infectious
encounters is suppressed. d: Mean number of potentially infectious encounters
$\overline{\varepsilon}_{\rm inf}$ as a function of $k$ for the two examples
from the main manuscript
($T_{\mathrm{inf}}=$3\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$ and
$T_{\mathrm{lat}}=2$ or $6\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$).
For $k\to\infty$, we recover the result for deterministic disease progression
(dashed lines), where the latent period induces a notable difference between
the two examples. For smaller $k$, the resonance effects remain relevant but
the difference decreases, see also Fig. S2.

In the main manuscript, we have focused on deterministic disease progression,
where latent and infectious periods had a precise duration. For a more
realistic view, we want to allow for the latent and infectious periods to vary
from case to case (Fig. S1). To that end, we draw the periods $T_{\rm lat}$
and $T_{\rm inf}$ from a gamma distribution, but we keep the mean duration
fixed111Since now gamma-distributed periods are stochastic, we can no longer
enumerate the statistics of $P(n_{\mathrm{inf}})$, but have to sample it.
Specifically, we sample $10^{6}$ disease realizations, where we first draw a
random realization of the disease progression $(T_{{\rm lat},i},T_{{\rm
inf},i})$, to then draw disease start times $s^{i}$ from the ensemble of all
encounters until we find an $s_{i}$ such that the disease progression is
within the remaining duration of the experiment, i.e., $s_{i}+T_{{\rm
lat},i}+T_{{\rm inf},i}\leq T_{\rm exp}$. Only once we have a valid disease
start time $s_{i}$, we count the number of subsequent encounters within the
infectious period as above. If, for any disease realization, we need to draw
more than 1000 disease start times until we find a valid one, we abort the
estimation for the set of parameters $(T_{\rm lat},T_{\rm inf},k)$. By this
procedure to first draw and fix a random realization of the disease
progression, we avoid a bias towards small periods that would occur due to the
finite period of the experimental data. . The case-to-case variability is then
parameterized through the dispersion parameter $k$ (Fig. S1a).

In particular, we fix the mean $\langle T_{i}\rangle$ by choosing the scale
$\theta=\langle T_{i}\rangle/k$ such that
$P_{k}(T_{i})=T_{i}^{k-1}e^{-T_{i}/\theta}/[\Gamma(k)\theta^{k}]$. The
dispersion parameter allows us to interpolate between a delta distribution
($k\to\infty$) and an exponential distribution ($k=1$), as commonly assumed in
computational epidemiology for mathematical tractability [1].

Clinically observed distributions of periods between disease states are
neither delta distributed nor exponential distributed and may be best
described by distributions with a clear peak but vanishing probability at
zero, such as log-normal distributions or gamma distributions with dispersion
parameters in between $(1,\infty)$. On the one hand, delta-distributed periods
seem like a convenient but unrealistic simplification. On the other hand,
exponentially distributed periods may appear more realistic, however, they
imply an artificially high probability of short durations, which in turn leads
to realizations where the infectious period either starts shortly upon
infection or has close-to-zero duration (example traces in Fig. S1a). In fact,
it has been argued already in the past that gamma distributions are more
realistic [24, 25, 26], see for example empirical distributions of latent
periods for COVID-19 [36], such that more realistic dispersion parameters
could be in the range $h\in[5,20]$ which is between delta and exponential. To
investigate how case-to-case variability affects the number of potentially
infectious encounters $\varepsilon_{\rm inf}$, we again consider the
probability distribution $P(\varepsilon_{\rm inf})$ (Fig. S1c), and revisit
the two examples ($T_{\rm
inf}=$3\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$ with either $T_{\rm
lat}=$2\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$ or $T_{\rm
lat}=$6\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$).

For the delta-disease ($k\to\infty$, top), the red ($T_{\rm
lat}=$6\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$) and blue ($T_{\rm
lat}=$2\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$) distributions
exhibit a peak at zero, they are broad, and have a long tail. Comparing red
and blue, we find that the latent period determines the height of the peak at
zero as well as the position of the bulk distribution, and, thereby,
determines the mean number of potentially infectious encounters (dashed
vertical lines). The mean values clearly differ. Again comparing with the
randomized encounter trains (yellow line), we find no peak at zero, a shorter
tail, and no dependence on the latent period (the respective randomized lines
fall on top of each other). Importantly, a peak at zero implies that the
infected individual does not pass on the infection, so that the disease
becomes more likely to die out if case numbers are low.

Changing to the non-deterministic disease progression ($k=10$, Fig. S1c
bottom), the distribution from randomized data is barely affected. However,
the red and blue distributions are more similar to each other, but also
broader and smoother than for the the delta-disease. This can be explained by
gamma-distributed periods acting as a smoothing kernel along both dimensions
of Fig. S1f, where variability in the infectious period directly affects
$\overline{\varepsilon}_{\text{inf}}$, while variability in the latent period
affects $\overline{\varepsilon}_{\text{inf}}$ through the resonance effect.
Consequently, we expect that with decreasing dispersion parameter $k$, the
mean number of potentially infectious encounters becomes independent of the
mean latent period. Indeed, when considering
$\overline{\varepsilon}_{\text{inf}}$ as a function of $k$, we find that the
estimates for our two examples approach each other, as $k$ decreases (Fig.
S1d).

Note that our analysis has a lower bound in $k$ once realizations of disease
progression (latent + infectious period) cannot find sufficiently many initial
encounters to fit into the finite duration of the experiment (4 weeks for
CNS). However, using other examples with smaller latent and infectious periods
(where we can acquire enough statistics), we show that the two extreme cases
meet for $k\approx 1$ (see Fig. S2).

###  III.2 Non-deterministic disease progression with fast disease stages

With the example of the main text, we were not able to sample the dispersion
relation for small $k$ values (more variability across disease realizations)
because the 28-day duration of the data becomes too short once the periods of
disease progression are close-to exponentially distributed ($k\to 1$).

To avoid this issue and to illustrate shorter timescale, we here compare with
another hypothetical example of a “faster” disease progression, where
$T_{\text{inf}}=$0.5\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$ and
$T_{\text{lat}}$ is either $1$ or $1.5$ days (Fig. S2). In this case, expected
periods are very short so that even realizations with periods that severely
exceed their expected value fit into the 28-day duration.

We find a few noteworthy aspects: First, the absolute value of potentially
infectious encounters $n_{\mathrm{inf}}$ is much lower for faster disease
progression. This is due to the much shorter infectious period. However, the
relative deviation from the randomized surrogate data is consistent between
the two examples. Second, we now find that the longer latent period ($1.5$ vs.
$1.0$) leads to fewer potentially infectious encounters (red vs. blue). This
is just a result of the chosen latent periods: choosing an even longer latent
period would result in a respective decrease of $n_{\mathrm{inf}}$ (e.g. $2.0$
vs $1.5$). Third, as expected, as $k\to 1$, the estimates of
$n_{\mathrm{inf}}$ overlap for different latent periods, because individual
disease realizations become very variable (cf. Fig. S1a).

Figure S2: Example with smaller latent and infectious periods provides
additional insight on low-dispersion regime of non-deterministic disease
progression. Here, we compare the example from the main manuscript (“slow”,
$T_{\text{inf}}=$3\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$ and
$T_{\text{lat}}$ either $2$ or $6$ days) with a hypothetical “fast” disease
progression
($T_{\text{inf}}=$0.5\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$$ and
$T_{\text{lat}}$ either $1$ or $1.5$ days) and show both the mean number of
potentially infectious encounters as a function of the dispersion parameter
$k$ of the gamma-distributed latent and infectious period (top) as well as
their distributions for selected $k$. Due to the finite duration of the
recording, the accessible low-dispersion regime is determined by the mean
latent and infectious period, because for low dispersion large periods quickly
exceed the finite duration. For faster disease progression (smaller latent and
infectious period), we observe modulations on smaller timescales and, in
addition, reach the low-dispersion regime of exponentially distributed periods
($k=1$) commonly assumed in epidemiological simulations. As one can see more
clearly for faster disease progression, the mean number of potentially
infectious encounters approach each other in this low-dispersion regime, which
can only be anticipated from the results for slower disease progression. This
implies that modulations will not be present for exponentially distributed
latent and infectious periods but only for more realistic non-exponentially
distributed ones with higher dispersion. In addition, one can see that
estimates of mean potentially infectious encounters for randomized surrogate
data in both cases of slow and fast disease progression are independent of
dispersion, confirming again that randomized encounter trains destroy temporal
structure and the subsequent modulations.

###  III.3 Simple point processes cannot fully reproduce temporal features

As we noted in the main text, each of the considered simple generative models
is insufficient to reproduce the full set of observed temporal features of
contact patterns. Here, we provide a more complete overview of processes and
the features they reproduce (Fig. S3). We identified three relevant features:

i) a time-varying, cyclostationary encounter rate $\rho(t)$ (first column)

ii) a heavy-tailed inter-encounter-interval distribution $P(\delta t)$ (second
column)

iii) encounter rates vary between individuals (third column)

Feature i is responsible for modulation of the conditional encounter rate
$\Psi(\tau)$ (fourth column). Combined with disease progression, this causes
modulations of $\overline{n}_{\mathrm{inf}}$ and distributions
$P(n_{\mathrm{inf}})$ that vary with $T_{\mathrm{lat}}$ and $T_{\mathrm{inf}}$
(last two columns). Feature ii is responsible for zero-spreading events, which
destabilize epidemic outbreaks and manifest in the peak of
$P(n_{\mathrm{inf}})$ for $n_{\mathrm{inf}}\approx 0$ (last column). Feature
iii is responsible for super-spreading events, which correspond to the tail of
$P(n_{\mathrm{inf}})$ and, thus, can cause a systematic shift of
$\overline{n}_{\mathrm{inf}}$.

Figure S3:  Simple generative processes are not capable of reproducing all
aspects of the full real-world contact statistics. An inhomogeneous Poisson
process can reproduce cyclostationary $\rho(t)$ but fails to reproduce
sufficiently heavy-tailed $P(\delta t)$. A Weibull renewal process reproduces
heavy-tailed $P(\delta t)$ but fails to reproduce the cyclostationary
$\rho(t)$. Generative processes need to reproduce variability between
individuals $P(n_{\textrm{train}}$) to produce long-tailed distributions
$P(n_{\mathrm{inf}})$. In order to reproduce the data in all the here
considered aspects, we had to designed a tailored Weibull renewal process with
homogeneous rates.

###  III.4 Control regarding continuous contribution of participants

In the main manuscript, we use the full published data set of the Copenhagen
Networks Study [17], covering the physical proximity data of 675 participants.
Upon closer inspection, there are periods both at the beginning and the end of
the experiment without entries for some of these 675 participants. Since
entries also occur for Bluetooth signals with unknown devices, this may
indicate irregularities in the contact behavior of some of the participants,
e.g., incomplete participation of individuals.

In order to make sure that our results are not affected by such boundary
effects, we reanalyzed the data and included only the contact trains of those
individuals for which any Bluetooth signal was recorded on both the first and
last day of the study (Fig. S4). Technically, we achieved this easily by
restricting our analysis to those IDs for which timestamps were recorded
within the first day (timestamp $<1\cdot 24\cdot 60\cdot
60$\text{\,}\mathrm{s}$$) and the last day (timestamp $>27\cdot 24\cdot
60\cdot 60$\text{\,}\mathrm{s}$$), reducing the data set to 533 contact
trains.

This control analysis fully supports our quantitative results from the main
manuscript (Fig. S4) such that we can rule our artifacts from boundary effects
of incomplete participation. In particular, we observe a matching weekly
structure of the encounter rate (Fig. S4a), a matching distribution of inter-
encounter intervals that can be fitted with a Weibull distribution (Fig. S4b),
and a matching conditional encounter rate (Fig. S4c). Consequently, both mean
potentially infectious encounters for deterministic disease progression (Fig.
S4d and f) as well as for non-deterministic disease progression (Fig. S4e and
g) match our main results.

Figure S4: Control regarding continuous contribution of participants. Here, we
excluded those trains that did not have any encounters during the first or
last day. This results in 533 instead of 675 trains. Panels match key figures
from the main manuscript. Results are consistent. The fitted Weibull
distribution has shape parameter $0.37$ and scale parameter $3161.95$.

###  III.5 Analysis of an alternative data set

To further challenge the robustness of our results, we repeated our analysis
on a completely independent data set. Here, we consider contact data recorded
at one of the office buildings of the French Institute for Public Health
Surveillance “InVS” [37]. This data is recorded with a different technique,
namely so-called near-field chips that only record signals in close proximity
($\lesssim$5\text{\,}\mathrm{m}$$) and avoid to threshold the Bluetooth
signal. Moreover, the temporal resolution of contacts is
$20\text{\,}\mathrm{s}$ as opposed to
$5\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}$ in the main manuscript. In
addition, the data is recorded for a different social group, namely adults
within an office building. Last, the data is recorded in another country
(France) by a different collaboration (SocioPatterns). The data set spans 2
weeks of recording 145 participants (two thirds of the staff agreed to
participate).

The analysis of this completely independent data set provides completely
consistent results to those presented in the main manuscript (Fig. S5). When
comparing the results, we have to highlight that the available statistics for
this data set are much smaller due to smaller duration and smaller number of
participants. However, we clearly see the expected weekly structure in the
encounter rate (which is here again dominated by working days because of
office hours), the distribution of inter-encounter intervals that is well
described by a Weibull distribution, as well as the typical conditional
encounter rate with a peak at 7 days. Consequently, also the results for
(delta) disease progression are consistent with our main findings on the
modulation of potential secondary infections (Fig. S5d). Results on non-
deterministic disease progression are confined to shorter latent and
infectious periods due to the shorter experimental duration. We conclude that
the additional data set fully supports our results in the main manuscript.

Figure S5:  Main results for other data set. Because the recording lasted
“only” two weeks, the duration of disease that can be sampled to estimate the
dispersion were limited to the fast disease progression. The fitted Weibull
distribution has shape parameter $0.22$ and scale parameter $675.03$.

###  III.6 Infections from outside the study group

To check the effect of contacts that could take place with people who were not
part of the study, we investigate disease onsets that do not directly follow
the contact patterns observed in the data (Fig. S6). In our analysis of the
main manuscript, an infection could only originate from an encounter with
another participant in the data set. Here, we keep the original encounter
trains (to evaluate potentially infectious encounters) but the disease onset
can occur at any time, due to a hypothetical contact with an external person.
We focus on the resulting distribution and the mean of potentially infectious
encounters $n_{\rm inf}$ (Fig. S6).

We distinguish the following different possibilities of disease onset:

internal

Onsets occur as in the main manuscript only at encounters recorded in the data
set. This naturally incorporates the spatio-temporal structure of encounters,
in particular their temporal inhomogeneity and their variability across
participants.

external i)

Onsets occur completely random, at random times for random participants. This
neglects both temporal inhomogeneity of encounters and their variability
across individuals.

external ii)

Onsets occur at random times with probability proportional to the encounter
rate for random participants with probability proportional their total number
of encounters. This incorporates both the (averaged) temporal inhomogeneity of
encounters and the variability across individuals.

external iii)

Onsets occur at random times with probability proportional to the encounter
rate but for uniformly random participants. This incorporates the (averaged)
temporal inhomogeneity of encounters but neglects the variability across
individuals.

external iv)

Onsets occur at uniformly random times for random participants with
probability proportional their total number of encounters. This neglects the
temporal inhomogeneity of encounters but incorporates the variability across
individuals.

Once an onset has been chosen, the disease progression is modeled as in the
main manuscript. We focus on the usual examples but use gamma-distributed
latent and infectious periods with $k=10$ to evaluate potentially infectious
encounters ($k=10$ aims to strike a balance between completely deterministic
and exponential disease stages [24, 25, 26]).

Comparing the different version of disease onset (Fig. S6), we can attribute
clear effects to both the temporal inhomogeneity of the onset time as well as
the variability of onset times across individuals. Please note that in all
cases the encounter statistics of the actual encounter trains did not change —
all we change is the statistics of the disease onset time. Please note further
that the results for different version of disease onset represent the extreme
scenario where all disease onsets originate from external sources.

Comparing the distributions of $n_{\text{inf}}$ for fixed $T_{\text{lat}}$
(Fig. S6a and Fig. S6c each), we notice that those distributions that best
resemble the shape of internal disease onset are those where external disease
onset statistics incorporate the variability across individuals (external ii
and iv). This can be explained by the fact that also for internal disease
onset more onsets occur for contact trains with more encounters, which in turn
increases the probability of higher $n_{\text{inf}}$ and thereby also the mean
$n_{\text{inf}}$. It appears that for the overall shape of the distribution,
as well as the leading order of its mean value, it is not necessary that
disease onsets occur with the same temporal inhomogeneity as true encounters
for the chosen infectious periods (this may change for very small infectious
periods though).

Comparing further the results of specific disease onsets for different
$T_{\text{lat}}$ (Fig. S6b comparing solid vs opaque symbols), we notice that
incorporating the temporal inhomogeneity into the disease onset (external
disease ii and iii) seems relevant for the modulation that causes differences
in $\overline{n}_{\text{inf}}$ due to latent periods. This can be explained by
reoccurring contact patterns: if infections are not more likely to occur
during times of many contacts (as for i and iv), then differences in the
number of potential infections are smeared out when averaging over multiple
realizations.

Concluding, the results for the extreme scenarios of external infections fully
support our conclusions from the main manuscript.

Figure S6:  Infectious encounters for external infections In our analysis in
the main manuscript, we preserved the temporal features by constraining
disease onsets to available encounters. a: Distribution of potentially
infectious encounters for non-deterministic disease progression ($k=10$). 6
days latent period. b: Comparison of the mean $\bar{n}_{\rm inf}$ between 6
days latent period (dark) and 2 days latent period (light). c: Same as a), but
2 days latent period.

###  III.7 Impact of the distance threshold for considered contacts

In the main manuscript, we filtered the physical proximity data to only
include contacts with a minimal Bluetooth signal strength
($\text{RSSI}\gtrsim$-80\text{\,}\mathrm{d}\mathrm{B}\mathrm{m}$$) as a proxy
for contact distances below $2\text{\,}\mathrm{m}$ [30]. This choice was
motivated in the main manuscript and moreover corresponds to a value that
optimizes the ratio between strong and weak links [30]. Here we show that our
results are robust under changes in the cut-off signal strength which sets a
scale, but does not alter the patterns observed. In particular, we consider
two cases of lower
($\text{RSSI}\gtrsim$-75\text{\,}\mathrm{d}\mathrm{B}\mathrm{m}$)$ and higher
($\text{RSSI}\gtrsim$-95\text{\,}\mathrm{d}\mathrm{B}\mathrm{m}$)$ threshold.
As can be seen from Fig. S7, this threshold changes the scale of the encounter
rate (Fig. S7a) but does not alter the resonance effect compared to random
encounters (Fig. S7b).

Figure S7:  To test which impact the distance at which a contact becomes
potentially infectious has on our results, we varied the RSSI threshold in the
preprocessing step. A higher threshold
($\text{RSSI}\gtrsim$-95\text{\,}\mathrm{d}\mathrm{B}\mathrm{m}$$) corresponds
to including more contacts (which previously would be considered too far
away), thus resulting in a higher number of potentially infectious encounters,
and (in absolute numbers) to a higher conditional encounter rate. Vice versa,
a lower threshold
($\text{RSSI}\gtrsim$-75\text{\,}\mathrm{d}\mathrm{B}\mathrm{m}$$,
corresponding to about 1 meter distance) leads to a lower conditional
encounter rate. In both cases, the overall functional shape of the conditional
encounter rate (featuring valleys at night, peaks at daytime and pronounced
peaks at seven days) remains intact. Consequently, also the number of
potentially infectious encounters relative to randomized remain mostly
unaltered.

###  III.8 Impact of the contact duration for considered contacts

In the main manuscript, we filtered contacts to only include those with a
sufficient duration of at least
$15\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{u}\mathrm{t}\mathrm{e}\mathrm{s}$.
This choice was motivated in the main manuscript, and it is particularly
suitable to avoid transient contacts that arise in passing. Still, different
choices are possible. Here, we show for shorter and longer minimum durations
that this choice only affects the scale of the encounter rate but does not
affect any results relating to the modulation of potentially infectious
encounters (Fig. S8).

Figure S8:  To test which impact the minimum duration after which a contact
becomes potentially infectious has on our results, we varied the required
duration for qualifying contacts during the preprocessing step. A shorter
required duration (e.g.
$5\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{u}\mathrm{t}\mathrm{e}\mathrm{s}$
vs.
$15\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{u}\mathrm{t}\mathrm{e}\mathrm{s}$
in the main manuscript) leads to more contacts qualifying and a higher
conditional encounter rate. Vice versa, a longer required duration (e.g.
$30\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{u}\mathrm{t}\mathrm{e}\mathrm{s}$)
leads to a lower conditional encounter rate. In both cases, the overall
functional form of the conditional encounter rate remains intact. Again,
although we see a higher absolute number of potentially infectious encounters,
the change relative to randomized encounter trains is mostly unaltered.