key: cord-0610645-02gsos1a
authors: Pausch, Johannes; Garcia-Millan, Rosalba; Pruessner, Gunnar
title: Noise can lead to exponential epidemic spreading despite $R_0$ below one
date: 2021-09-01
journal: nan
DOI: nan
sha: c41ec68edfe17a28b7af717fffc5303df69bd64a
doc_id: 610645
cord_uid: 02gsos1a

Branching processes are widely used to model evolutionary and population dynamics as well as the spread of infectious diseases. To characterize the dynamics of their growth or spread, the basic reproduction number $R_0$ has received considerable attention. In the context of infectious diseases, it is usually defined as the expected number of secondary cases produced by an infectious case in a completely susceptible population. Typically $R_0>1$ indicates that an outbreak is expected to continue and to grow exponentially, while $R_0<1$ usually indicates that an outbreak is expected to terminate after some time. In this work, we show that fluctuations of the dynamics in time can lead to a continuation of outbreaks even when the expected number of secondary cases from a single case is below $1$. Such fluctuations are usually neglected in modelling of infectious diseases by a set of ordinary differential equations, such as the classic SIR model. We showcase three examples: 1) extinction following an Ornstein-Uhlenbeck process, 2) extinction switching randomly between two values and 3) mixing of two populations with different $R_0$ values. We corroborate our analytical findings with computer simulations.

Branching processes are widely used to model evolutionary and population dynamics as well as the spread of infectious diseases. To characterize the dynamics of their growth or spread, the basic reproduction number R0 has received considerable attention. In the context of infectious diseases, it is usually defined as the expected number of secondary cases produced by an infectious case in a completely susceptible population. Typically R0 > 1 indicates that an outbreak is expected to continue and to grow exponentially, while R0 < 1 usually indicates that an outbreak is expected to terminate after some time.

In this work, we show that fluctuations of the dynamics in time can lead to a continuation of outbreaks even when the expected number of secondary cases from a single case is below 1. Such fluctuations are usually neglected in modelling of infectious diseases by a set of ordinary differential equations, such as the classic SIR model. We showcase three examples: 1) extinction following an Ornstein-Uhlenbeck process, 2) extinction switching randomly between two values and 3) mixing of two populations with different R0 values. We corroborate our analytical findings with computer simulations.

