Self-driven criticality in a stochastic epidemic model
SSelf-driven criticality in a stochastic epidemic model
Gil Ariel and Yoram Louzoun
Department of Mathematics, Bar-Ilan University, Ramat Gan 52000, Israel
We present a generic epidemic model with stochastic parameters, in which the dynamics self-organize to a critical state with suppressed exponential growth. More precisely, the dynamics evolveinto a quasi-steady-state, where the effective reproduction rate fluctuates close to the critical valueone, as observed for different epidemics. The main assumptions underlying the model are that therate at which each individual becomes infected changes stochastically in time with a heavy-tailedsteady state. The critical regime is characterized by an extremely long duration of the epidemic.Its stability is analyzed both numerically and analytically.
Compartmental epidemics models, and their many variations or derivatives have been proven useful in understand-ing, analyzing and predicting real epidemic outbreaks [1, 2]. However, fitting such models to the observed patient ordead counts has proven challenging [3]. One of the main reasons for that is that compartmental models are essentiallyexponential, at least locally in time [1, 2, 4–6], while observed data are often not. Exponential dynamics emergesince at any given time in the evolution of the epidemics, the equation for the dynamics can be linearized around itscurrent state, suggesting an exponential growth or decay of variables (except at particular time points such as the localmaximum of infected individuals). The predicted exponential growth/decay motivates the commonly used notion ofthe effective reproduction rate, R t , or equivalently, doubling time [1, 2]. The effective reproduction rate describesapproximately the instantaneous exponential rate of change in the number of infected, hospitalized, deceased or othertypes of individuals. Fitting R t to real data is not straight-forward, and several methods have been proposed andapplied [7–9].Often, the epidemic dynamics seem to be at, or close to the critical state R t = 1 [6, 10–16]. Consequently, thenumber of new cases per day (after some smoothing to eliminate weekly periodicities) is constant or linear. Indeed,several authors studied the dynamics of epidemics, both in compartmental and network models, assuming that theepidemic is poised at the critical threshold between exponential growth and decay [4, 5, 14, 17]. This dynamicalpattern can be explained by assuming that governments modify social distancing rules or, alternatively, people adapttheir behavior according to the perceived spread of the epidemics to fine-tune R t [11, 16]. However, such negativefeedback typically takes effect on long time scales, possibly up to years [16]. Alternatively, Stollenwerk and Jansen[12, 13] suggested a sand-pile-type model that exhibits self-organized criticality. The model assumes that the epidemicspreads on a square 2D lattice. Criticality is due to a vanishingly small rate of mutation to a deadly strain. Recentversions assume that the epidemic spreads to neighbors only once the viral load is above a threshold (thus theanalogy to a sand-pile) [10] or that lattice sites are partially isolated communities (cliques) [15]. Although thesemodels are applicable for specific examples (e.g. meningococcal disease [12, 13]), they are not as generic as standardcompartmental models. Moreover, there is no evidence supporting their main assumptions in general.We suggest a generic mechanism driving epidemics towards criticality. The main new assumption is that theinfectivity of each individual (i.e. the probability of each individual to get infected upon the meeting with an infectedindividual) is a time-varying correlated random variable. Three dynamical regimes or phases are identified, dependingon parameters: (i) Exponential growth, (ii) Exponential decay, and (iii) a quasi-steady critical state - A novel regimein which the dynamics naturally evolves to a steady-state that is close to critical, i.e., R t fluctuates close to 1. In thisregime, the exponential dynamics are suppressed, resulting in a time-series data that is polynomial, typically appearingto be constant or linear over long time scales. The critical state does not require fine-tuning of parameters to haveparticular prescribed values. The only condition is for a large enough tail in the infectivity probability distribution.For simplicity, we concentrate on the SIR modeling point of view, assuming that the evolution of the disease ineach individual follows the three states (S)usceptible -¿ (I)nfected -¿ (R)ecovered [1, 2]. In addition, we assume awell-mixed population [1]. These assumptions underlie considerable simplifications, and realistic epidemic modelingrequires more detailed models [1, 18, 20]. However, our main goal is to show that the new self-organized criticalregime is generic and does not depend on specific modeling details. Hence, we restrict our focus to the simplest, mostuniversal approach. The mechanism proposed here applies similarly to more complex models. Similarly, for the sakeof simplicity, we assume that the variability is only on the infected side, and not on the infecting side and a finitecorrelation time.The SIR model dynamics are determined by three parameters: The population size N , the infectivity rate λ , andthe recovery rate, taken with no loss of generality to be 1. Our only departure from the standard SIR model is inrelaxing the assumption that the infectivity rate of an individual is constant. Instead, we assume it follows a stochasticprocess, i.e., it is changing randomly in time. In other words, the infectivity rate of individual k is replaced by astochastic process λ k ( t ). Thus, each infected individual can infect each susceptible individual with a rate proportional a r X i v : . [ n li n . AO ] J a n to the susceptible infectivity, λ k ( t ).Different biological and social mechanisms can be proposed to explain a random evolution of the individual suscep-tible infectivity rate. For example, individual differences in health, behavior, social distancing measures, and more.Indeed, previous models considered the effects of quenched heterogeneity (i.e., assume that individuals have differentproperties, but they are fixed in time) [19–22]. Several authors considered deterministic time dependence of λ k ( t ),for example due to seasonality or decline in vaccination rates [16, 23, 24]. Stochasticity has also been considered,either through a randomly evolving network [25], by additive and multiplicative noise [23, 26, 28, 30] or by stochastic(typically normal) fluctuation of parameters [27, 29]. The main conclusion of these approaches is that assuming pop-ulation growth (possibly with a random but strictly positive growth rate), stochastic SIR models admit a steady-statesolution. This solution does not exist if the population size N is fixed. Intuitively, the idea underlying these modelsis clear: at the steady-state, the average population growth compensates for the average rate of new infections. Thissteady-state is reached on the time scale of the population growth, typically considerably longer compared to theepidemic time scale.Here, we assume that infectivity rates λ k ( t ) are independent, identical stationary stochastic processes with twoimportant properties: i. univariate marginal distributions that have a power-law tail, ii. a finite correlation time.Note that, while individual infectivity rates vary, the statistics of all individuals are identical and constant in time. Inparticular, the process is autonomous, i.e., there is no explicit time dependence. To be specific, the model is definedas follows.The probability that a susceptible individual k becomes infected in the period [ t, t + dt ) is p k (S → I) = λ k ( t ) i ( t ) dt + o ( dt ). The probability that an infected individual k recovers in the period [ t, t + dt ) is p k (I → R) = dt + o ( dt ).There are different ways to implement a stochastic process satisfying the two conditions described above. We choosea particular realization which is convenient to simulate and admits some analytical analysis. The main idea is todefine the characteristic infectivity times µ k ( t ) = 1 /λ k ( t ) instead of the rates. We assume that µ k ( t ) are independentstationary processes with gamma univariate marginal distribution, µ k ( t ) ∼ Γ( α, β ) (1)and an exponential autocorrelation function C ( t, s ) = (cid:104) µ k ( t ) µ k ( s ) (cid:105) − (cid:104) µ k ( t ) (cid:105) (cid:104) µ k ( s ) (cid:105) = e −| t − s | /τ . (2)It is generally accepted that realistic infectious periods are not exponential but are better approximated by a gammadistribution [31]. Note that an ensemble of exponential random variables in which the mean has a gamma distributiondoes not necessarily obey a gamma distribution, so our assumption is not the same as in [31]. This difference isnot essential, as the gamma model is only used to approximate a more complex reality. Recall that for a gammadistribution with shape parameter α and rate β , the average is (cid:104) µ k ( t ) (cid:105) = α/β . In addition, the instantaneousdistribution of infectivity rates λ k ( t ) = 1 /µ k ( t ) has an inverse gamma distribution, characterized by a tail that decayswith a power α + 1.Writing the SIR model as a continuous-time Markov process, we simulate the dynamics using the Guilesspie al-gorithm [34]. Synchronous simulations give similar results. The fraction of Susceptible, Infected and Recoveredpopulations are denoted s ( t ), i ( t ) and r ( t ), respectively, with s ( t ) + i ( t ) + r ( t ) = 1. In simulations, we fix the shapeparameter 1 < α <
2, which implies that infectivity rates have a finite mean, but infinite variance. The rate parame-ter β is chosen to have a given average transition time (cid:104) µ k ( t ) (cid:105) . Initial conditions are 10 infected individuals and norecovered. The initial distribution of infection times is taken to be the steady-state distribution Γ( α, β ). Thus, theprocess for µ k ( t ) is stationary. See the appendix for further implementation details.Fig. 1a shows the time evolution of s , i and r for a population of 10 individuals. Parameters are α = 1 . (cid:104) µ k ( t ) (cid:105) = 1 . β (cid:39) .
81 and (cid:104) λ k ( t ) (cid:105) (cid:39) .
7) and τ = 4. Following a short exponential transient that follows theSIR model with either equal or quenched infectivities dynamics, most individuals that started with a high infectivityrate are infected, and the rate of new infections decreases (fig. 1b). As a result, the fraction of infected individualsdecreases, but on a considerably longer scale compared to previous models with non-critical dynamics.Fig. 1b depicts our main result, showing that the effective reproduction rate R t stabilizes to a value that fluctuatesaround the critical value R t = 1. An intuitive explanation is a negative feedback mechanism: Initially, the averageinfectivity rate is larger than 1. As s is close to 1, the reproduction rate R t > λ k (0) are more likely to be infected duringthis exponential growth phase. As a result, given the heavy tail of the distribution of λ k , the average infectivity ratein the remaining susceptible population decreases, decreasing R t . On the other hand, if the rate of new infections i (cid:48) ( t )is sufficiently large (compared to 1 /τ ) so that R t <
1, then the number of infected individuals decreases exponentially.Since there are not many new infections, and because the distribution of infectivity rates returns to the steady-stateat a characteristic time τ , the average infectivity rate increases back towards its unperturbed steady-state average,which is greater than one. Note that an essential requirement for such a mechanism is a heavy-tailed distributionof the infectivities, otherwise, the reduction of the right-hand tail would not affect the average infectivity for largeenough populations. In addition, if the correlation time τ is too small or large, then the infectivity rate will not returntowards its steady-state value.Overall, if the increase and decrease in the average infectivity rate occur on comparable time scales, then thedynamics self-organize to the critical value. In the SIR model, the fraction of susceptible individuals s ( t ) decreasesmonotonically with time. Denoting averages over susceptible individuals by (cid:104)·(cid:105) S , the reproduction rate is given by R t = (cid:104) λ k ( t ) (cid:105) S s ( t ). Thus, criticality is maintained through an increase in the average infectivity of each susceptible,through a stronger bias for infections in the tail of the distribution (dashed blue curve in Fig. 1b).The duration of the critical state can be orders of magnitude longer than in the annealed analog. However, it doesnot correspond to a steady-state due to the finite population. Accordingly, we term this regime a quasi-steady criticalstate. To understand the different regimes that emerge in this model, we estimate the model phase diagram. Fig. 2ashows the three dynamical regimes (exponential growth, exponential decay or quasi-steady criticality) as a functionof the average infectivity time and the shape parameter. We refer to the exponential decay case as an instance inwhich the overall recovered population (when the infected population reaches 0) is smaller than 1%. Otherwise, theepidemic spreads (with exponential or critical growth). The regime is labeled critical if R t is between 0.95 and 1.05for at least half of the duration of the epidemic (i.e., until i ( t ) = 0). Otherwise, the regime is labeled as exponentialgrowth. Note that these definitions contain multiple arbitrary parameters, which may affect the precise borders ofthe different regimes. However, the main properties of the phase diagram are not affected by the details of the formaldefinitions.We see that the critical regime inhabits a large section of the phase diagram. Not surprisingly, it is highly dependenton α . When alpha approaches 1, the dynamics are critical for a very wide range of beta values. For high alpha values,the model converges to the classical SIR model distinction of exponential increase for high infectivities (lower timefor infection), and exponential decrease for low infectivities (Fig. 3).Fig. 2b shows the dependence of the dynamics on the characteristic time-scale of the stationary gamma process, τ .All other parameters are the same as in fig. 1. The final value of r , indicating the overall fraction of the populationwhich was infected, decreases monotonically with τ . However, the duration of the epidemic depends non-monotonouslyon τ . As the critical state is characterized by long epidemic duration, the figure shows that the stability of the R t = 1state increases for τ ≥ τ = 20 −
50. The sharpness of the drop in the duration suggests aphase transition. Finally, fig. 4 shows the dependence on the population size N and on the initial number of infectedindividuals I (0). We note that the duration of the epidemics strongly depends on N and is approximately proportionalto ln N .To further gain insight into the stability of the critical regime, we derive self-consistent equations for the existenceof a steady-state. We take advantage of a particular realization of stationary gamma processes for the case in whichthe shape parameter α is equal to half an integer. From [38, 39], the solution to the following stochastic differentialequation, dµ t = − τ (cid:18) µ t − αβ (cid:19) dt + 2 √ βτ √ µ t dW t , (3)satisfies (1) and (2), where W t is the standard Brownian motion. We assume that the population is at a steady-state,in which the rate of individuals leaving state S (i.e., they become infected) is constant. In other words, µ t is a diffusionprocess that is killed at a rate that is proportional to µ − t . To uphold the steady-state, individuals are reintroducedto the susceptible population at a constant rate. Thus, one can write the Fokker-Plank equation associated with theeffective steady-state distribution p ( µ ). The analysis is detailed in the appendix. The steady-state distribution maybe written as an asymptotic series p ( µ ) = p ( µ ) + 12 τ i critical q ( µ ) + O ( i ) , (4)where p ( µ ) is the density of the Gamma distribution Γ( α, β ) and i critical is the fraction of infective individuals duringthe steady-state. Then, q ( µ ) satisfies an inhomogeneous confluent hypergeometric equation whose solutions can bewritten in terms of a generalized gamma function. With α = 3 /
2, we find that a steady state can be obtained for alarge range of parameters (See supplementary material). In contrast, for α = 2 and α = 5 /
2, a steady-state cannotbe obtained for long durations that are at least of order ln N .The work presented here is the first step in understanding how temporal fluctuations in individual propertiesaffect the dynamics of epidemics. Some fundamental questions are still open. For example, deriving a rigorouslimiting equation is essential for a thorough understanding of the self-organized principle we suggest. In addition,linking the self-organizing principle to similar known physical processes, in particular the well-studied self-organizedcriticality [10, 12, 13, 15] and the stable steady-state inferred in stochastic models (with population growth) [16, 23–26, 28, 30] is also of importance. Relating the model proposed here to realistic epidemic data requires introducingmore complicated models to take into account several details that have been proven important to realistic modeling[31], such as interaction networks, heterogeneity (e.g. by age) and coupling between infectivity rates and the epidemicdynamics. These issues are beyond the scope of the current manuscript.To summarize, we showed that a small change in the well-studied SIR model in which the infectivity rates (or rathertimes) of each individual are changing stochastically in time can qualitatively change the dynamics. For simplicity,we assumed a minimal model with a well-mixed population, statistically identical individuals, and a stationary,independent random process for the infection time of each individual. Our results demonstrate that the effect ofincluding temporal fluctuations in the infectivity times is drastic. In particular, the dynamics admit a new regimethat is not present in previous compartmental models, in which the dynamics evolve spontaneously into a quasi-steadycritical state in which the reproduction number fluctuates close to 1. Without a doubt, fitting and understandingrealistic epidemic data require further generalizations to more complicated models. Yet, this work takes the first stepsin understanding the implications of the self-organized quasi-stable critical regime introduced here, which can haveanalog consequences in other systems, such as population dynamics and kinetics of chemical reactions. Acknowledgments
We thank Rainer Klages and Jeremy Schiff for comments and discussions. [1] Hethcote HW. SIAM review , 599 (2000).[2] Murray JD. Mathematical biology: I. An introduction . Springer Science and Business Media (2007).[3] Akira O, Namekawa Y and Fukui T. Progr. Theo. Exp. Phys., ptaa148 (2020).[4] Radicchi F and Bianconi G. arXiv :2007.15034 (2020).[5] Ben-Naim E and Krapivsky PL. Phys. Rev. E , 050901 (2004).[6] Rhodes CJ, Jensen HJ and Anderson RM. Proc. Royal Soc. London. Series B: Biological Sciences , 1639 (1997).[7] Ma J. Infectious Disease Modelling , 129 (2020).[8] Manrique-Abril FG, Tellez-Pinerez C and Pacheco-Lopez M. F1000Research , 868 (2020).[9] Hotz T, Glock M, Heyder S, Semper S, Bohle A and Kramer A. arXiv :2004.08557 (2020).[10] Contoyiannis Y, Stavrinides SG, Hanias MP, Kampitakis M, Papadopoulos P, and Potirakis S. arXiv :2004.00682 (2020).[11] Gans JS. National Bureau of Economic Research Working Paper Series, w27632 (2020).[12] Stollenwerk N, and Jansen VA. Physics Letters A , 87 (2003).[13] Stollenwerk N. In Recent Advances in Applied Probability, 455. Springer, Boston, MA (2005).[14] Grassberger P.Mathematical Biosciences , 157 (1983).[15] Ion S and Marinoschi GA. Discrete Cont. Dyn. Sys. B , 383 (2017).[16] Brett T, Ajelli M, Liu QH, Krauland MG, Grefenstette JJ, van Panhuis WG, Vespignani A, Drake JM and Rohani P.PLoS comput. biol. , e1007679 (2020).[17] Horstmeyer L, Kuehn C and Thurner S. Phys. Rev. E , 042313 (2018).[18] Epstein JM. Modelling to contain pandemics. Nature , 687 (2009).[19] Lachiany M and Louzoun Y. Phys. Rev E , 022409 (2016).[20] Britton T, Ball F and Trapman P. Science , 846 (2020).[21] Muchnic L, Yom-Tov E, Levy N, Rubin A and Louzoun Y. arXiv :2008.00445 (2020).[22] Miller JC. Phys. Rev. E , 010101 (2007).[23] Kloeden PE and Potzsche C. Math. Meth. App. Sci. , 3495 (2015).[24] Kloeden PE and Kozyakin VS. MESA , 159 (2011).[25] Ochab JK and Gora PF. Eur. Phys. J. B, , 373 (2011).[26] Tornatore E, Buccellato SM and Vetro P. Physica A: Stat. Mech. App. , 111 (2005).[27] Dureau J and Kalogeropoulos K. Biostatistics , 541 (2013).[28] Caraballo T and Colucci R. Commun. Pure and Applied Analysis , 151 (2017).[29] Faranda et al. Chaos , 051107 (2020).[30] Jiang D, Yu J, Ji C and Shi N. Math. Comput. Model. , 221 (2011).[31] Wearing HJ, Rohani P and Keeling MJ. PLoS Medicine , e174 (2005).[32] Arnold L. Random dynamical systems . Springer, Berlin, Heidelberg (1995).[33] Klages R. Europhys. LEtters , 796 (2002).[34] Allen LJS. Infectious Disease Modelling , 128 (2017).[35] Chen-Charpentier BM and Stanescu D. Math. Comput. Model. , 1004 (2010).[36] Kloeden PE and Platen E. Numerical solution of stochastic differential equations . Vol 23. Springer Science and BusinessMedia (2013).[37] Wolpert RL and Ickstadt K. Biometrika , 251 (1998).[38] Wolpert RL. Research notes, 2011-10-05, Draft 69 (2011).[39] Wolpert RL. Lecture notes (2016).
FIG. 1. Example simulation results with a population of 1M individuals, average infectivity rate 2.7 and shape parameter α = 1 .
3. Time units are the inverse rate of recovery. (a) The time evolution of the fraction of susceptible, s ( t ), infected, i ( t ) andrecovered r ( t ) individuals on a semilog scale. (b) After a short relaxation period, the effective reproduction rate R t fluctuatesclose to the critical value 1. As R t remains constant while s decreases, the average infectivity rate has to increase. FIG. 2. phase diagram. (a) The regime of the epidemic (exponential growth, exponential decay, or self-organization to thecritical state) for different average infectivity times and shape factor α . The critical regime occupies a large fraction of thephase space and does not require fine-tuning of parameters. (b) Dependence on τ . The figure shows the final values of r ,indicating the overall fraction of the population which was infected (solid black line), and the duration of the epidemic (dashedred line) as a function of τ . As the critical state is characterized by long epidemic duration, we see that the stability of the R t = 1 state increases for τ ≥ τ = 20 −
50, suggesting a phase transition.
FIG. 3. Results for annealed dynamics, τ → ∞ . All other parameters are the same as in fig. 1. (a) s ( t ), i ( t ) and r ( t ) on asemilog scale. The duration of the epidemic is an order of magnitude shorter compared to the critical regime. Time units arethe inverse rate of recovery. (b) Both the reproduction number r ( t ) and the average infectivity rate decrease monotonically intime. (c) The phase diagram shows that criticality is only obtained at the exact boundary between the exponential growth andexponential decay regimes. FIG. 4. (a) Dependence on N . Both the overall fraction of the population infected and the epidemic duration increaselogarithmically with N . However, while the former varies only by about 20%, the duration changes by a factor of five. (b)Results do not depend on the initial number of infected individuals for all I (0) ≥ Appendix: Implementation details
Recall that for each k , µ k ( t ) is a stationary process with univariate marginal distribution µ k ( t ) ∼ Γ( α, β ) andexponential auto-correlation function exp( −| t − s | /τ ) for µ k ( t ) and µ k ( s ). We approximate µ k ( t ) by a piece-wiseconstant jump process with step ∆ t using the discrete AR(1) process presented in [37–39]. To be precise, let ρ =exp( − ∆ t/τ ) and take µ k (( i + 1)∆ t ) = ρµ k ( i ∆ t ) + ζ i , (5)where, G i ∼ Γ( α,
1) (6) N i ∼ Pois((1 − ρ ) G i /ρ ) (7) ζ i ∼ Γ( N i , β/ρ ) . (8)Here, Pois( L ) denotes a Poisson distribution with mean L .The SIR model is simulated in steps of ∆ t . During the i ’th step, we assume that S → I transitions have an exponentialdistribution with constant rates λ k = 1 /µ k ( i ∆ t ). I → R transition rates equal 1. Within each time period ∆ t , the SIRdynamics is modeled using the Guilesspie algorithm as follows. For simplicity, we consider the first period, i = 0 and t = 0. Denote by A S ( t ) and A I ( t ) the sets of susceptible and infected individuals at time t , respectively.1. Let Λ = 1 N | A I ( t ) | (cid:88) k ∈ A S ( t ) λ k , (9)where | A | denotes the number of elements in a set A .2. Let dt S2I ∼ Exp(1 / Λ), an exponentially distributed random variable with mean 1 / Λ. Let dt I2R ∼ Exp(1 / | A I | ). dt = min { dt S2I , dt
I2R } .3. If t + dt > ∆ t , then no transitions occur until the end of the period [0 , ∆ t ] and the process jumps to t = ∆ t .The current time-step ends.4. Otherwise, a transition occurs. If dt I2R < dt
S2I , then a recovery event occurs. Let I denote an element from A I ( t ) chosen at random with uniform distribution. Then, individual I recovers: A I ( t + dt ) ← A I ( t ) − { I } and t ← t + dt . Return to 1.5. Otherwise, an infection event occurs. Let u ∼ U (0 , , I = arg max I (cid:40) (cid:80) k ∈ A S ( t ) ,k ≤ I λ k (cid:80) k ∈ A S ( t ) λ k < u (cid:41) . (10)Then, individual I becomes infected: A S ( t + dt ) ← A S ( t ) − { I } , A I ( t + dt ) ← A I ( t ) ∪ { I } and t ← t + dt . Returnto 1.Note that the algorithm makes use of the no-memory property of exponential random variables. Thus, cutting off thedistribution for dt if t + dt > ∆ t and redrawing it does no change the distribution of dt (unless Λ or | A I ( t ) | change). Appendix: Self-consistent conditions
We take advantage of a particular realization of stationary gamma processes, i.e., stochastic processes with univariatemarginal distribution and an exponential auto-correlation function. Following [38-39], if the shape parameter α isequal to half an integer, then one can write such a stationary gamma processes as a stochastic differential equation.Note that this realization is different than the AR(1) process defined in Eq. (4), i.e., it may have different higher( ≥
3) moments. Specifically, for α = n/ β >
0, the solution of the stochastic differential equation dX t = − τ (cid:18) X t − n β (cid:19) dt + 2 √ βτ (cid:112) X t dW t , (11)0satisfies [38,39], X t ∼ Γ( n/ , β ) (12) (cid:104) X t X s (cid:105) − (cid:104) X t (cid:105) (cid:104) X s (cid:105) = e −| t − s | /τ (13)The Fokker-Plank equation associated with (11) is ∂p ( x, t ) ∂t = 2 τ ∂∂x (cid:20)(cid:18) x − n β (cid:19) p ( x, t ) (cid:21) + 2 βτ ∂ ∂x [ xp ( x, t )] . (14)Our goal is to model the distribution of infectivity times (here denoted X ) in the susceptible population. We assumethat the population is at a quasi steady state, in which the rate of individuals leaving the state S (i.e., they becomeinfected) is constant. In other words, X t is a diffusion process that is killed at a rate that is inversely proportional tothe infectivity time, 1 X t (cid:82) ∞ x − p ( x, t ) dx . (15)In order to uphold the steady state, or in other words, renormalize p ( x, t ) to have unit integral, individuals arereintroduced to the susceptible population with a constant rate. Overall, the Fokker-Plank equation associated withthe effective steady-state process is ∂p∂t = 2 τ ∂∂x (cid:20)(cid:18) x − n β (cid:19) p (cid:21) + 2 βτ ∂ ∂x [ xp ] − ρ (cid:34) x (cid:82) ∞ x − p ( x ) dx − (cid:35) p ( x ) , (16)where ρ ≥ ∂p/∂t = 0,2 τ (cid:20)(cid:18) x − n β (cid:19) p (cid:21) (cid:48) + 2 βτ [ xp ] (cid:48)(cid:48) − ρ (cid:34) x (cid:82) ∞ x − p ( x ) dx − (cid:35) p ( x ) = 0 . (17)Taking ρ = 0, we verify that the density of a gamma distribution with rate parameter β and α = n/ p ( x ) = 1Γ( α ) β α x α − e − βx , (18)is indeed a steady state solution. Here, Γ( · ) is the Gamma function. Multiplying (17) by τ / (cid:15) = τ ρ/ L p ( x ) − (cid:15)β (cid:34) x (cid:82) ∞ x − p ( x ) dx − (cid:35) p ( x ) = 0 , (19)where L f = (cid:104)(cid:16) βx − n (cid:17) f (cid:105) (cid:48) + [ xf ] (cid:48)(cid:48) , (20)is the forward Kolmagorov operator associate with (11). Assuming (cid:15) (cid:28)
1, we expand p ( x ) in an asymptotic series, p ( x ) = p ( x ) + (cid:15)q ( x ) + O ( (cid:15) ) . (21)Substituting into (19) and expanding in powers of (cid:15) , L q ( x ) − β (cid:34) x (cid:82) ∞ x − p dx − (cid:35) p + O ( (cid:15) ) = 0 , (22)where we used L p = 0. In the following, the high order O ( (cid:15) ) term is neglected.As α = n/ < α < α = 3 /
2. Theleading order term is, p ( x ) = 2 √ π β / x / e − βx . (23)1Substituting p into (22), we obtain a second order ordinary differential equation for the leading order perturbationterm q ( x ), βq + (cid:18) βx + 12 (cid:19) q (cid:48) + xq (cid:48)(cid:48) = 2 √ π β / (cid:18) √ x − β √ x (cid:19) e − βx . (24)Changing variables x = t/β , we write y ( t/β ) = q ( x ). The equation for y reads, y + (cid:18) t + 12 (cid:19) ˙ y + t ¨ y = 2 √ π β (cid:18) √ t − √ t (cid:19) e − t . (25)This is an inhomogeneous confluent hypergeometric equation (a.k.a. Kummer’s equation). It is easily verified that √ te − t is a solution of the homogeneous equation. The second solution involves the imaginary error function whichdiverges at infinity and is therefore not normalizable. Following the variation of parameters approach, we look for asolution in the form of y ( t ) = √ te − t [1 + βz ( t )] . (26)Substituting into (25), z ( t ) satisfies (cid:18) − t (cid:19) ˙ z + t ¨ z = 2 √ π (cid:18) t − (cid:19) . (27)Denoting w = ˙ z , t ˙ w + t (cid:18) − t (cid:19) w = 2 √ π (1 − t ) . (28)We split the right hand side into two cases t ˙ w + (cid:0) − t (cid:1) w = 1 t ˙ w + (cid:0) − t (cid:1) w = t (29)Using Matlab’s symbolic tool box, the solutions are w ( t ) = −√ πt − / e t erfc( √ t ) w ( t ) = − t − √ πt − / e t erfc( √ t ) (30)where, erfc is the complementary error function and the constants were chosen such that w and w vanish in thelimit t → ∞ . Overall, w ( t ) = 2 √ π t − t − / e t erfc( √ t ) . (31)Integrating, z ( t ) = 2 √ π ln t − (cid:90) t s − / e s erfc( √ s ) ds + C. (32)Substituting into y and then q , we obtain, q ( x ) = (cid:112) β √ xe − βx (cid:34) βC + 2 √ π β ln( βx ) − β (cid:90) βx s − / e s erfc( √ s ) ds (cid:35) . (33)The constant C is determined to enforce normalization.1 = (cid:90) ∞ p ( x ) dx = (cid:90) ∞ p ( x ) dx + (cid:15) (cid:90) ∞ q ( x ) dx + O ( (cid:15) ) . (34)As p is normalized, (cid:82) qdx = 0. Denote h ( y ) = y / e − y (cid:90) y s − / e s erfc( √ s ) ds. (35)2Expanding h for small y , the limit y → ∞ . Approxi-mating numerically, (cid:82) ∞ h ( y ) dy (cid:39) − . C (cid:39) √ π ( − . − γ − − β , (36)where γ is the Euler-Mascheroni constant. Substituting into (33), q ( x ) = (cid:40) β x = 0 p ( x ) (cid:104) D + ln( βx ) − √ π (cid:82) βx s − / e s erfc( √ s ) ds (cid:105) x > D (cid:39) − . q ( x ) is continuous and integrable ( (cid:82) ∞ q ( x ) dx = 0), the correction for theaverage infectivity rate (cid:82) ∞ x − q ( x ) dx , which also appears in (19), diverges close to x = 0 (because q (0) > (cid:15) , (cid:82) x − pdt = (cid:82) x − p dx + (cid:15) (cid:82) x − qdx + O ( (cid:15) ), asthe first integral is finite, while the second diverges.In order to overcome this difficulty, we need to take into account the finite population size. In our assumed quasisteady state, the population size is the current number of susceptible individuals, S . For simplicity, we apply a hardcutoff, and assume that the smallest value of x is P (1 / ( S + 1), where P ( x ) = (cid:82) x p ( t ) dt , the cumulative distributionfunction of p ( x ), (cid:90) m p ( x ) dx = 1 S + 1 . (38)Neglecting terms of order (cid:15) , the leading order term is m = β − S − / . Applying this cuttoff , we need to calculate (cid:90) ∞ m x − [ p ( x ) + (cid:15)q ( x )] dx. (39)Using known power law expansions of the error function, exponential integral ( (cid:82) t − e − t ) and lower incomplete gammafunction, (cid:90) ∞ m x − p ( x ) dx = 2 β + O ( S − / ) (40) (cid:90) ∞ m x − p ( x ) ln( βx ) dx = − γ + 2 ln 2) β + O ( S − / ) (41) (cid:90) ∞ m x − p ( x ) h ( βx ) dx = − √ π (cid:18) − γ + 23 ln S (cid:19) + O ( S − / ) (42)In particular, the integral over x − q ( x ) diverges logarithmically in S . Overall, (cid:90) ∞ p ( x ) dx = 2 β − β(cid:15) (1 . . S ) . (43)We now need to address the consistency of the approximation for the critical quasi steady state for which it is applied.Recall that we assume that, during the critical state, the fraction of infected individuals, i , is constant. Since therecovery rate is one, the duration of the critical state is at most 1 /i , i.e., T < /i . This duration is at least O (ln N ),even in exponential growth and decay regimes. Fig. 4a shows this assumption indeed holds. Therefore, i < C/ ln N .In SIR, the fraction of infected individuals also determines the rate at which susceptible individuals vanish (i.e.,become infected). Hence, in our case, ρ = i . Overall, we find that (cid:15) < C τ N . (44)Hence, the expansion is consistent as long as the proportionality constant C is small enough. For example, the slopeof a linear fit to Fig 4a yields C (cid:39) . k becomes infected is i ( t ) /µ k ( t ). Therefore, the population-averaged rate3is i ( t ) (cid:104) /µ k ( t ) (cid:105) S , where (cid:104)·(cid:105) S denotes averaging over the susceptible individuals. On average, the rate of change in s ( t )is dsdt = − (cid:104) /µ k ( t ) (cid:105) S s ( t ) i ( t ) . (45)Similarly, the population-average i satisfies didt = (cid:104) /µ k ( t ) (cid:105) S s ( i ) i ( t ) − i ( t ) = [ (cid:104) /µ k ( t ) (cid:105) S s ( i ) − i ( t ) . (46)Thus, a steady state would imply a critical reproduction rate, R t = (cid:104) /µ k ( t ) (cid:105) S s ( t ) = 1 . (47)We assume that s (cid:39)
1. As described above, the steady state distribution of µ k ( t ) has density p ( x ) = p ( x ) + (cid:15)q ( x ) + O ( (cid:15) ) with (cid:15) = iτ /
2. Hence, the reproduction rate at criticality is given by (43) R t (cid:39) (cid:90) ∞ x − p ( x ) dx = 2 β − βτ (1 . . N ) i critical . (48)Setting R t = 1 and β = α/ (cid:104) µ (cid:105) = 3 / (2 (cid:104) µ (cid:105) ) (recall that (cid:104) µ (cid:105) is the average infectivity time) yields, i critical = 4 (cid:18) − (cid:104) µ (cid:105) (cid:19) τ (1 . . N ) . (49)Of course, the critical i is positive, which sets a lower bound for (cid:104) µ (cid:105) , (cid:104) µ (cid:105) < . (50)In addition, recall our previous assumption that (cid:15) = O (1 / ln S ). Taking (cid:15) = iτ / (2 ln S ) (cid:39) δ and S (cid:39) N , we find thatthe expansion is consistent if i critical < δ/ ( τ ln S ). Comparing with (49), and assuming ln N (cid:29)
1, we obtain an upperbound for (cid:104) µ (cid:105) , (cid:104) µ (cid:105) > . δ. (51)Overall, we find a consistency condition for the existence of the steady state at α = 3 / . δ < (cid:104) µ (cid:105) < , (52)Comparing with the α = 1 . . < (cid:104) µ (cid:105) < .
7, we see that δ (cid:39) . α = 2. Following the same protocol, we find that p ( x ) = p ( x ) + (cid:15)q ( x ) + O ( (cid:15) ) p ( x ) = β xe − βx q ( x ) = p ( x )( γ − βx ) (53)Unlike the α = 3 / (cid:90) x − [ p ( x ) + (cid:15)q ( x )] dx = β (1 − (cid:15) ) . (54)As a result, assuming that s (cid:39) R t = (cid:104) /µ (cid:105) S = β (cid:18) − τ i (cid:19) = 1 . (55)We readily see that with β <
1, (58) cannot be satisfied. With β >
1, a solution may exist. However, i = O (1), whichimplies that the duration of the steady state is independent of N and therefore short.4Lastly, we take α = 5 /
2. We find that p ( x ) = p ( x ) + (cid:15)q ( x ) + O ( (cid:15) ) p ( x ) = 43 √ π β / x / e − βx q ( x ) = p ( x )( D + ln βx ) (56)where D (cid:39) − . (cid:90) x − [ p ( x ) + (cid:15)q ( x )] dx = β (1 − (cid:15) ) . (57)As a result, assuming that s (cid:39) R t = (cid:104) /µ (cid:105) S (cid:39) β (cid:18) − . (cid:15) (cid:19) = 1 . (58)With β < /
2, (58) cannot be satisfied. With β > /