key: cord-0203457-6ksi5xiy
authors: Malz, Daniel; Pizzi, Andrea; Nunnenkamp, Andreas; Knolle, Johannes
title: Seasonal epidemic spreading on small-world networks: Biennial outbreaks and classical discrete time crystals
date: 2020-07-02
journal: nan
DOI: nan
sha: 0b362bc7e9278b56c2740ab7bc10e91c036f835a
doc_id: 203457
cord_uid: 6ksi5xiy

We study seasonal epidemic spreading in a susceptible-infected-removed-susceptible (SIRS) model on smallworld graphs. We derive a mean-field description that accurately captures the salient features of the model, most notably a phase transition between annual and biennial outbreaks. A numerical scaling analysis exhibits a diverging autocorrelation time in the thermodynamic limit, which confirms the presence of a classical discrete time crystalline phase. We derive the phase diagram of the model both from mean-field theory and from numerics. Our work offers new perspectives by demonstrating that small-worldness and non-Markovianity can stabilize a classical discrete time crystal, and by linking recent efforts to understand such dynamical phases of matter to the century-old problem of biennial epidemics.

Already its inventor, David Bernoulli, recognized the use of epidemic modelling to guide public health decisions by advocating inoculations to prevent the spread of smallpox [1] . Naturally, understanding and preventing disease has always been of great interest, reflected in the correspondingly vast body of literature [2] . Pioneered in the early 20 th century [3] [4] [5] , a modern formulation of epidemic dynamics uses coupled ordinary differential equations for the number of susceptible (S), exposed (E), infected (I), and recovered (R) individuals. Among the many rich dynamical phenomena that can be captured by S(E)I(RS) models is the emergence of a subharmonic response to a periodic variation in parameters, given for example through seasonality. In particular, biennial outbreaks have been observed in real-world measles case data more than a century ago [3] . SIR models to reproduce this phenomenon therefore have been an early target for computer-aided simulations [6] . Indeed, such a behavior is not specific to measles, but also appears in other diseases such as croup (cf. Fig. 1 ). To this day, sophisticated versions of these models, accounting for factors such as seasonality, immunity, cross-immunity, and the effect of distancing measures, are used to model outbreaks, including the COVID-19 pandemic [7] .

The elegant and simple description via global variables has had tremendous success, but misses potentially important effects due to spatial fluctuations and to the probabilistic nature of the transmission. This is more appropriately captured by probabilistic cellular automata on networks, where vertices (cells) represent individuals and edges potential infection pathways, which in recent years has shown great promise [8] . Among the key challenges is to understand the limits of the validity of global variables and to characterize phenomena that cannot be described otherwise. For example, the COVID-19 crisis has cast a spotlight on the shortcomings of using global variables to describe pandemics, which are both characterized by large spatial fluctuations and have fat tails [9] . A natural choice is to study infection dynamics on small-world [15] . A qualitatively similar behavior appears in our seasonal SIRS model on a small-world graph, in which the fraction of infected vertices pi also shows biennial major outbreaks (bottom, for φ = 5 × 10 −3 , k = 2 and N = 10 6 ). On an abstract level, this subharmonic response breaks time-translation symmetry, which can be considered a crystalline phase in time (bar with alternating +−, corresponding to major and minor outbreaks, respectively).

graphs [10] , as they exhibit the famous small-world property that characterizes social networks [11] , for which the average distance between any two vertices is several orders of magnitude smaller than the size of the network [12] . To give an example of using such a description, SI(R) models on small-world graphs naturally predict the well-known phenomenon of herd immunity [10] . Using percolation theory, one can show that in small-world graphs there is a critical level of immunity above which no extensive connected cluster of susceptible vertices exists, precluding outbreaks [10] . Furthermore, small-world networks have the added benefit that they allow for analytic understanding in the large-system-size limit [13, 14] . The graph is constructed starting from a periodic one-dimensional lattice in which each vertex is connected to its 2k nearest neighbors. In a 'rewiring' procedure, with probability φ, the ends of each short-range edge (dashed) can be moved to a random location in the graph, creating a shortcut. Each vertex is either susceptible (S, white), infected (I, black), or recovered (R, grey), that is, immune. Disease spreading occurs according to the shown three update rules, that are further detailed in the main text. (b) Phase diagram in the plane of the small-world parameter φ and winter duration Tw, comparing numerical results for N = 10 4 and k = 2 with the mean-field prediction. The trivial and DTC phases are characterized by annual or biennial outbreaks, respectively. We fit Tw = αφ ν to the first eight and last eight values (dashed) and find (α, ν) = (0.12 ± .01, −1 ± .01) and (.65 ± .02, −.65 ± .01). The mean-field prediction (solid red) shows striking agreement with numerics for small φ, but cannot explain the deviation at large φ, which we attribute to a number of effects that are missed by mean-field theory, most notably the discreteness of space and time.