A branching process [1, 2] is a stochastic process in which individuals randomly create copies of themselves or become extinct. This model has a wide range of applications and individuals might represent offspring [3] [4] [5] , particles [6] [7] [8] , active neurons [9] [10] [11] [12] [13] , or infected individuals [9, [14] [15] [16] among others [17] . At the centre of many investigations of branching processes is the statistics of spells of activity, which are called avalanches [6] [7] [8] [9] [10] [11] [12] [13] in most applications or outbreaks in the context of infectious diseases. Avalanches, or outbreaks, are defined as the activity in the branching process which is initiated by a single individual and lasts until all subsequent individuals have become extinct. In the context of particles, an avalanche starts with a single particle and ends when the system is empty. In the context of infectious diseases, an outbreak starts with a single infected individual and ends when no infected individual is left (here we use the words infected and infectious without distinction).

A branching process is subcritical if the expected size of an outbreak decays exponentially. It is called supercritical if the expected size of the outbreak grows exponentially. Otherwise, if the size approaches a constant value in time, the branching process is said to be critical. We will use this criterion, asymptotically constant expected outbreak size, as the definition of the critical point throughout this work. The critical point generally divides a parameter region resulting in asymptotically exponential growth from one resulting in asymptotically exponential decline. For particle avalanches as well as for outbreaks of infectious diseases, it is of great interest to identify simple parameters that indicate whether the process is supercritical or not. One of those parameters is the basic reproduction number R 0 .

The basic reproduction number R 0 is defined as the expected number of secondary individuals that are created from a single individual [18] [19] [20] [21] [22] . More explicitly in a branching process, an existing individual waits until a branching event occurs. In many models, the waiting time between branching events is fixed [1, 23, 24] , however in the present work, we consider only exponentially distributed waiting times [8, 11, 25] . At a branching event, an individual is replaced by its offspring, which is a random number K ∈ N 0 of individuals. The case K = 0 corresponds to the extinction of the parent individual. In general, the offspring distribution can be defined by specifying the branching probabilities p 0 , p 1 , p 2 , . . . ∈ [0, 1] such that

i.e. the probability P (K = k) that an individual has K = k offspring equals p k . Although our analytical results hold for any distribution of K, in our simulations we use only the binary offspring distribution with K ∈ {0, 2}. All individuals are independent and there is no bound to the number of individuals in the system, which can be interpreted as an unbounded population of susceptible individuals. In this setup, the basic reproduction number R 0 is defined as the expected number of offspring,

In particular, R 0 is dimensionless and does not indicate how quickly or slowly an avalanche/outbreak evolves. There are many difficulties in deriving R 0 from data [21, 26, 27] and researchers have defined several similar arXiv:2109.00437v1 [q-bio.PE] 1 Sep 2021 quantities related to R 0 [28, 29] . In addition, more detailed models of infectious diseases take other characteristics such as age, immunity, behaviour or the evolution of the disease itself into account, which make the definition of R 0 more difficult [21, [26] [27] [28] 30] . Rather than including more detailed aspects into our models, we restrict ourselves to the basic model outlined above and keep the discussion at the level of stochastic processes. For example, we assume that infected individuals are also infectious.

In many real-world occurrences of branching processes, the environment and the process itself are imperfect in the sense that they fluctuate in time, for example because individuals and the environment change the conditions for disease transmission [31] . Such fluctuations will be affecting the branching process over time and are not easily dealt with analytically. Of particular importance are fluctuations that affect the population as a whole. As our calculations show, basic approximations of branching processes with noise can be misleading by predicting subcritical dynamics where a more detailed calculation reveals supercritical behaviour. It is the main aim of this article to highlight such, often counter-intuitive, phenomena, which are easily missed by traditional modelling of epidemics, which draw on a coarse-grained set of equations, such as the classic SIR model [32] .

The article is organized as follows: In Sec. II, the branching model without noise is presented as a Master Equation. It forms the basis of the models with noise that follow. In Sec. III, we introduce a branching process coupled to an Ornstein-Uhlenbeck process. Although a mean-field approximation predicts a critical point at R 0 = 1, our detailed analysis supported by simulations reveal a shift of the critical R 0 to values smaller than 1. In Sec. IV, the branching process is coupled to a stochastic process called telegraphic noise, which implements a random switch between two different extinction rates of the branching process. In this system, it is much less obvious what a suitable definition of R 0 would be. Neither of the two R 0 values associated with the two extinction rates, nor a simple weighted average of them predict R 0 = 1 to be the critical point. An exact calculation, reveals that the critical R 0 is smaller than 1. In Sec. V, two branching processes with two different R 0 are coupled. We show that neither of the two R 0 nor a linear combination of them correctly predict the critical point of the system. We conclude in Sec. VI. The detailed analytical results are based on field-theoretic approaches which are presented in the appendix.

At the center of our models is the Master Equation for a branching process [8] . In order to use a consistent language, we refer to individuals being created or becoming infected in branching events. They are then present in the population. We will say that they disappear from the population when an extinction occurs. We hope that the reader can translate this language to their application by replacing individual with particle, signal, ... and population with system, neural network, ... according to their requirements. Let the population size N (t) denote the number of individuals present at time t, with initial condition N (0) = 1, and P (N (t) = n) the probability that there are n individuals present at time t. Then, the probabilistic dynamics of the basic branching process are described by the master equation [8] 

where p k is the offspring distribution and s is the overall event rate. In particular, there are two probabilistic components: 1) s is the parameter for an exponentially distributed waiting time until a branching event occurs, and 2) at a branching event of an individual, the offspring distribution p k determines by how many individuals it is replaced, i.e. the original individual does not continue to exist alongside the k new individuals. Alternatively, one could say that k − 1 new individuals are created while the original individual continues to be present. The independence of the processes implies that the product sp k equals the event rate for an individual to be replaced by k other individuals. In the case k = 0, we denote by := sp 0 the extinction rate, i.e. the rate of the exponential distribution that determines the time when an individual disappears from the population without producing any offspring. In the present work, we consider branching dynamics for which fluctuates in time.

To clarify the rôle of the extinction rate , we rewrite Eq. (3) as

In our approach, there is no modelling of a healthy population and there is no saturation of a population with infected individuals. In particular, there is no upper bound to the number of infected individuals. This is because we want to identify the effect of the noise on the extinction rate without getting lost in the too many details and model parameters.

The time-homogeneous branching process described by the Master Equation ( 3) has been studied before [1, 8] . The temporal evolution of the expected number of infected individuals,

follows from (3),

with R 0 as defined in Eq. (2) . The solution E [N (t)] = exp(s(R 0 − 1)t) illustrates the important role of the basic reproduction number R 0 : if R 0 < 1, the expected number of individuals decreases over time (subcritical case), while it increases exponentially if R 0 > 1 (supercritical case). The case R 0 = 1 is called the critical point as the expected number of individuals stays constant in time.

The model described by Eq. (3) does not rely on the law of large numbers, i.e. it is valid for small numbers and even for one individual. The equation also allows deriving non-trivial dynamics for higher moments. For example, the equation that governs the variance V[N (t)] of the number of infected individuals is

Here the appearance of E (K − 1) 2 implies that the second moment of the offspring distribution affects the dynamics too -not just R 0 . The model in Eq. (3), and therefore also Eqs. (4), (6) and (7), assumes a static setup, i.e. branching always happens with the same rates and the same offspring distribution. However, this assumption may be unrealistic for many applications. What if event rates and offspring distributions fluctuate in time? This is the topic of the next two sections, after which we also consider the case where two populations with different basic reproduction numbers R 0 interact.

As a first example of noisy branching processes, we couple the extinction rate in Eq. (4) to an Ornstein-Uhlenbeck (OU) process

where the rate y(t) is governed by

Here, η(t) is a Gaussian white noise with mean E [η(t)] = 0 and correlator E [η(t)η(t )] = 2Dδ(t − t ). The dimensionless parameter λ is the coupling strength and β is the return rate. The persistence time β −1 induces temporal correlations, or a memory, in the noise y(t) in a similar way to active fluctuations in the motion of active Ornstein-Uhlenbeck particles [33, 34] , Eq. (18) . In principle, (t) may become negative, which would render the process ill-defined. In numerical simulations, we can guard against that by monitoring the value of and replacing it by 0 whenever it becomes negative. For the parameters considered below, this is exceedingly rare, affecting a single realisation in well over 10 10 . An example trajectory is shown in Fig. 1 . Eq. (9) implies that the steady state distribution of y is the Gaussian distribution

