{'fraction_non_alphanumeric': 0.03578290915219956, 'fraction_numerical': 0.006573403000168549, 'mean_word_length': 4.892442149170722, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
How contact patterns destabilize and modulate epidemic outbreaks

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.


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, realworld 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 realworld data can be partly remedied by generating surro- * JZ and FPS contributed equally gate 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. 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.

II. RESULTS
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 meters and longer than 15 minutes, 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 τ is the time elapsed since an infection, the disease can be either in a latent phase τ ∈ [0, T lat ), during which individuals are infected but not infectious, or in an infectious phase τ ∈ [T lat , T lat + T inf ), during which the disease can be transmitted with probability p 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 inf (Fig. 2a). We estimate the statistics of n inf by sampling hypothetical infections that may start at any encounter of the data set. For the case of fixed T lat and T inf , this can be solved by a simple enumeration.

Human contact patterns destabilize epidemic outbreaks
We start by comparing the distributions P (n inf ) obtained for two examples of disease progression with the same T inf but different T lat (Fig. 2b). We find that human contact patterns lead to high probabilities of very small n inf (zero-spreading events) and very large n inf (superspreading 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 inf ), we construct a data-driven branching process (Fig. 2c). Here, each infection first generates n inf encounters according to P (n inf ), each of which is infected independently with probability p inf resulting in binomialdistributed secondary infections. By summing over all 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 inf , are estimated from encounter trains for a disease model with latent period T lat and infectious period T inf . b: P (n inf ) depends on T lat for human encounter trains (red, blue) but not for randomized encounter trains (yellow lines), which underestimate both zerospreading as well as super-spreading events. c: Data-driven branching model where n inf are drawn from the empirical distribution P (n inf ) and get infected with constant probability p 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 inf in (b). possible n inf , this yields the offspring distribution For branching processes, one can calculate the asymptotic extinction probability with the help of the probability generating function, π(θ) = ∞ x=0 p(x)θ x : The extinction probability is the smallest θ * that solves θ = π(θ) [19], which we solve numerically because P (n inf ) is an empirical distribution.
For our chosen examples of T 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 where n inf is the average number of potentially infectious encounters. (Note that this representation accounts for the surprising differences in n 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 ≤ 1, P ext = 1 such that asymptotically no outbreak occurs, and for R 0 > 1, P ext < 1 such that there is a non-zero chance for an outbreak. In addition, we observe that using P (n inf ) from human encounter trains leads to a higher extinction probability than using P (n inf ) from randomized encounter trains. Whereas the extinction probability for randomized trains is independent of T lat , we find for human encounter trains that T lat determines the extinction probability by changing the probability of low n 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 As noted above, we observe that n inf depends on T lat for human encounter trains. This is at first glance surprising, because for memory-less Poisson processes n inf is proportional to T inf but independent of T lat . Hence, in the following we systematically vary T lat and T inf to study how the interplay between human contact patterns and disease progression affects n inf , which via Eq. (2) also affects R 0 .
For a fixed T inf = 3 days, we observe a systematic dependence of n inf on T lat for human encounter trains that is not present for randomized encounter trains (Fig. 3a). For human encounter trains n inf can be higher (T lat < 1 and 4 < T lat < 8) or lower (1 < T lat < 4) as compared to the randomized case, which suggests a periodic modulation.
To further study the dependence on T inf , we consider n inf relative to randomized (Fig. 3b). This accounts for the trivial increase of n inf with increasing T inf . For small T inf < 1 day, we find daily modulations as a function of T lat , with regions below-randomized (blue) and aboverandomized (red). This effect diminishes with increasing T inf but remains as daily sub-modulation. For larger T inf , we find extended triangular regions that are belowrandomized (blue) or above-randomized (red). The inter- Real-world contact patterns modulate the pace of epidemic spread as a function of latent period. a: Absolute n inf are a periodic function of T lat for human encounter trains but not for randomized (yellow). b: Relative n inf /n rand inf reveals periodic modulations on both daily and weekly scale in the full (T lat , T inf )-plane. c: Ψ(τ ) 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 n inf by combinations of T lat and T 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 Ψ(τ ) and infected with probability p inf for R0 average secondary infections. e: Initializing the model with I0 = 100 random infections in [−T lat −T inf , 0), the (average) number of infections grows exponentially. The pace of epidemic spread decreases with T lat but is modulated for non-Poissonian Ψ(τ ). f: The growth rate λ for time-independent encounter times in a fixed T inf is expected to decrease trivially with T lat . If not constant, Ψ(τ ) modulates λ and causes regimes of slower-than-random (blue) or faster-than-random (red) pace of epidemic spread.
faces between these regions are located at T lat + T inf = 7 days and at T lat = 7 days, which suggests that these are weekly modulations. Summarizing, we find periodic modulations of n inf as a function of T lat , with periods on the order of days for small T inf and on the order of weeks (with daily sub-modulations) for larger T inf .
In order to understand how changes in T lat result in periodic modulations of n inf , we introduce the conditional encounter rate Ψ(τ ) 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, Ψ(τ ) measures the rate of potentially infectious encounters, but neglects fluctuations across hypothetical infections (cf. Fig. 2b). We find that Ψ(τ ) features a peak at 0 (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 Ψ(τ ) marks periods of statistically high and low encounter rates after infection, and thereby implies a timedependent secondary attack rate. From Ψ(τ ), we can obtain an estimate of n inf for a particular disease pro-gression (T lat , T inf ) by integrating within the infectious period (Fig. 3c, shaded areas). For the two previous examples with fixed T inf = 3 days, we reproduce a lower n inf for T lat = 2 days (blue area covering the valley) and a larger n inf for T lat = 6 days (red area covering the 7-day peak). We thus find that the alignment of the infectious period with regions of low or high Ψ(τ ) determines n inf .
Since n 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 Ψ(τ ). During the infectious period, encounters again have a probability p 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 lat − T inf , 0]. The time evolution of daily new cases shows clear exponential growth (Fig. 3e), where the growth rate λ trivially decreases with the generation time and, thus, T lat (Fig. 3f). However, this expected decrease of λ for Poisson contacts (yellow) is modulated by Ψ(τ ) for human contact patterns, which results in slower-thanrandom (blue) or faster-than-random (red) growth.
Summarizing, the dependence of n inf on T lat implies that human contact patterns modulate the basic reproduction number and thereby the growth rate of an epidemic outbreak.
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 ρ(t) averaged across individuals and weeks is timedependent but cyclostationary; ρ(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 n inf (and Ψ(τ ), see Supplemental ii) The distribution of inter-encounter intervals P (δt) has high probability for small δt and a heavy tail of nonvanishing probability for large δt (second column). Because this tail corresponds to long periods without any encounters, it causes the high probability of n inf ≈ 0 (last column) that is responsible for the destabilization of epidemic outbreaks. P (δt) is dominantly shaped by the clustering of human contacts and can be well approximated by a Weibull distribution [7,20]. Accordingly, a Weibullrenewal processes (last row) well reproduces P (δt) and P (n inf ), but it does not have a time-varying ρ(t) and cannot reproduce the period modulations of Ψ(τ ) and n 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 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 superspreading events, it inevitably causes periods of zero contacts that we find to increase the probability of zerospreading 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 zerospreading 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 lat , T 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 n inf with T lat is smeared out with increasing variability in the period durations. Thus, in the un- FIG. 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 ρ(t) (first column), heavy-tailed P (δt) (second column), and heterogeneous P (ntrain) (third column), reproduces core characteristics such as modulation of n inf (fourth column, cf. Fig. 3 and Supplementary Material for additional figures for Ψ(τ )) and increased probability of zero-spreading in P (n inf ) (fifth column, cf. Fig. 2). center: An inhomogeneous Poisson process with only time-varying, cyclostationary ρ(t) and heterogeneous P (ntrain) reproduces the modulations in n inf , but underestimates the probability of low n inf due to a lack of clustering and the corresponding large inter-cluster times in P (δt). bottom: A Weibull renewal process with contact clustering from P (δt) and heterogeneous P (ntrain) reproduces the high probability of rare outcomes in P (n inf ), but cannot reproduce the modulations in n inf due to a lack of cyclic temporal structure. If, additionally, one relaxes the constraint of heterogeneous P (ntrain) and considers statistically identical trains, then both surrogate processes underestimate the probability of large n inf related to super-spreading (see Supplementary Material). realistic (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 interencounter 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 lat = 0, which is caused by the high frequency of small interencounter 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 colocation 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 From the lists of contacts, we construct a pointprocess-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 timesorted list of all contact starting times and can formally be written as 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 RSSI ≈ −80 dBm corresponding to a distance of about 2 m ± 0.5 m. Since the data provides a maximal RSSI per time window, we consider RSSI < −80 dBm to indicate interactions to be further apart than 2 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 ≥ 0) and have RSSI −80 dBm. The data set covers a duration of t max = 28 days, with a temporal resolution of ∆t = 5 min, for 675 encounter trains.
Average time-dependent encounter rate ρ(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 ∆t. Assuming statistical independence between weeks and between participants, we determine the average time-dependent encounter rate ρ(t) by averaging the number of encounters in a time windows of size ∆t = 1 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 δt: To study temporal clustering and contact bursts, we measure the interval δ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 (δ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 Ψ(τ ): To investigate how contact patterns interact with disease spread, we measure the encounter rate Ψ(τ ) upon a hypothetical infection at τ = 0. To construct Ψ(τ ), 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, τ = t − s i = 0, until τ = τ max (we typically chose τ max = 10 days) or, if t max − s i < τ max , until τ = t 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: exposedinfectious-recovered. The duration T lat within the exposed state is called latent period and the duration T 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 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 max (s i + T lat + T inf ≤ t max ), and if true estimate n inf as the number of subsequent encounter s j for which T lat < s j − s i < T lat + T 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 inf ), we construct a discrete-time data-driven branching pro-cess (Fig. 2c). In a first step, each infection causes X ∼ P (n inf ) potentially infectious encounters. In a second step, each of these encounters can cause a secondary infection with probability p inf , such that the number of secondary infections is binomial, Y ∼ B(X, p inf ). From Z t infections in generation t, we thus obtain Z t=1 = Zt i=1 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 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 train,i / n train , where · is the average across trains.
To reproduce a time-dependent encounter rate, we employ thinning [31]: Given a hidden process with rate ρ i = max t [ρ(t)], we accept events at time t with probability p(t) = ρ(t)/ max t [ρ(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, To reproduce the empirical distribution of interencounter intervals, we construct a Weibull renewal process: Inter-encounter intervals are drawn from a Weibull distribution with scale parameter λ and shape parameter k. The Weibull distribution was parameterized by a fit to the data yielding (k, λ) = (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 ρ i = [λ i Γ(1 + 1/k i )] −1 , such that we can simply choose k i = k and λ i = λ/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, λ i = λ/w i , and time-dependent acceptance probability p(t) with time-average p(t). The mean rate of each process is ρ i = p(t)w i /λΓ(1 + 1/k). Since we cannot fit (k, λ) of the hidden process, we further constrain the parameters with the mean rate from data, i.e., ρ(t) = ρ i = p(t)/λΓ(1 + 1/k), with w = 1 by construction. Since ρ(t)/p(t) = max t [ρ(t)], we thus find λ = [max t [ρ(t)]Γ(1 + 1/k)] −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 (δt 0.5 days) of model and empirical P (δt), finding k ≈ 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 , . . . , x g } the variance of an observableÔ = f (X) using a systematic resampling approach [33]. Jackknife samples O j are generated by systematically leaving out data, e.g.,Ô j = f (X) with X = {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 biasreduced estimators of the mean and variance are given by [34] where h j = ( g j=1 m j )/m j , andÔ = f (X) is the naive estimate. For blocks of equal size, m j = m, we have h j = g and this simplifies tô 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., n inf or P (n 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 [13,17].

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 D), as well as an analysis of an additional, independent data set (Sec. III E). 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 B), and (ii) an analysis that verifies that our conclusions remain valid for infections from outside of the study group (Sec. III F).

A. Non-deterministic disease progression
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 lat and T inf from a gamma distribution, but we keep the mean duration fixed 1 . The case-to-case variability is then parameterized through the dispersion parameter k (Fig. S1a).
FIG. 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 → ∞) distributions. For small k, T lat and T 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 ε inf for k → ∞ (top) and k = 10 (bottom). When randomizing trains, the probability of zero-infectious encounters is suppressed. In particular, we fix the mean T i by choosing the scale θ = T i /k such that P k (T i ) = T k−1 i e −Ti/θ /[Γ(k)θ k ]. The dispersion parameter allows us to interpolate between a delta distribution (k → ∞) 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, ∞). 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 ∈ [5,20] which is between delta and exponential. To investigate how case-to-case variability affects the number of potentially infectious encounters ε inf , we again consider the probability distribution P (ε inf ) (Fig. S1c), and revisit the two examples (T inf = 3 days with either T lat = 2 days or T lat = 6 days).
For the delta-disease (k → ∞, top), the red (T lat = 6 days) and blue (T lat = 2 days) 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 ε inf , while variability in the latent period affects ε 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 ε 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 ≈ 1 (see Fig. S2).

B. 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 → 1).
To avoid this issue and to illustrate shorter timescale, we here compare with another hypothetical example of a "faster" disease progression, where T inf = 0.5 days and T 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 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 inf (e.g. 2.0 vs 1.5). Third, as expected, as k → 1, the estimates of n inf overlap for different latent periods, because individual disease realizations become very variable (cf. Fig. S1a).
FIG. 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 inf = 3 days and T lat either 2 or 6 days) with a hypothetical "fast" disease progression (T inf = 0.5 days and T 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.
C. 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 ρ(t) (first column) ii) a heavy-tailed inter-encounter-interval distribution P (δt) (second column) iii) encounter rates vary between individuals (third column) Feature i is responsible for modulation of the conditional encounter rate Ψ(τ ) (fourth column). Combined with disease progression, this causes modulations of n inf and distributions P (n inf ) that vary with T lat and T inf (last two columns). Feature ii is responsible for zero-spreading events, which destabilize epidemic outbreaks and manifest in the peak of P (n inf ) for n inf ≈ 0 (last column). Feature iii is responsible for super-spreading events, which correspond to the tail of P (n inf ) and, thus, can cause a systematic shift of n inf .
FIG. 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 ρ(t) but fails to reproduce sufficiently heavy-tailed P (δt). A Weibull renewal process reproduces heavy-tailed P (δt) but fails to reproduce the cyclostationary ρ(t). Generative processes need to reproduce variability between individuals P (ntrain) to produce long-tailed distributions P (n 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.

D. 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 · 24 · 60 · 60s) and the last day (timestamp > 27 · 24 · 60 · 60s), 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.

E. 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 ( 5 m) and avoid to threshold the Bluetooth signal. Moreover, the temporal resolution of contacts is 20 s as opposed to 5 min 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 nondeterministic 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.  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 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 inf for fixed T 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 inf and thereby also the mean n 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 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 n 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.
In the main manuscript, we filtered the physical proximity data to only include contacts with a minimal Bluetooth signal strength (RSSI −80 dBm) as a proxy for contact distances below 2 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 (RSSI −75 dBm) and higher (RSSI −95 dBm) 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). 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 (RSSI −95 dBm) 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 (RSSI −75 dBm, 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.
In the main manuscript, we filtered contacts to only include those with a sufficient duration of at least 15 minutes. 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). 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 minutes vs. 15 minutes in the main manuscript) leads to more contacts qualifying and a higher conditional encounter rate. Vice versa, a longer required duration (e.g. 30 minutes) 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.