transient phenomena or can rather persist up to infinite times is a very interesting and highly non-trivial one. The subharmonic response observed in structureless mean-field models is generally spoiled when moving to a network, where noise tends to randomize the phase of the response and to ultimately destroy any system's autocorrelation with its far past. Against this expectation, a 'classical discrete time crystal' (DTC) is a dynamical phase of matter that maintains a persistent subharmonic response in the thermodynamic limit, in spite of noise and perturbations. DTCs have originally been studied in the context of quantum many-body localization [16] [17] [18] [19] and have in the last few years attracted a tremendous amount of interest [20] [21] [22] [23] [24] [25] . Recent work has extended the field to classical stochastic Markovian processes in one [26, 27] and two dimensions [28] .

Here, we investigate a seasonal SIRS model on small-world graphs and show under which circumstances it gives rise to biennial outbreaks, supporting numerics with an analytical meanfield theory that captures the salient features of the model's rich phenomenology. With a scaling analysis, we show that in the thermodynamic limit, biennial outbreaks exhibit longrange order in time, heralded by a diverging autocorrelation time. We attribute the robustness of the DTC to the interplay of small-worldness and, crucially, non-Markovianity. Our work links the well-known phenomena of biennial outbreaks to the theory of classical DTCs and thus unveils striking real-world ramifications of these remarkable dynamical phases of matter.

The remainder of this paper is organized as follows. In Section II we introduce the model, detailing the structure of the small-world network and the update rules describing seasonal epidemic spreading, and present its dynamical phase diagram. In Section III we gain understanding of the system by deriving an analytical mean-field theory. A thorough characterization of the model is presented in Section IV, whereas the timecrystalline nature of the biennial epidemics is investigated in Section V by means of a scaling analysis to assess its stability. Finally, in Section VI we discuss the results and conclude with an outlook for future research.

We consider a stochastic microscopic SIRS model for seasonal epidemic spreading on the celebrated small-world networks introduced by Watts and Strogatz [12] , sketched and described in Fig. 2a . The shortcut density φ parameterizes the 'small-worldness', which allows us to interpolate in a controlled manner from a linear 1d chain (φ = 0) to random graphs (φ = 1). Each vertex can be in one of three states: susceptible (S), infected (I), immune/recovered (R), and the system evolves according to the following rules:

(I) Autumn: vertices are exposed to the disease with a probability q and, if susceptible, contract it, S→I.

(II) Winter: in each of total T w steps, the infection spreads from infected vertices to susceptible vertices they are connected to, S→I.

(III) Summer: the infected vertices recover and become immune, I→R, whereas immune vertices lose their immunity and become susceptible again, R→S.

These rules are repeated annually, which defines the fundamental periodicity of the model. In the first year we take all vertices to be susceptible, but the qualitative results shown here are robust to perturbations of the initial condition. Underlying this model is the idea that, in a social network, long-distance travellers act as shortcuts for disease spreading across distant regions of the globe, as parameterized by the small-world parameter φ. Similarly, we can imagine that in autumn the vertices are exposed to the infection because of travellers returning from hidden regions (e.g., faraway countries) that are not accounted for explicitly in the network, so that q ∼ φ. In the following, we make the choice q = φ for concreteness. We stress that our goal here is not to model real-world data, but instead to study a minimal toy model to explore the statistical mechanics of epidemic spreading.