and therefore the time average equals E τ [y] = 0. The noisy branching process Eq. (4) with Eq. (8) is described by the following master equation:

The additional contribution y(t) to the extinction rate (t) is affecting all individuals equally, so that rather than being reduced by the law of large numbers, the effect of y(t) grows linearly with the population size.

What is the effect of this noise on the dynamics of the branching process? A mean-field approximation predicts that this perturbation does not have any impact, because in the mean-field approach all occurrences of y(t) are replaced by its mean E τ [y] = 0 and therefore the mean-field expected number of secondary cases equals R 0 , Eq. (2). In other words, mean-field theory predicts that R 0 = kp k = 1 results in a critical process.

However, closer inspection reveals that with OU noise, the offspring distribution P (K = k) = p k , Eq. (1), has effectively become time-dependent. To see this, we regard the branching process as a collection of simultaneous, independent, exponentially-distributed waiting processes -one process for each K ∈ N 0 . Because they are independent, the waiting time until the first of these events occurs is exponentially distributed with rate s + λy(t). Which of the processes actually occurs first can be answered probabilistically by calculating the ratio of the rate of that process divided by the sum of the rates of all simultaneous processes:

which satisfies normalisation, k≥0 P (K = k) = 1, and describes the effective offspring distribution given a rate y = y(t). Hence, the time-dependent expected number of secondary cases equals

with R 0 = kp k , Eq. (2). The effect of y(t) in Eq. (13) is not symmetric about 0, as can be seen by expanding s/(s+λy) = 1−(λy/s)+(λy/s) 2 −. . . Taking the expectation over y, suggests E [K(t)] = R 0 (1+(λ/s) 2 E y 2 +. . .).

over the stationary distribution of y, Eq. (10), effectively the time-average of E [K(t)], can be calculated in closed form

where D + denotes the Dawson function, defined in (B7). Demanding the expected offspring number E y [E [K(t)]] to be unity at the critical point, produces

which generally is less than unity, as can also be gleaned from R 0 = 1 − (λ/s) 2 D/β + . . . using the expansion discussed after Eq. (13) and E y 2 = D/β from Eq. (10). However, (15) needs further scrutiny as it is based on the wrong assumption, as we will explain now. The process is at the critical point when the average number of offspring spawned per reproductive event is unity. How- (15) is an average over the stationary distribution of y, assuming that the same number of reproductive events take place at any such value. Because of the temporal correlations in y(t), however, episodes of high extinction typically occur when the population size is small anyway. As a result, fewer reproductive events are affected by high extinction rates than by low extinction rates.

Using a Doi-Peliti field theory, which is derived and explained in Appendix B 1, we can calculate the expected population size directly, producing the final result

The basic reproduction number R 0 that results in asymptotically constant expected population size, ∞ > lim t→∞ E [N (t)] > 0, according to (16) , is the one that makes the exponent vanish for all t, namely 

of a branching process driven by an Ornstein-Uhlenbeck noise estimated from numerical simulations (symbols) based on 10 7 iterations. The parameters of the noise in each panel are as follows: λ = 0.05 (top row); λ = 0.1 (bottom row); β = 0.5, D = 0.125 (left column); β = 1, D = 1 (middle column); β = 1, D = 2 (right column). For each set of noise parameters, we can choose different values of R0 such that the population has supercritical (red symbols, black lines), critical (orange symbols, grey lines), and subcritical dynamics (blue symbols, pink lines). Our numerical results are in agreement with Eq. (16) . In particular, we find values R0 < 1 such that the population displays supercritical behaviour due to the external noise.

Remarkably, in any non-trivial setup, this critical R 0 is less than unity. The average number of secondary infections of an isolated individual not subject to noise needed to sustain an outbreak, is thus less than unity. The explanation for this counter-intuitive result is similar to the reason why Eq. (15) is based on the wrong assumptions: Because the noise is correlated, in general population sizes experiencing low extinction rates are larger than those experiencing large extinction rates. The noise correlator of the the Ornstein-Uhlenbeck process (8) is

whose characteristic time β −1 is a measure of the persistence of active fluctuations. Although the noise has vanishing mean, its effect is biased towards larger population sizes. In addition, Eq. (14) does not incorporate the change in frequency with which events take place overall -at times of high extinction rates, more events take place than at times of low extinction rates. Our field theoretic result Eq. (17a) provides a systematic expansion of the critical R 0 in orders of λ and is subtly different from the ad hoc result (15), as the denominator of the leading order correction in (15) is βs 2 rather than β 2 s in Eq. (17a).

