Superspreaders and High Variance Infectious Diseases

Superspreaders and High Variance Infectious Diseases Yaron Oz 1, Ittai Rubinstein 2, and Muli Safra 2 1 Raymond and Beverly Sackler School of Physics and Astronomy, Tel-Aviv University, Tel-Aviv 69978, Israel 2 Blavatnik School of Computer Science, Tel-Aviv University, Tel-Aviv 69978, Israel (Dated: September 6, 2020) A well-known characteristic of pandemics such as COVID-19 is the high level of transmission heterogeneity in the infection spread: not all infected individuals spread the disease at the same rate and some individuals (superspreaders) are responsible for most of the infections. To quantify this phenomenon requires the analysis of the effect of the variance and higher moments of the infection distribution. Working in the framework of stochastic branching processes, we derive an approximate analytical formula for the probability of an outbreak in the high variance regime of the infection distribution, verify it numerically and analyze its regime of validity in various examples. We show that it is possible for an outbreak not to occur in the high variance regime even when the basic reproduction number R0 is larger than one and discuss the implications of our results for COVID-19 and other pandemics. PACS numbers: 87.10.+e I. INTRODUCTION The classic SIR models provide an epidemiology frame- work for studying the spread of a disease [1]. The ba- sic reproduction number R0 in these models is the mean value of secondary infections caused by an infected indi- vidual. It determines the threshold R0 > 1 for an out- break. Alternatively, it determines the fraction of the population that will be infected before herd immunity is reached. In view of the importance of this parameter, major measures (such as lockdowns) are taken in order to reduce the value ofR0. The estimation for the COVID-19 pandemic, for example, is R0∼ 2− 3. The structure underlying the epidemic spreading is that of a complex heterogeneous network, where a small number of the nodes act as hubs while the majority of nodes have few contacts (for a review and references therein see [2]). Indeed, not all people cause a similar number of secondary infections and there is clear empiri- cal evidence for high levels of transmission heterogeneity in the infection spread (see e.g. [3–5]). The analysis in [4] for the COVID-19 pandemic suggests that between 5% to 10% of infected individuals are responsible for 80% of secondary infections. This may be due to differences in the number of contacts, in protective equipment, in levels of hygiene, in time of diagnosis or biological effects such as tendency to cough and sneeze. Individuals with high secondary infection rate are com- monly referred to as superspreaders. This is encoded in the degree distribution of the epidemic spread network corresponding to the infection distribution. While homo- geneous random networks such as the Erdos-Renyi model exhibit a statistical homogeneity of the nodes and the degree distribution is peaked around the average value, heterogeneous networks such as the scale free models re- veal a power law structure of the degree distribution and nodes with very large degree. The infection distribution is taken not over a random individual, but rather over a random infected individual, i.e. it is weighted according to the a priori probability of each individual to be infected. For instance, an individual in contact with many people has a higher likelihood both to be infected and to infect others and this is reflected in the degree of the corresponding node in the network. Studying the phenomenon of superspreaders, which seems to follow the Pareto-type Principle [6], as well as its implications on the spread of the disease is crucial when devising and implementing control policies [6, 7]. In order to analyze the impact of the superspreaders on the epidemic spread we have to consider the effect of the variance and the higher moments of the infection dis- tribution. The main goal of this paper, is to study a question of utmost importance when facing pandemics such as COVID-19, namely:“what is the probability that a disease will disappear without a major outbreak?” An outbreak is often referred to as a sudden rise in the number of infected individuals. In this paper, how- ever, we define an outbreak with reference to the total fraction of infected individuals in the long term and not at any specific point in time. Thus, we consider that an outbreak has not occurred if the disease has disappeared with a negligible herd immunity. Note, that we will anal- yse the natural evolution of the disease irrespective of the measures—social and others—taken to reduce R0. We will work in the framework of Galton-Watson branching processes (for a review see e.g. [8]), and use it to predict the probability of an outbreak as a func- tion of the infection distribution, that is the probability distribution for an individual to infect a given number of people. We derive an approximate analytical formula for the probability of an outbreak in the high variance regime of the infection distribution, verify it numerically in various examples, compare it to COVID-2 data and discuss its implications for the COVID-19 pandemic. In particular, we will show that it is possible for an outbreak not to occur in the high variance regime even when the basic reproduction number R0 is larger than one. This phenomenon has been observed in numerical simulations . CC-BY-NC 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 8, 2020. ; https://doi.org/10.1101/2020.09.06.20189258doi: medRxiv preprint NOTE: This preprint reports new research that has not been certified by peer review and should not be used to guide clinical practice. 2 [9]. II. THE HIGH V ARIANCE REGIME The infection distribution specifies, for each natural number k, the probability of an infected individual to infect k others. We denote by R0 and V the mean and variance of the number of people infected. When R0 < 1, it is well established that the disease would disappear on its own, while when R0− 1 is not small compared to the varianceV , one can use deterministic models such as SIR that provide an accurate description. Let us thus focus on the high variance regime: 0<R 0− 1≪V . (II.1) Our main result can be stated as follows. The probability that a disease will disappear without herd immunity is: Pr =γn, (II.2) wheren is the current number of infected individuals and γ in the regime (II.1) can be approximated as: γ≈ 1−Q , (II.3) where Q = 2(R0− 1) R2 0 +V−R0 . (II.4) Below, the corrections to the approximate formula (II.3) and (II.4) are bounded by higher powers of the ratioQ as well as the higher moments of the infection distribution. In section III we formulate the main result precisely and prove it. However, before delving into the proof let us consider some of its qualitative implications, compare it to pandemic data and numerically verify its accuracy. First, the larger the variance V compared to R0− 1, the higher the probability for the disease to disappear before herd immunity is reached. Thus, the fate of the disease does not depend only onR0−1. Second, the fewer infected individuals, the higher the probability for the disease to disappear and, consequently, the less stringent the pandemic measures that must be taken, even when R0 > 1. Third, the effective dependence on the variance is V n . Let us numerically compare our approximate analyti- cal formula to the exact γ for the re-scaled infection dis- tribution of COVID-2 [3]. The latter is based on fitting pandemic data to a distribution obtained by sampling a Poisson distribution whose mean is sampled from a Gamma distribution, which we will call Gamma-Poisson distribution. Since R0 of COVID-2 is high, we define the infection distribution for lower values of R0 by re-scaling the original one, that is, we fix the shape of the distribu- tion that is determined by a parameterk and re-scale the parameterθ that determines the scale of the distribution. FIG. 1: Our approximate formula 1 −Q vs. the exact Galton-Watson coefficient γ for the re-scaled infection distribution of COVID-2 [3]. The latter fits data to a distribution obtained by sampling a Poisson distribution with mean being sampled from a Gamma distribution. The parameter k≈ 0.19 for COVID-2 controls the shape of the Gamma distribution, while the parameter θ controls its scale. We constructed our data by fixing k and scale θ to give different values of R0 =kθ. We reach R0≈ 1.6 for the lower value of k. The results are depicted in figure 1 and, as expected, we see that lower R0 implies better accuracy. We use our formula to estimate the probability to avoid an outbreak for the COVID-2 and COVID-19 pan- demics as a function of R0 and the number of infected n—requiring an estimate of the ratio V R2 0 . Based on [3] we set V = 5R 2 0 for COVID-2. As noted above, the analysis in [4] estimates that the ph value (the percent- age of the infected population responsible for 80% of all secondary infections) is around 5% − 10%. Assuming Gamma-Poisson distribution, ph = 10% implies k = 0.1 and V = 10R2 0, and lower ph values correspond to even higher variance. We plot the results in figures 2 and 3: for given values of R0 andn, the higher the variance the higher the probability of avoiding an outbreak. While it is clear that γ cannot be determined precisely byR0 andV alone and the information about the higher moments of the infection distribution is necessary, our numerical analysis reveals that for certain distributions that are often being employed for real world pandemics the accuracy of (II.3) is mostly determined by the value of R0− 1 as depicted in figure 4 and figure 5. In figure 4, the Poisson distributions has λ values in the range 1.0 to 1.1, the geometric distributions has p values in the range 0.43 to 0.5, the Poisson10x distribu- tion is obtained by selecting a Poisson distribution with 10≤λ≤ 20 value with probability 10% or the zero dis- tribution with probability 90%, and the Truncated Power Law distributions has a cut-off at 100 with powers rang- ing from 2.1 to 2.375. In figure 5 we consider the ratio between the logarithms since this determines the ratio between the values of n that would give a specific prob- ability to avoid an outbreak. . CC-BY-NC 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 8, 2020. ; https://doi.org/10.1101/2020.09.06.20189258doi: medRxiv preprint 3 FIG. 2: The probability to avoid an outbreak when V = 5R2 0 (COVID-2) as a function of the basic reproduction number and the number of infected individuals. FIG. 3: The probability to avoid an outbreak when V = 10R2 0 (COVID-19) as a function of the basic reproduction number and the number of infected individuals. FIG. 4: The exact Galton-Watson coefficient γ vs. our approximate formula 1−Q for various infection distributions. The line y =x corresponds to γ = 1−Q. The accuracy of the formula depends mostly upon the value of R0, which explains the different deviations of the distributions from the y =x line: R0 for the Poisson, Geometric, Poisson10x and the truncated power law distributions are in the ranges [1, 1.1], [1, 1.32], [1, 2], [1, 1.72], respectively. III. FORMAL ST A TEMENTS AND PROOFS We define the infection distribution by a sequence of real variablesak, whereak is the probability that a carrier infects k individuals and is removed. The normalization condition is: ∑ k ak = 1. (III.1) Denote by Mi the ith moment of the infection distri- FIG. 5: A comparison of our formula to the exact value of γ for various distributions as a function of R0, and the larger R0 the larger the deviation. We tested our results for several distributions: Gamma-Poisson distributions for COVID-2 [3], Poisson and Truncated Power Law distributions. bution: Mi = ∑ k akki, M 1 =R0, M 2 =V +R2 0 , (III.2) and the quantity η by: η = 1 R2 0 +V−R0 ∑ i≥3 (3Q 2 )i−2 Mi i! . (III.3) Let p(x) be the polynomial p(x) = ∑ k akxk−x . (III.4) The proof of our result (II.2) and (II.3) consists of prov- ing three statements: • Pr =γn where 0≤γ 1 then the probability that the disease will disappear without herd immunity is γn whereγ satisfiesp(γ) = 0 and 0≤γ 0 s.t. if η < η0 then p(x) has a root within the interval [1− (1−cη)Q, 1− (1 +cη)Q] Claim III.3 (Single Root). p(x) has exactly one root in the interval [0, 1). Combining the above assertions, we see that if: . CC-BY-NC 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 8, 2020. ; https://doi.org/10.1101/2020.09.06.20189258doi: medRxiv preprint 4 • R0 > 1 (condition for Claim III. 1) • η <η0 (condition for Claim III. 2) •Q < 1 (1+cη0) (the root in Claim III. 2 is positive) then we arrive at our main result: 1−γ∈ [1−cη, 1 +cη]Q . (III.5) Claim III. 1 is a standard analysis of Galton-Watson processes [8], which we will now briefly review for com- pleteness. One views the number of sick individuals as a Markov process, where at each point we pick a sick individual, add the number of people infected by him and remove him. As above, ak is the transition probability from a state with n sick people to a state with k added infected people and one removed: n→n +k− 1 . (III.6) Denote byf(n) the probability that no major outbreak will occur at any timet> 0 if we have att = 0n infected people, and define γ =f(1)∈ [0, 1]. For the disease to die out, every branch that begins from one of the n infected individuals at t = 0 should disappear. Since we neglect the interaction between the infected individuals, these are independent random vari- ables and f(n) =f(1)n =γn. Using time independence and the total probability, one gets the recursion relation: f(n) = ∑ k akf(n +k− 1) . (III.7) Setting n = 1 in (III.7) we have: ∑ k akγk−γ = 0 , (III.8) that is, γ is a root of the polynomial p(x) (III.4). Finally, in order to complete the proof of Claim III. 1, we have to show that γ̸= 1. This is not surprising, as we are dealing with the R0 > 1 regime and setting γ = 1 would make the probability of an outbreak 1 − 1n = 0 regardless of the number of infected at t = 0. In order to prove the claim, we have to show that the probability of an outbreak converges to 1 as n→∞, but this is easy to see (for instance, by applying Chebyshev’s inequality on the probability that n sick will infect less than R0+1 2 n individuals). Consider next Claim III. 2. It is convenient to denote γ = 1 +δ and analyze the roots of p(x) : p(1 +δ) = ∑ k ak(1 +δ)k− (1 +δ) = 0 . (III.9) Expanding (III.9) and using (III.1) and (III.2) we get: p(1+δ) = (R0−1)δ+ 1 2 ( R2 0 +V−R0 ) δ2+ corrections. (III.10) From (III.10) we get the approximate formula (II.3) where the corrections are bounded by: corrections ≤ ∑ i≥3 Mi i! δi . (III.11) Let η0 = 1 10 and c = 5. We are interested in the case where δ∈ [−1−cη,−1 +cη]Q ⊆− [3Q 2 ,Q 2 ] . (III.12) Therefore: p(1 +δ) = = (R0− 1)δ + 1 2 ( R2 0 +V−R0 ) δ2± (R2 0 +V−R0)δ2η , (III.13) where we denote X = Y ±Z iff|X−Y| < Z. It is straightforward to see that when δ = −(1 + 5η)Q we have p(1 + δ) ≥ 0, while when δ = −(1− 5η)Q, we have p(1 + δ) ≤ 0. Combining these results with the Intermediate Value Theorem, we conclude the proof of Claim III. 2. In order to prove Claim III. 3, consider the second derivative ofp(x): p′′(x) = d2p(x) dx2 = ∑ k≥2 k(k− 1)xk−2 , (III.14) andp′′(x)> 0 forx> 0. Thus, p(x) is convex in R+, and must have at most two non-negative roots. Using (III.1) we see that x = 1 is one of these non-negative roots. Furthermore,x = 1 is not a local minimum of p(x), since p′(1) = ∑ k kak− 1 =R0− 1> 0 , (III.15) and in particular it cannot be the global minimum for p(x) in x ∈ R+. This implies that p(x) must have a negative value. ∀x > 1: p′′(x) > 0 implies that p′(x) > p′(1) > 0 and hence p(x) > 0. Therefore, p reaches its minimum in the R+ region at some point b, 0 ≤ b < 1. From the Intermediate Value Theorem, there is a point c, 0≤ c < b <1 such that p(c) = 0, and it is clearly unique, concluding our proof. IV. DISCUSSION AND OUTLOOK We have carried out an analysis of the stochastic spread of a disease in the high variance regime of the infection distribution. This allowed us to study an im- portant characteristic of the COVID-19 and other pan- demics where not all infected individuals spread the dis- ease at the same rate and superspreaders are responsible . CC-BY-NC 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 8, 2020. ; https://doi.org/10.1101/2020.09.06.20189258doi: medRxiv preprint 5 for most of the infections. We derived an approximate analytical formula (II.2 and II.3) for the probability to avoid an outbreak in the high variance regime (II.1) and estimated its accuracy numerically and analytically. We found out that R0− 1 is the main control parameter for the higher moment corrections. Curiously, for all the distributions that we analyzed we found that γ≤ 1−Q, giving us an upper bound on the approximation. We compared the formula to infection distribution data and discussed its implications for the COVID-2 and COVID- 19 pandemics. Our analysis reveals the general coarse-grained struc- ture of the infectious diseases irrespective of the detailed graph or network structure of the disease spread. We studied the natural evolution of the disease under the assumption that the infection and recovery are time- independent random variables. There are several rea- sons to consider the time dependence of R0, V and the higher moments, an obvious one being the measures, so- cial and other, taken to reduce them. A less obvious one is related to the time-dependent details of the disease’s evolution structure. There can be a major change due to a reduction in the number of superspreaders that are removed, which leads to interesting insights about the disease spread, such as reaching herd immunity faster than previously assumed [10]. Acknowledgements We would like to thank Nir Kalkstein for valuable dis- cussions on the importance of the high variance to the spread of the disease, as well as Baruch Barzel for com- ments on the manuscript. The work is supported in part by the Israeli Science Foundation center of excel- lence. The work is supported in part by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 835152), as well as by ISF 2013/17 and BSF 2016414. [1] W. O. Kermack and A. G. McKendrick, ”A contribution to the mathematical theory of epidemics”, Proceedings of the Royal Society A, Volume 115, Issue 772 (1927). [2] R. Pastor-Satorras, C. Castellano, P. Van Mieghem and A. Vespignani, ”Epidemic processes in complex net- works”, Rev. Mod. Phys. 87, 925 (2015). [3] J. O. Lloyd-Smith, S. J. Schreiber, P. E. Kopp and W. M. Getz, ”Superspreading and the effect of individual variation on disease emergence”, Nature 438, 355–359 (2005). [4] D. Miller et al.,”Full genome viral sequences inform pat- terns of SARS-CoV-2 spread into and within Israel”, doi:https://doi.org/10.1101/2020.05.21.20104521. [5] T. Britton, F. Ball, P. Trapman, ”A mathematical model reveals the influence of population heterogeneity on her immunity to SARS-Cov-2”, Science 10.1126/sci- ence.abc6810 (2020). [6] M. E. J. Woolhouse et al., ”Heterogeneities in the trans- mission of infectious agents: Implications for the design of control programs”, Proc. Natl Acad. Sci. USA 94, 338 (1997). [7] R. Pastor-Satorras and A. Vespignani, ”Immunization of complex networks”, Phys. Rev. E. 65, 036104 (2002). [8] D. R. Insua, F. Ruggeri, M. P. Wiper, ”Bayesian Analysis of Stochastic Process Models”, Wiley (2012). [9] B. Barzel, private communication. [10] Y. Oz, I. Rubinstein and M. Safra, to appear. . CC-BY-NC 4.0 International licenseIt is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 8, 2020. ; https://doi.org/10.1101/2020.09.06.20189258doi: medRxiv preprint