Depending on the value of the parameters T w and φ, at long times the system can behave in two fundamentally different ways, featuring outbreaks that occur either annually or biennially. As we will show later, these two behaviors constitute two genuine dynamical phases of matter: The phase with annual outbreaks preserves the time translation symmetry of the model, and can therefore be referred to as the 'trivial' phase. In contrast, if outbreaks occur biennially, the discrete time translation symmetry is broken, and the system enters the DTC phase. The phase diagram, derived analytically and confirmed numerically, is shown in Fig. 1d . Intuitively, the transition may be understood in terms of the effect of immune vertices. After a long winter (T w large) and/or a fast outbreak (φ large), in which a large proportion of vertices has been infected, the established herd immunity suppresses the outbreak in the following year. Two years later, the immunity has decayed, and a large outbreak occurs again. At the core of these effects is clearly the immunity, which prevents a given vertex to be infected for two consecutive years, a memory effect that we refer to as 'non-Markovianity'. Since the infection spreading is quicker on a smaller world, the critical φ c dividing the two phases decreases with the winter duration T w . We fit αφ ν to the data (shown as dashed lines in Fig. 2b and obtain the critical exponents ν = −1 ± .01 for small φ, and ν = .65 ± .01 for large φ. At small φ the agreement with mean-field theory is striking, whereas for large φ the discreteness of the model plays an increasing role, which cannot be accounted for by mean field.

To gain analytical understanding of the trivial phase, the DTC, and the transition between them, we derive a meanfield theory. To do so, we first find analytical expressions for the epidemic dynamics throughout a single winter, and then describe how the epidemic evolves from one year to the next.

Pioneered by Newman et al. [13] , an accurate mean-field description can be found in the limit of large graphs N → ∞ and small small-world parameter φ, by treating space and time as continuous, such that the epidemic dynamics can be described by ordinary differential equations. In contrast to the approach in Ref. [13] , here we have to account for the presence of immune vertices that act to suppress outbreaks, as this underpins the biennial response. We do so while retaining the analytic solubility of the model by making the approximation that immune vertices are randomly distributed. Thus, they can be accounted for by the probability P S|Ī that a non-infected vertex is susceptible (rather than immune)

where p i (t) and p r are the fractions of infected and immune vertices, respectively. Since immunity only changes in summer (cf. Fig. 2a) , p r is constant throughout the winter.

In the following we mean by 'infected region' a contiguous local chain (connected by short-range edges) of infected vertices. In contrast to the concept of connected clusters, we therefore count two regions connected by a long-range edge as two separate regions. For instance, there are 3 regions in Fig. 2a . Let us denote the density of infected regions by ν i (t). In every time step each of the 2N ν i infection fronts (boundaries between infected and susceptible regions) advances by k steps, enveloping 2N kν i vertices. To obtain the rate at which sites are infected, we have to multiply this by the probability that the enveloped vertices are susceptible,

In the spreading of infection fronts, the number of infected regions may change according to two mechanisms. One, a new infected region is spawned if an infection front crosses a shortcut between two susceptible vertices. The rate for this happening is the product of the vertices exposed in one times step, 2N kν i , and the density of shortcuts among non-infected regions, 2kφ(1 − p i ). The factor of 1 − p i arises, as the remaining shortcuts necessarily connect two uninfected regions. Two, infected regions can merge, which decreases ν i at a rate of 2kν 2 /(1 − p i ). This factor is easiest understood as the product of the rate at which the total susceptible regions shrink (2kν) times the density of boundaries to infected regions in the remaining susceptible fraction (ν/(1 − p i )). Optionally, the distribution of gap sizes can be considered explicitly [13] . Crucially, both these processes will be effective only if the involved non-infected vertices are susceptible, which happens with probability P S|Ī (t), such that

Eqs. (2) and (3) can be integrated analytically (cf. SI), yielding

where the notation f (p r , t) emphasizes the dependence on the immune fraction p r , and where the inverse timescale θ reads

As a sanity check, in the limit t → 0 we obtain p i = φ(1 − p r ), recovering the initial fraction of infected vertices, whereas in the limit t → ∞ the fraction of infected vertices asymptotically approaches p i → 1 − p r .

The solution to epidemic spreading in a single year provides us with a formula for the total fraction of infected vertices at the end of one year p

r is the immune fraction during the n-th winter. The summer transition from infected to immune means that p

, and therefore we obtain the discrete 'one-year map'

where we have omitted the dependence of f on the winter duration T w , which is henceforth treated as a fixed parameter.

Applying f twice, we obtain the two-year map f (2) p

which is helpful in describing the biennial epidemics. Note that to leading order in φ, only the product of T w and φ appears in f , which readily explains why at small φ, the critical small-world parameter scales as φ c ∼ 1/T w (cf. Fig. 2(b) ).

To conclude this section, we briefly discuss the shortcomings of our mean-field theory. Most importantly, while the variable ν i contains the information on the number of infected regions, no such information is stored for the immune population, which only enters in the form of the probability P S|Ī . This discounts the effect that clusters of immune sites may stop infection fronts, and also enclose a susceptible cluster, preventing it from being infected. Indeed, the mean-field model does not account for any fluctuations or other random structure that may emerge, such as fluctuations in the number of initially infected vertices, the relative positioning that may or may not be conducive to fast outbreaks, the number of long-range links, the size of the connected cluster of the randomly sampled small-world graph, etc. Ultimately, the mean-field theory is validated by the good agreement with numerics, which is shown in the next Section.

Throughout a single winter, the fraction of infected vertices p i (t) follows a characteristic logistic-like growth, with an initial superlinear growth followed by saturation to an asymptotic value at long times. In Fig. 3(a) we verify that the mean-field agrees remarkably well with the full numerics in the absence of immune vertices (p r = 0). The multi-year dynamics is investigated in Fig. 3(b) , which, depending on the small-world parameter φ, exhibits two qualitatively different behaviors. If φ is smaller than a critical φ c , we observe annual outbreaks, whereas if φ > φ c , major outbreaks occur biennially, and p i (t) oscillates with a period of two years.

We observe that the accuracy of the mean-field predictions deteriorates after the first winter. To understand why this is the case, we note that a vertex that gets randomly infected at the beginning of the winter typically lies in the bulk of a susceptible region. In the early stages of the winter the disease thus propagates essentially unperturbed by the presence of the immune vertices. Only at a later stage of the winter some infection fronts get stopped by immune regions. This effect depends on the spatial structure and distribution of immune vertices and therefore cannot be captured by mean field. Nonetheless, as we are about to see, mean-field theory is still remarkably accurate in capturing and explaining the phase transition between annual and biennial epidemics.

Within mean-field, in fact, the long-time dynamics can be understood using the tools of discrete dynamical maps [29] . A simple inspection as shown in Fig. 3 (c) reveals that the equation x = f (x) is always fulfilled by one and just one fixed point (FP), that we call x 2 . The stability of this FP is determined by the first derivative f (x 2 ). Since the infected fraction decreases with increasing immune fraction, f (x) < 0, and the only possibilities are

When x 2 is stable, p (n) n→∞ − −−− → x 2 , and the system lies in the trivial phase with annual outbreaks at long times. To understand what happens when x 2 is unstable, we turn to the two-year map f (2) in Eq. (7), which is also investigated in Fig. 3(c) . As it is easy to check, x 2 is a FP also of f (2) . When the FP x 2 is unstable, f (2) (x 2 ) = [f (x 2 )] 2 > 1, indicating the emergence of two new stable FPs x 1 and x 3 for f (2) . When x 2 is unstable, at long times we therefore get p Plotting the FPs of the two-year map f (2) versus the smallworld parameter φ, in Fig. 3(d) we obtain a bifurcation diagram. At T w = 30, the critical small-world parameter that separates the trivial and the DTC phases is φ c = 4.24(8) × 10 −3 . Near criticality, the mean-field scaling x 3 − x 1 ∼ (φ − φ c ) ξ is characterized by the critical exponent ξ = 1/2 (see inset). Numerics suggests that the actual critical exponent may be larger, although large fluctuations close to the transition, a finite-size effect, prevent us from drawing clear conclusions.

Varying the winter duration T w , this method allows us to calculate the full phase diagram shown in Fig. 2b . The mean-field model predicts the phase separation to lie at φ c = 0.126918 T −1 w (dashed line), which shows striking agreement with numerics for small φ, which predicts φ c = 0.12065T −1 w . As in previous work [13] , the approximations made to derive the mean-field solution deteriorate at larger values of φ, which explains the deviation. Considering no immune vertices in the network, we observe a logistic-like growth, with an initial superlinear rise (quickest for smaller worlds, that is larger φ) and eventual saturation when all vertices get infected, as it is also accurately reproduced by mean-field. (b) Epidemics over multiple years can have two qualitatively very different behaviors. For φ < φc (left) and φ > φc (right), outbreaks occur annually and biennially, respectively. In the latter case, corresponding to a small-enough-world, the system responds with a period twice that of the underlying seasonality, pointing to a DTC. Mean-field (solid line) is compared to numerics for N = 10 4 (markers). (c) Mean-field fixed-points analysis. The one-year map f only has one FP x2. If φ < φc, x2 is stable and the end-of-winter infected fraction p (n) i eventually reaches x2. On the contrary, if φ > φc, x2 becomes unstable and p (n) i eventually oscillates at every year between x1 and x3, the emerging stable FPs of the two-year map f (2) . The mean-field dynamics is highlighted with a cobweb in blue. (d) Bifurcation diagram of the FPs of f (2) within mean-field (red lines) and end-of-winter infected fraction p (n) i in the even and odd years at long times for N = 10 5 (blue markers). Solid and dashed lines stand for stable and unstable FPs, respectively. In mean-field, the DTC order parameter x3 − x2 scales with critical exponent 1/2 (inset). Note, in (b) and (d) the mean-field prediction is in this case qualitatively correct but quantitatively inaccurate, due to its inability to account for the clusterized structure of the immune vertices. Mean-field is surprisingly successful in locating the critical φc.

It is impossible to determine whether the biennial response is a genuine dynamical phase of matter from mean-field theory alone. The reason is that microscopic fluctuations typically destroy correlations over time, due to the presence of random phase flips (defects in the crystal). A defect in this setting means that two consecutive years have minor (or major) outbreaks, such that the biennial response is broken.

To signal a DTC, the subharmonic response of the system needs to show a diverging autocorrelation time with system size, and this divergence should be robust both to noise and to the perturbation of the model's parameters [26, 28] . To analyze this, we extract the defect density from calculated time traces and observe their scaling with system size (cf. Fig. 4 inset). In the trivial phase, we find that defects occur at random irrespective of system size. This is expected, as there is no subharmonic response in the first place. With increasing winter duration T w , the system crosses into the DTC phase and we observe clear exponential scaling with system size (see inset in Fig. 4) . The exponent λ(T w ) characterizing this scaling shows a clear departure from zero as the system enters the DTC phase, e.g., at T w ≈ 13 in Fig. 4 .

Revisiting the century-old problem of biennial epidemics from the perspective of statistical mechanics, we have studied a microscopic model of seasonal epidemic spreading on a smallworld network, giving rise under appropriate circumstances to a biennial outbreak behavior, which is well-described by a fully analytical mean-field theory. Surprisingly, the subharmonic response of the stochastic dynamical system is stable in the thermodynamic limit against noise and perturbations, and we Previous research highlighted that classical DTCs are typically prevented by deconfined domain walls in short-range 1d models [26] , but can be stabilized by either larger dimensionality [30] or long-range interactions [27] . These works also show how the essential features of classical DTCs can be captured by probabilistic cellular automata (PCA) [26, 27] , in which the system is described by two-state variables. This raises the more general question which aspects (dimensionality, range of interactions, contractive dynamics,..) of these models is responible for stabilizing the time crystalline phase. In this context, on the one hand our work explores small-world graphs, in which a few random links shortcut across a 1d chain, accounting for the small-world property of real-world social networks. On the other hand, a crucial role is played by immunity, which could equivalently be understood as a non-Markovian memory effect that precludes a vertex from being infected two years in a row. Our work therefore identifies non-Markovianity and small-worldness as two key ingredients that can stabilize timecrystallinity against the proliferation of noise-induced defects.

We conclude outlining some possible directions for future research. First, it would be desirable to improve the mean-field theory further such that it accounts explicitly for the number ν r of immune regions. Such a theory could start from differential equations analogous to Eqs. (2) and (3), but introducing the parameter ν r would most likely come at the price of losing analytical solubility. Second, the model studied here is very simplistic compared to models used to forecast actual epidemics. Additional ingredients may reveal a richer phase diagram or novel phenomena. As a few examples, one could for instance consider the survival of the infection through the summer (rather than a random autumn exposure), account for different epidemic spreading rules in various regions of the network (e.g., as to mimic different seasonalities in the two hemispheres), include more states of the vertices (e.g., infected asymptomatic) or strains of virus, or account for randomized and/or longer immunity times, which could result in the emergence of n-DTCs with major outbreaks every n years. Interestingly, the last point relates to early work that argued that only n = 2 could produce a stable subharmonic response [31] , in contrast to what our work would suggest. A particularly timely avenue is to investigate the effects of periodic interventions (such as an intermittent reduction of φ as to mimic travel restrictions), and whether these can lead to a subharmonic response, too. Third, more complex, multistate variable models support even richer dynamics, such as the prime-numbered response in predator-prey models on 2d graphs [32] . While in our model the subharmonicity is set by the duration of immunity, this raises the possibility of time crystals in which the periodicity is generated dynamically.

Finally, an open question that persists is whether models of the kind studied here may support continuous time symmetry breaking. Indeed, in the absence of periodic driving, this model is known to produce a periodic response [33] , and future investigations should clarify whether the system's autocorrelation time diverges with the system size. The possibility of continuous time symmetry breaking, proven impossible for ground states of time-independent quantum Hamiltonian systems [34] , opens a new avenue of research in a classical stochastic setting.

Supplementary Information for "Seasonal epidemic spreading on small-world networks: Biennial outbreaks and classical discrete time crystals" Daniel Malz, Andrea Pizzi, Andreas Nunnenkamp, and Johannes Knolle These Supplementary Information are devoted to technical derivations and complimentary results and consist of two sections, S1 and S2, in which we provide details on the mean-field theory and on the numerical routines, respectively. 

To solve Eqs. (S1) and (S2), we take their ratio, reading

which is solved by

The constant C is set by initial conditions. At the beginning of the winter, a fraction q of the population is exposed to the disease and, being a fraction µ r of the population susceptible, we have µ i (0) = 1 − qµ r . The initial density infected regions ν i (0) is instead computed as follows. A given vertex will be an edge between an infected region and a non-infected one with probability 2qµ r (1 − qµ r ). Because each continuous non-infected region is associated with two such edges, we get ν i (0) = qµ r (1 − qµ r ). The constant C therefore reads

Plugging Eq. (S4) back into Eq. (S1) we get

from which, after completing the square and integrating both sides, we get

where we defined the positive constants

Solving the integral in Eq. (S7), we get

where x(t) = 2µi(t)−α γ . Inverting Eq. (S9) for x(t) we obtain

from which finally we obtain the fraction of infected people p i (t) = 1 − µ i (t) for a given fraction of immune people p r = 1 − µ r throughout a winter

After straightforward manipulations, Eq. (S11) is rearranged as

(S13)

In the case q = φ, Eq. (S13) recovers Eq. (4) in the main text.

Here we describe the various numerical routines used to calculate the data shown in the main text.

To calculate the infected fraction as a function of time, as for example displayed in Fig. 3a , we follow these steps:

1. Randomly generate a small-world graph [12] . This is done by taking a periodic chain of N vertices, each connected to its 1 st , 2 nd , · · · , k th neighbors (resulting in a coordination number 2k), and moving each edge end to a random location with a probability φ. For the results shown in the main text, we take k = 2 everywhere and vary φ and N according to the description in the figures' captions.

2. To initialize the dynamics, each of the N vertices is infected with a probability q (Autumn in Fig. 1c) . The parameter q is in principle independent, but, since the autumn infections can for example be thought of as arising from long-range connections to other parts of the world, it is meaningful to consider q ∼ φ. For concreteness, in the main text we considered q = φ.

3. The infection dynamics is run for T w time steps, which yields a certain final state. The fraction of infected vertices at each point in time is denoted by p i (t), and this is what we show for instance in Fig. 3a,b. 4. In "summer", all the vertices that were in the immune state are put back in the susceptible state, whereas the infected vertices are transferred to the immune state, such that they do not participate in the infection dynamics of the following year. Points 2-4 are repeated for every year, generating time traces such as in Fig. 3b .

To each year in a multi-time trace we can associate the total number of infected vertices in that year, which measures the size of the outbreak. Corresponding sample data are shown in Fig. S1a . Depending on system parameters, we observe that these values cluster either around a single or two values, which is easily appreciated by looking at histograms of the data (Fig. S1b) . We extract the dominant peaks from the histograms and by majority vote over many averages determine whether a given set of parameters lies in the trivial or DTC phase (Fig. S1c) , which we assemble to form the phase diagram (Fig. S1d) . 

To obtain a bifurcation diagram, we again take the data of a large number of time traces, but now fit either a single or two Gaussians to the total distribution of all infected fractions (Fig. S1e) . The obtained mean and variance are represented by a point and error bars in the bifurcation diagram (Fig. S1f) .

We again take all time traces corresponding to a given set of parameters {N, φ, T w } and now calculate the differences in the total infected fraction of adjacent years. In a DTC, the differences should alternate between positive and negative differences. We count the number of occurrences of differences of the same sign next to each other (Fig. S1g) and divide by the total number of years to obtain the defect density. This is calculated for N = {5, 6, · · · , 30} × 10 2 which reveals how the defect density changes

Essai d'une nouvelle analyse de la mortalite causee par la petite verole

Infectious Diseases of Humans

On epidemic disease in England-The evidence of variability and persistency of type

A contribution to the mathematical theory of epidemics

The Interpretation of Periodicity in Disease Prevalence

Recurrent outbreaks of measles, chickenpox, and mumps

Projecting the transmission dynamics of SARS-CoV-2 through the postpandemic period

Epidemic processes in complex networks

Tail risk of contagious diseases

Epidemics and percolation in small-world networks

An Experimental Study of the Small World Problem

Collective dynamics of 'smallworld' networks

Mean-Field Solution of the Small-World Network Model

Mean-Field and Anomalous Behavior on a Small-World Network

Anderson, Pediatric hospitalizations for croup (laryngotracheobronchitis): biennial increases associated with human parainfluenza virus 1 epidemics

Phase Structure of Driven Quantum Systems

Discrete Time Crystals: Rigidity, Criticality, and Realizations

Equilibration and order in quantum Floquet matter

Modeling spontaneous breaking of time-translation symmetry

Observation of discrete time-crystalline order in a disordered dipolar many-body system

Observation of a discrete time crystal

Discrete time-crystalline order in cavity and circuit qed systems

Period-n discrete time crystals and quasicrystals with ultracold bosons

Long-Range Prethermal Phases of Nonequilibrium Matter

Classical discrete time crystals

Bistability and time crystals in long-ranged directed percolation

Classical stochastic discrete time crystals

Complex and Adaptive Dynamical Systems

Classical stochastic discrete time crystals

Stability of temporally periodic states of classical many-body systems

Prime number selection of cycles in a predator-prey model

Small World Effect in an Epidemiological Model

Absence of Quantum Time Crystals

We observe clear exponential behavior, to which we fit exponentials. The fitted exponents and errors are shown in a second diagram as function of φ for a specific value of T w (Fig. S1i)

Acknowledgements. We acknowledge support from the Imperial-TUM flagship partnership. D. M. acknowledges funding from ERC Advanced Grant QENOCOBA under the EU Horizon 2020 program (Grant Agreement No.