To test Eq. (17a) numerically, we have performed Monte-Carlo simulations to estimate the critical R 0 as shown in Fig. 5 . Given the other parameter values, a fairly small range of λ is available, as otherwise (t) might stray in to negative territory. The perturbative result Eq. (17a) is in excellent agreement with the numerics. Phase diagram of branching process coupled to an Ornstein-Uhlenbeck process for s = 1 and different values of D and β. We estimated numerically the critical R0 for some values of λ (symbols), see Fig. 4 . Critical values separate subcritical (below) and supercritical regimes (above). Solid lines (labelled as "theory") indicate the critical R0 for a given λ as approximated by Eq. (17a). These theoretical curves are first-order approximations in orders of Dλ 2 /β 3 , which explains the deviation of the numerical estimates from the theory. We expect larger deviations between true values and our approximation (17a) for larger ratios of D/β 3 . The mean-field approximation in (15) (dashed line) is common to the three sets of values since D/β is the same in all three cases.

As a second example of how noise can shift the critical point in unexpected ways, we consider a branching process in which the extinction rate switches spontaneously between two values. We call this random switching telegraphic noise [36] and write,

where the binary random variable T switches between the two values on > 0 and off = 0. An example trajectory is shown in Fig. 2 . Analogously to Eq. (13), we can immediately deduce the time-dependent expected number of secondary cases

with R 0 as defined by the distribution p k , Eq. (2). The waiting times between switching events are exponentially distributed with rates µ off to go from on to off , and µ on to go from state off to on ,

The switching rates µ on and µ off induce temporal correlations in the noise T (t) in the same way as the active fluctuations in the motion of run-and-tumble particles [37, 38] , Eq. (27). The master equation describing the telegraphic noise only is

where P (T ; t) is the probability distribution of T at time t. From Eq. (22) we derive the expected value of T at time t > 0:

In this set-up, the expected number of secondary cases E[K], Eq. (20) , switches between two values. If one of them predicts supercritical behaviour and the other one predicts subcritical behaviour, we cannot immediately deduce the criticality of the population that randomly switches between them. One way of determining the critical point is to demand that the rate with which offspring are produced in branching events equals that with which they go extinct. In each branching into k particles, k − 1 offspring are produced, so that the production rate of particles is

using Eq. (2). This production rate is unaffected by the state of the system. The extinction, on the other hand, depends on the state. As the transitioning times are exponentially distributed, the system spends on average 1/µ on amount of time in the T = off = 0 state and 1/µ off amount of time in the T = on state. This implies that at an arbitrary point in time, the system is in state T = off = 0 with probability µ off /(µ on + µ off ) and in state T = on with probability µ on /(µ on + µ off ). The effective extinction rate is therefore

(25) Equating this with Eq. (24) produces the criterion for the critical point,

However, Eq. (26) does not correctly predict the critical point as generally a larger population is affected by small extinction rates than by large extinction rates, because µ on and µ off are finite, so that the system lingers in either state. The mean field theory is expected to describe only the case of µ on , µ off s correctly, when the telegraphic noise changes so quickly that population size and state of the noise become uncorrelated. In its steady state, the telegraphic noise has a Pearson correlation coefficient ρ of

which indicates that correlations become irrelevant as µ on and µ off become large compared to s, the other event rate of the system. The derivation of Eq. (27) is presented in Appendix C 1. As in Sec. III, the expectation of the number of offspring averaged over all branching events is not a simple average of Eq. (20) , as it lacks a weighting by population size.

In order to capture whether outbreaks are supercritical or not, we inspect again the expected number of infected individuals over time E [N (t)]. For the calculation of E [N (t)], we use a Doi-Peliti field theory, derived in Appendix C, with initial condition N (0) = 1 and T (0) = on . We determine in which parameter region E [N (t)] grows over time (supercritical phase) and which it decreases (subcritical phase). The boundary between the two regions defines the critical hypersurface

which is shown in Fig. 7 . The direct comparison between the mean-field approach, Eq. (26) , and the exact result, Eq. (28), in Fig. 7 shows that the mean-field approximation predicts subcritical behaviour in regions where the dynamics are actually supercritical, i.e. in the regions between solid and dashed lines. As expected, Eqs. (26) and (28) coincide when µ on , µ off s, in which case they reduce to on /s = (R 0 − 1)(1 + µ off /µ on ).

We verified the shifted critical line using Monte-Carlo simulations, shown in Fig. 6 . 

While in the previous sections the branching process was coupled to different noises via a dynamic change in the extinction rate, here we study a type of noise introduced by the interaction between different populations. In the context of infectious diseases, we can think of population subgroups with different susceptibility, perhaps as a matter of lifestyle, behaviour or underlying health condition. As before, we leave various interpretations and applications to the reader and focus on analyzing the dynamics of an example process.

We consider a branching process that is coupled to another branching process. Individuals from two populations A and B with branching probabilities p kA and p kB , Eq. (1), respectively, change from one population to the other with transmutation rates µ A (from A to B) and µ B (from B to A),

An example trajectory is shown in Fig. 3 . The master equation of this process involves the joint probability P (N A , N B ; t), where N A and N B are the number of individuals of populations A and B respectively at time t. The master equation is made of three blocks describing each subprocess: two independent branching processes for sub-populations A and B, modelled in (3) and a coupling term that models the interaction (29) between the two populations,

Here, the term ∂ t P A (N A , N B ; t) is given by (3) replacing N by N A , p k by p kA and P (N ; t) by P (N A , N B ; t); the term ∂ t P B (N A , N B ; t) is given by (3) replacing N by N B , p k by p kB and P (N ; t) by P (N A , N B ; t); and the term ∂ t P AB (N A , N B ; t) captures the transmutation of indviduals in (29) ,

This coupling term describes how an individual of population A joins population B with rate µ A and how individuals from B convert to A with rate µ B . Fig. 9 for µA = 1 and µB = 2). The basic reproduction number R0A is above 1 and increases from left to right panels. The basic reproduction number R0B is adjusted in each panel to illustrate supercritical processes (top row), critical processes (middle row), and subcritical processes (bottom row). The supercritical cases shown in the top row are incorrectly predicted to be subcritical by the mean-field theory. Numerical estimates (symbols), based on 2 · 10 5 trajectories, are in good agreement with exact predictions (lines) in Eq. (D6). Fig. 8 for µA = µB = 1). The basic reproduction number R0A is above 1 and increases from left to right panels. The basic reproduction number R0B is adjusted in each panel to illustrate supercritical processes (top row), critical processes (middle row), and subcritical processes (bottom row). The supercritical cases shown in the top row are incorrectly predicted to be subcritical by the mean-field theory. Numerical estimates (symbols), based on 2 · 10 5 trajectories, are in good agreement with exact predictions (lines) in Eq. (D6).

To derive the dynamics of one of the two populations, say A, we marginalise the joint probability P (N A , N B ; t) by summing over M , which gives the probability that population A has N A individuals at time t,

By summing over N B in (30) we obtain

which shows that, from the perspective of sub-population A, its dynamics can be cast into a branching process (3) with slightly adjusted extinction rate and an additional influx, akin to spontaneous creation. The first two terms on the right hand side of Eq. (33) are indeed identical to the branching process in Eq. (3), the term parameterised by µ A corresponds to a spontaneous extinction and the last term, parameterised by µ B , is reminiscent of a spontaneous creation. However, the rates of the gain and loss terms of this creation differ and depend on E [N B (t)|N A (t)], which is a deterministic function of the stochastically varying size N A of sub-population A.

It is the only term that links the dynamics of the two sub-populations. In particular, it encapsulates the conservation of individuals by transmutation. The branching dynamics of sub-population B disappears from the dynamics of sub-population A otherwise.

If the branching processes of both populations are supercritical, we expect the coupled populations to remain supercritical, irrespective of the transmutation, as it conserves the total population size and cannot reduce it. Similarly if both processes are subcritical, as the transmutation cannot increase the population size either. However, the overall dynamics is not straightforward if the two populations lie in different criticality regimes. Without loss of generality, we assume in the following that basic reproduction numbers are R 0A > 1 and R 0B < 1, both defined by Eq. (2) with p k replaced p kA and p kB respectively. Is the joint population of A and B super-or subcritical?

A naive approach would be to consider how much time an individual spends as part of population A before joining population B and vice versa. As the waiting time between transmutations is exponentially distributed with parameters µ A and µ B respectively, an individual spends on average 1/µ A time in population A and 1/µ B time in population B. Thus, a weighted average of the two R 0 values is given by

and demanding that this weighted average is unity at the critical point determines the critical hypersurface as

which is shown in Fig. 10 as dashed lines. This ad hoc approximation ignores some important details of the interaction between the two sub-populations: Firstly, many more individuals are initiated in sub-population A, a bias not accounted for by the time-averaging taken above. Secondly, whenever a individual resides in A, many more branching events, namely those of its many offspring, will be characterised by R 0A . As in Sec. IV, only in the limit of large µ A , µ B with constant µ A /µ B , can we expect Eq. (35) to be correct. In order to find the critical point where the average total population size starts displaying exponential growth, we use a Doi-Peliti field theory, Appendix D. Given transmutation rates µ A and µ B , the border between the subcritical and the supercritical regime is the critical line defined by

shown in Fig. 10 

Branching processes are used to model a variety of avalanche-like dynamics ranging from neuronal activity to infectious diseases. One of the key points of interest is the prediction of the criticality of the dynamics, i.e. whether activity might diverge and continue forever or whether it will die out eventually. Here, we show in several examples that global noise in the branching dynamics induce correlations in the entire population which influence the growth or decline of the population over time.

The three examples are i) Sec. III, a branching process in which the subprocess of extinction follows an Ornstein-Uhlenbeck process, ii) Sec. IV, a branching process in which the subprocess of extinction displays telegraphic noise, i.e. its extinction rate randomly switches between two values and iii) Sec. V, two coupled branching processes where individuals can convert from one branching process and its parameters to another branching process with a different set of parameters.

In each case, we saw that mean-field arguments fail, even when they seem to capture much of the dynamics. In each case, correlations and fluctuations are important to be captured correctly, something that coarsegrained descriptions, such as SIR models, routinely ignore. Instead, we made use of more sophisticated, fieldtheoretical techniques to characterise the branching processes either perturbatively, as in the case of Ornstein-Uhlenbeck noise (Secs. III and B) or exactly (Secs. IV, V, C and D). This approach allowed us to determine effective critical values of R 0 that, surprisingly, are less than unity.

In all three examples, the effect of the noise is that the critical point dividing divergence from decline is shifted in unexpected ways. Mean-field approaches predict subcritical behaviour in some parameter regions where in fact supercritical dynamics occur. The reason for the failure of the mean-field theory is that it is based on averages that do not take the population size into account. When the extinction rate fluctuates, typically large populations are exposed to low extinction, and typically small populations are exposed to high extinction rates, unless the change of the extinction is fast compared to the process. And yet, at the critical point, as gains and losses are balanced, we expect the typical population sizes to be identical in the low and the high extinction rate regimes. In other words, according to this argument, the meanfield theory may fail to characterise the average number of offspring produced, but it should still predict the critical point correctly. However, critical or not, small, symmetric fluctuations in the time spent in a state of either high or low extinction rate have a disproportionate effect on the population size: Spending additional time ∆t in a state of low extinction creates many more offspring, given the initially large population size, than go extinct by spending the same additional time ∆t in the state of high extinction, or that are failed to be created when the time spent in the low extinction rate is reduced by ∆t.

A similar bias is visible in the expectation of the exponential E[exp(z t )] = 1 + E[z 2 t ] + . . . > 1 of a symmetric random walk z t with E[z t ] = 0. Somewhat ad hoc, the population size of the branching process with Ornstein-Uhlenbeck noise, Sec. III and Sec. B, may in fact be approximated by assuming that it is the exponential of the time integral of the instantaneous effective mass r [8] , and correspondingly averaged,

as the total, effective, instantaneous mass is r + λy t with r = s(1 − R 0 ), Eq. (A1), and E[y(t)] = 0. Using E[y(t)y(t )] = (D/β) exp(−β|t − t |) of the Ornstein-Uhlenbeck process [35] , this integral produces indeed the correct first order correction, Eq. (B15), (38) and, apparently, all higher order corrections of the population size, Eq. (B14). Given the ad hoc nature of this expression, we cannot confirm its validity to all orders, nor do we know whether correlation functions are correctly captured [cf. Eqs. (6) and (32) in 11].

In this paper we consider a branching process with an external noise as a basic model of infectious disease spreading. This model of unmitigated epidemic assumes an infinite population of susceptible individuals and hence it does not account for elements such as immunity or saturation. However, our results show an important aspect in disease spreading that, to our knowledge, has not been accounted for before [39] : fluctuations present in the transmission of the infection can shift the critical value of R 0 below 1. Therefore, we may observe an epidemic outbreak despite R 0 being smaller than 1. This needs to be taken into account when designing interventions aimed at containing an epidemic outbreak.

Moreover, we expose the failure of mean-field theories, which are widely used in epidemic modelling, to predict the critical point, see discussion above. The main reason for this failure is that mean-field theories and deterministic models, such as the classic SIR model, are not designed for capturing fluctuations and noise. These may provide useful approximations in the mist of an ongoing epidemic crisis [40] [41] [42] [43] [44] but lack the capacity to account for randomness. Therefore, incorporating the stochastic nature of disease spreading in epidemiologic models will provide a better understanding and prediction of the evolution of an outbreak [31] . We leave for future work the study of the effect of immunisation and the study of interventions that are able to contain the outbreak as well as the calculations of observables such as the peak prevalence (maximum value of infected individuals over time) and the final size of the outbreak (proportion of population infected at any time during the epidemic).

In practice, the failure of the mean-field theory means that a faithfully defined critical R 0 value which is obtained by ignoring fluctuations and correlations is unreliable. Even basic, handwaving arguments fail -we had initially expected the noise to push the critical value of R 0 above unity, because additional fluctuations might terminate a branching process by wiping out the last individual of a small, highly volatile population of an otherwise near-critical branching process. Instead, as long as the noise has a finite correlation time, so that the system lingers in a state of higher or lower extinction rate, there is an intrinsic bias towards larger populations. The most striking consequence of this bias is the critical basic reproduction number R 0 , in our naive definition Eq. (2), becoming less than unity. There is no unique definition of R 0 , yet in the case of unbiased Ornstein-Uhlenbeck noise, Sec. III, we cannot think of a redefinition of R 0 that renders its critical value unity.

Future research will focus on finding a more suitable observable or set of observables which constitute a sufficient predictor for the criticality of the noisy branching dynamics.

where the fields φ and φ are functions of time t and the binomial coefficient is zero if j > k. We assume that the system is initialised with one individual at time t 0 = 0, N (t 0 ) = 1, throughout. To calculate an observable O in the dynamics of the branching process, we perform the path integral

which generally involves the bilinear part in (A1a) and the nonlinear couplings q j . However, the observables that we are concerned with in this paper, do not involve the couplings q j , so we do not consider them beyond this point. The Gaussian model of the branching process follows from the bilinear part of (A1a), A 0 [φ, φ] = − dt φ (∂ t + r) φ, which gives the bare propagator from the Gaussian model [8] ,

where δ(ω + ω ) = 2πδ(ω + ω ) and

In (A3), the frequencies ω and ω are reciprocal to times t and t 0 under the Fourier transform convention,

and similarly for φ, with dω = dω/(2π). Our general approach in the three examples illustrated in Appendices B, C and D, where the extinction rate sp 0 is modulated by an external noise, is the following. We first derive the action that governs the dynamics of the external noise A Y with Y either an Ornstein-Uhlenbeck process (Y=OU), a telegraphic noise (Y=T) or a branching process (Y=BP). Then, we derive the action A int that describes the interaction between the noise Y and the branching process. Merging the three parts, the action that encapsulates all concurring processes is

In each of the actions of the three subprocesses there are, or may be, bilinear terms. We include those bilinear terms in the Gaussian model A 0 and group the rest in the perturbation A 1 , so that the overall action A is written as

To calculate an observable, such as the expected number of infected individuals E [N (t)] = φ(t) φ(t 0 ) , we then perform a perturbative expansion about the Gaussian model,

where the fields ψ, ψ represent the external noise.

Appendix B: Ornstein-Uhlenbeck Process

While some noisy processes can be described by a Doi-Peliti field theory (Sec. IV and V), where fields capture the time-dependent density of a degree of freedom, others are easier described using a Langevin equation and the response field formalism [45, 46] , where fields represent the degree of freedom itself. The Langevin equation of the Ornstein-Uhlenbeck process is

where y(t) is the degree of freedom of the Ornstein-Uhlenbeck process, β is the inverse persistence time, and η(t) is a Gaussian white noise with mean E [η] = 0, and correlator E [η(t)η(t )] = 2Dδ(t − t ).

Using the Janssen-DeDominicis-Martin-Siggia-Rose (JDMSR) response field formalism [45] [46] [47] [48] , the Langevin equation of the Ornstein-Uhlenbeck process, Eq. (B1) can be transformed into the action [46] 

which defines the field theory for the OU process. The field ψ represents the values of the random variable in the OU process. The auxiliary field ψ is not related to the ψ field. It is not to be confused with a Doi-shifted creator field. It is introduced in the JDMSR response field formalism purely to enforce the system to obey the OU Langevin equation, (B1), [46] . The bare propagator of the external noise iŝ 

which contains the Hilbert transform of a Gaussian. The Hilbert transform of a function f is defined as [49] H 

is the Hilbert transform of a Gaussian,

we obtain the result in (14) .

We couple the Ornstein-Uhlenbeck process in (B2) to the extinction rate of the branching process in (A1a), as described by the master equation, Eq. (11) . From (11), we derive the interaction term

which produces the nonlinear coupling

The overall action is then

which combines the Doi-Peliti and JDMSR formalisms, drawing on Eqs. (A1), (B2) and (B9). The subprocess governed by the action A BP + A int describes a branching processes with a factor ψ that modulates the extinction process. The expectation of any observable O of this subprocess is represented by a n-point correlation function O(ψ) BP that is a function of ψ,

The OU process renders ψ a random variable with a probability distribution given by (10) . In order to take this distribution into account, the expectation O(ψ) BP needs to be considered as an observable for the OU process,

In order to determine whether the dynamics are superor subcritical, we calculate the expected number of infected individuals N (t) = φ(t) φ(0) and look for its exponential growth and decay respectively. The extended action (B11) introduces the following loop corrections to this propagator,

The one-loop correction in (B14a) iŝ

so that the propagator is

The pre-factor of e −rt indicates that, in this approximation, the critical point remains at 0 = r = s(1 − E[K]).

Calculating only this one-loop approximation is indeed insufficient to see a shift in the critical point. Calculating, on the other hand, all loop corrections in (B14) is doable in principle, but calculating a closed form expression for each of them and doing appropriate bookkeeping is an arduous task, in particular for entangled loops such as the ones in (B14c). Instead, we find that the Dyson sum, which includes all loop corrections of the form (B14a)-(B14b), and excludes entangled loops, known as the non-approximation, provides a better approximation to the expected number of infected individuals E [N (t)]. This is,

where from Eq. (B17b) to Eq. (B17c), we calculated a geometric sum of the one-loop correction from Eq. (B15). The exponential growth rate in Eq. (B17d) gives the critical point approximated by the Dyson sum,

which forms a hypersurface in r-λ-D-β space and which is transformed into the critical hypersurface in Eq. (17a).

"on", and χ, χ for n particles in state "off" with action A T ,

The values on and off , Eq. (19) , can be chosen arbitrarily. For simplicity, we set off = 0.

To couple the telegraphic noise T (t) to the branching process such that its value is added to the death rate, we use the master equation from the previous model (branching with OU noise), Eq. (11), where we replace λy by T (t). The value of the telegraphic noise T (t) is represented by on ψ † ψ in the field theory. The creator field is then Doi-shifted ψ † = ψ + 1, which leads to the following term in the action

As a result, the total extinction rate switches between values sp 0 and sp 0 + on in intervals that are exponentially distributed with rates µ on and µ off , respectively. The overall action of the process is then

(C4)

In this field theory we need the bare propagator of the driving noise, (C6)

In order to identify the critical point, we calculate the expected number of individuals,

While the first term in (C7b) is simply

the second term in (C7b) is more complicated and can be represented in Feynman diagrams as a Dyson sum, φ(t) φ(0) ψ(0) = + + + · · · (C9a) = − on e −iωt (−iω 1 + µ on ) (−iω + r)(−i(ω − ω 1 ) + r)((−iω 1 + µ off )(−iω 1 + µ on ) − µ off µ on )

− on e −iωt (−iω + r + µ on ) (−iω + r) (−iω + r + µ off )(−iω + r + µ on ) − µ off µ on + on (−iω + r + µ on ) dω ,

where from Eq. (C9b) to Eq. (C9c), a geometric sum over loop corrections is calculated. Using the abbreviations, 

In fact, to determine the critical point, we do not need to calculate E[N (t)] explicitly because the boundary between supercritical and subcritical regimes is marked by a change of sign of the imaginary part of the complex poles of φ(t) φ(0) ψ(0) in Fourier space, i.e. the poles of the ω integral in Eq. (C9c). The equation for the critical hypersurface is then given by the equation that set the imaginary part of the ω-poles equal to zero: on + r(r + µ off + µ on ) r + µ on = 0,

which transforms into Eq. (28), using Eq. (A1b).

The Pearson correlation coefficient ρ XY of two random variables X and Y is defined as [53] 

In the case of the telegraphic noise T , we are interested in X = T (t) and Y = T (t ). If we assume that t > t > 0, then E[T (t)T (t )|T (0) = on ] can be written in the field theory as E[T (t)T (t )|T (0) = on ] = 2 on ψ(t)ψ † (t )ψ(t )ψ † (0) .

Once the Doi-shift ψ † = ψ + 1 is performed, the only remaining non-vanishing term can be calculated as

= 1 (µ off + µ on ) 2 µ on + µ off e −(µ off +µon)t µ on + µ off e −(µ off +µon)(t−t ) .

Changing the order of t and t such that t > t amounts to swapping each t for a t and vice versa in Eq. (C15c).

In the steady state, only the difference between t and t enters and a single observation of T will be a Bernoulli experiment, in which T = on is drawn with probability µ on /(µ on + µ off ). Hence its expectation equals E[T ] = on µ on /(µ on +µ off ) and its variance equals V[T ] = 2 on µ on µ off /(µ on + µ off ) 2 , consistent with the second moment Eq. (C15c) for t − t → ∞ and t → ∞. Combining Eq. (C15c) with the exepectation and variance of T , we find the Pearson correlation coefficient in Eq. (27) .

The two populations of branching species A and B are represented by fields φ, φ and ψ, ψ respectively, and both follow the dynamics in the branching action (A1a) with coefficients r A , q jA and r B , q jB , respectively. The action of the coupled branching process then includes the branching actions A BP [φ, φ] and A BP [ψ, ψ], plus a third term that describes the interactions in (29) between the two species,

This introduces additional mass terms, so that the bare propagators read=

as well as the interaction vertices µ B and µ A (D3)

The overall action of the coupled branching processes is then

Since the interaction terms in the action are all bilinear, we can include them in the Gaussian model A 0 . Then, to find the critical point, all we need are the propagators,

The first moments of the particle numbers are then derived using inverse Fourier transforms:

The propagators, Eq. (D5), readily encode the critical point of the coupled branching process. Both contour integrals have two poles,

which are purely imaginary. In the subcritical regime, their imaginary part is negative, while in the supercritcal regime at least one of the poles has a positive imaginary part. Thus the critical hypersurface is determined as

which is transformed into (36) .

The Theory of Branching Processes

Branching processes

Evolutionary Dynamics

A Short History of Mathematical Population Dynamics

Neutron Fluctuations: A Treatise on the Physics of Branching Processes

Random Processes in Nuclear Reactors

From neuronal spikes to avalanches -effects and circumvention of time binning

Branching Processes in Biology

Incorporating human movement data to improve epidemiological estimates for 2019-ncov

Some challenges in monitoring epidemics

Stochastic Processes in Physics and Chemistry

Noise-induced transitions

Theory Exp. 2021

Critical dynamics

Grundzüge einer allgemeinen Theorie der linearen Integralgleichungen

Analyse mathématique sur les probabilitś des erreurs de situation d'un point (acad. r. sci

We thank Andy Thomas for his brilliant technical support, Nanxin Wei and Guillaume Salbreux for fruitful discussions. J.P. was supported by an EPSRC Doctoral Prize Fellowship. R.G.M.'s work was funded by the European Research Council under the EU's Horizon 2020 Programme, grant number 740269.

All authors were involved in all aspects of the conceptualization and investigation as well as the writing and editing of the manuscript. JP acquired partial funding for this project through an EPSRC Doctoral Prize Fellowship. Data visualization and other graphics were created by RGM and JP. RGM and GP acquired computing resources through Imperial College London.Appendix A: Branching process field theory We recall from [8] the Doi-Peliti field theory for the continuous-time branching process with arbitrary timeindependent offspring distribution, where K ∈ {0, 1, . . . } is the number of offspring produced at a branching event with probability P (K = k) = p k , Eq. (1). The waiting times between two branching events is exponentially distributed with rate s. The action of this field theory isAppendix C: Telegraphic noiseTo cast the telegraphic noise in a field theory, we consider two species, "on" and "off". The state of the telegraphic noise is given by the number of particles m and n of species "on" and "off" respectively. Particles transmute between the two species according to (21) such that the total number of particles is conserved. In our case, the total number of particles is m + n = 1 since, initially, the telegraphic noise is T (t 0 ) = on .Using a bra-ket notation, m, n| and |m, n , and ladder operators a, a † , b, and b † with commutators [a,and a|m, n = m|m − 1, n , a † |m, n = |m + 1, n , b|m, n = n|m, n − 1 , and b † |m, n = |m, n + 1 , the master equation (22) can be turned into an equation for the probability generating function |M(t) := {m,n} P (m, n; t)|m, n [51] ,Building on work by Peliti [52] , Eq. (C1) can be turned into a field theory for fields ψ, ψ for m particles in state