key: cord-0919360-0r41e75c
authors: Jacyna, G.; Thompson, J. R.; Koehler, M.; Slater, D. M.
title: Analytical Model of COVID-19 for lifting non-pharmaceutical interventions
date: 2020-08-04
journal: nan
DOI: 10.1101/2020.07.31.20166025
sha: da28d9bb7fc7bcc86543c2be2a04ffafba39344c
doc_id: 919360
cord_uid: 0r41e75c

In the present work, we outline a set of coarse-grain analytical models that can be used by decision-makers to bound the potential impact of the COVID-19 pandemic on specific communities with known or estimated social contact structure and to assess the effects of various non-pharmaceutical interventions on slowing the progression of disease spread. This work provides a multi-dimensional view of the problem by examining steady-state and dynamic disease spread using a network-based approach. In addition, Bayesian-based estimation procedures are used to provide a realistic assessment of the severity of outbreaks based on estimates of the average and instantaneous basic reproduction number R0.

Throughout this section we use precise definitions of certain 86 properties of epidemiological models based on Meyers et al., 87 Bettencourt et al. and Chowell et al. [9, 10, 11] . 88 • Transmissibility T is the average probability that an 89 infectious individual will transmit the disease to a suscep- Note that some studies refer to an effective reproductive number 104 R e . In the definitions above, R 0 depends on the degree distri-105 bution so there is no need to make this distinction. R 0 and R e 106 can be thought of as interchangeable terms. 107 

Here we outline the procedure for characterizing steady-state 109 disease spread on complex networks using bond percolation. 110 The analysis is based on the work of Pastor-Satorras and Vespig-111 nani [1] and Newman [2, 3, 4, 5] . Generalizations to the theory 112 are made to include both uniform and directed social distancing. 113 Given the transmission rate r i,j between node i and node j 114 of a network graph and the infection time τ , the transmissibility 115 is:

as δt → 0. Typically, r i,j = r, where r and τ are independent 117 random variables. The average transmissibility is then:

118 T = 1 − e −rτ P r (r)P τ (τ )dr dτ,

where P r (r) and P τ (τ ) are the respective probability density 119 functions (pdfs). For simplicity, it is assumed that P r (r) = 120 δ(r − r 0 ) and P τ (τ ) = δ(τ − τ 0 ) so that T = 1 − e −r 0 τ 0 .

For a randomly chosen vertex, let p k denote the probability 122 that this vertex has k edges. Define G 0 (x) as the generating 123 function for the degree distribution of this vertex:

This function is similar to the characteristic function; that is, 125 given a sum of N independent and identically-distributed (i. i. d.) 126 random variables, the generating function is G N 0 (x).

Three types of social networks are examined. The Erdös-

Renyi network has a Poisson degree distribution of the form:

where λ is the mean. The exponential network is used as a 130 proxy for an urban network and has a degree distribution of the 131 form: 

where α and κ are parameters and Li α (x) = ∞ k=1 x k /k α is the 137 poly-logarithm function. 138 An important result by Feld is that the degree distribution 139 of the first neighbor of a vertex is not the same as the degree 140 distribution of vertices as a whole [12] . There is a greater chance 141 that an edge will be connected to a vertex of high degree, in 142 fact, in direct proportion to its degree. Let q k denote the degree 143 distribution of a vertex at the end of a randomly chosen edge.

Then:

excluding the randomly-chosen edge, where z = k kp k . The 146 corresponding generating function for this distribution is:

where G 0 (x) is the derivative of G 0 (x).

The transmissibility along each edge is taken into account 149 by interpreting T in Eq. (2) as a probability (0 ≤ T ≤ 1). 150 The probability that m out k edges is active is binomially dis-151 tributed, so:

where, similarly, G 1 (x; T ) = G 1 (1 − T + xT ). 

where G 0 (1) = G 1 (1) = 1. A phase transition occurs when 193 c → ∞ or, from Eq. (13), when T c = 1/G 1 (1) which implies:

In addition to the average outbreak size c , the distribu- 

This is a fixed point problem. The average outbreak size not as-213 sociated with the pandemic can also be determined. Repeating 214 the derivation leading to Eq. (13) and noting that H 0 (1; T ) =

The risk of an individual infection is determined by not- the probability that a vertex of degree k is present. Define the 228 generating function F 0 (x) as:

8 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint Similarly, F 1 (x) = F 0 (x)/z. This is a generalization of the 230 generating functions G 0 (x) and G 1 (x). Since F 0 (1) = 1 and 231 F 1 (1) = 1, 1 − F 0 (1) is the probability that a randomly chosen 232 vertex has no active edges. From Eq. (12), it follow that [3]:

The average outbreak size prior to a pandemic with social 234 distancing is similar to Eq. (13):

A phase transition leading to a pandemic occurs at the criti- :

The average outbreak size not associated with the pandemic is 244 analogous to Eq. (16) with G 0 (x) and G 1 (x) replaced by F 0 (x) 245 and F 1 (x), respectively. Noting that F 0 (v; T ) = F 0 (1) − P ,

.

The risk of an infection during a pandemic is identical to Eq. (17) All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . the effective reproduction number is R 0 = T /T c = bR 0 . This means that the effective reproduction number R 0 decreases with 260 increased uniform social distancing. The outbreak size distri-261 bution similarly follows from Eq. (19).

The fraction of the population affected by the pandemic fol-263 lows directly from Eq.(21). Letting v = (v −1+b)/b and noting

where T eff = bT . This is identical to Eq. 

where Here, it is assumed that:

where all individuals are distanced with contact degree greater 275 than K max . This implies that:

and F 1 (x) = F 0 (x)/z. For a Poisson network:

For an exponential network:

10 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint For a power-law network:

The average outbreak size prior to a pandemic is given by to K max by:

It is interesting to note that for the power-law network, if a small The degree distribution p k for each network is calibrated to 295 have the same critical transmissibility T c . Since T c = 1/G 1 (1) 296 from Eq. (14), it follows that:

In this section we expand the procedure to characterize time- tices to be in specific disease states. This is typically a diffi-307 cult problem because it involves higher-order moments of prob-308 abilities that can only be approximated using moment-closure 309 11 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . that a vertex with degree k is susceptible (s k ) is simply:

Similarly, the probability that a vertex with degree k is infected

where γ = 1/T r is the recovery rate. Finally, the probability 329 that a vertex of degree k has recovered (r k ) is:

A general solution to Eqs. (32)-(34) is of the form [5]:

where s 0 = s k (0). Given u(0) = 1, Eq. (35) can be numerically 

where G 0 (x) is the moment generating function for p k . (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint 1 + (γ/β) log(u) − s 0 G 1 (u) = 0. If it is assumed that u ≈ 1, 340 then:

where the pandemic size P = w(∞) and

is the average recovery probability. Assuming β/γ 1, the 

where v (t) = k q k x k and q k = bq k . Letting w (t) = k q k r k , 358 it follows that dw /dt = γv (t). Substituting this expression 359 into Eq. (38) and simplifying:

where u(t) = exp(−β(w − w 0 )/γ) and s k1 are the fractional 361 number of susceptible individuals with k contacts remaining at 362 the end of the first period (t = t 1 ). Since

Finally, since dw /dt = −(γ/βu)du/dt and dw /dt = γv (t):

13 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint

The total average infection probability v T (t) needs to include 

In 372 particular:

where u 1 = u(t 1 ). These infected individuals recover at a rate 374 γ since, by definition, they are not in contact with each other.

So:

Note that at

, so continuity is maintained. 

For the third period when social distancing is lifted, the 

Let v (t) = k q k x k denote the average infection probabil-387 ity, where x k = n k − s k − r k . Now, s k has the form:

whereŝ k2 = bs k2 + (1 − b)s k1 , s k2 = s k1 u k 2 , and u 2 = u(t 2 ). The 

where w 0 = −(γ/β) log(u 2 )−b(γ/β) log(u 1 ). Combining results:

Note that v (t 2 ) = v (t 2 ), as required. The fractional number 393 of susceptibles is given by:

14 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. 

where η is the delay in testing after becoming infectious and, 400 as before, v (t) = k q k x k . Rather than assuming that infec- immunity to the disease.

A general solution to Eq. (50) can be determined as follows.

Integrating the first expression produces s k (t) =ŝ k3 u k , where:

Eq. (51), it follows that:

Averaging the second expression in Eq. (50):

(53) Noting that G 1 (x) = k q k kx k−1 :

Therefore, Eq. (53) can be rewritten as: 

Note that Eq. (55) is initialized by setting v (t 3 ) = v (t 3 ).

15 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. 

The average outbreak case probability C(t) is determined by 433 averaging Eq. (57):

where u(t) is given by Eq. (35) and k p k ks 0 u k = s 0 u(t)G 0 (u). 

Now, returning to Eq. (33), averaging, and integrating the re-438 sults from t to t + ξ, assuming that ξ 1:

where b ξ (t; T ) = exp γξ(s 0 T u(t)G 1 (u(t)) − 1) and T = β/γ. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint

Define the instantaneous basic reproduction number R 0 (t) 473 as R 0 (t) = T u(t)G 1 (u(t)), which follows from Eq. (60), so that cedures. Therefore, the density is approximated using sampling 487 procedures.

In order to illustrate the procedure, the following simple 489 problem is considered. Suppose one is required to evaluate the 490 N th moment of p(x|z):

Assume a proposal density q(x|z) that is relatively easy to sam- 

the N th moment can be approximated by weighting samples x (i) 497 from a proposal density q(x|z) by a set of importance weights 498w i .

The same procedure can be used to estimate the posterior density. In this case w i = p(R J(i) |∆N J )/q(R J(i) |∆N J ) for a suitably chosen proposal density q(R J |∆N J ). Now, using Bayes theorem:

18 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. .

where it is assumed that p(∆N J |∆N J−1 , R J ) = p(∆N J |∆N J−1 , R J ).

In addition, assume that R J is Markov so that p(R J |∆N J−1 ) = 501 p(R J |R J−1 )p(R J−1 |∆N J−1 ). From Eq. (69):

(70) The proposal density can be factored into the following form 503 using the product of conditional densities:

. Using the definition of w i above:

(71)

For a simple particle filter examined here, q(R

where R Cov(R,Ṙ) = q ξ 3 /3 ξ 2 /2

19 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint where q is a variance-like term. It then follows that p(R J |R J−1 ) 525 is a zero-mean Gaussian density with covariance Cov(R,Ṙ).

The likelihood function is computed from Eq. (64). There are outbreaks that can occur outside the main pan-559 demic cluster, although they are relatively small. This is not 560 predicted from SIR or Susceptible, Exposed, Infected, or Re- (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . For uniform social distancing, it is assumed that the prob-581 21 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint similar trends although its asymptotic values are smaller than 596 the Poisson network because the population is not well-mixed.

The power-law network shows little variation in its pandemic 598 onset or its peak (at R 0 = 6) for a 20% increase in social dis-599 tancing from the baseline. Even at 40% social distancing, the 600 percent affected population is only somewhat reduced, although 601 its pandemic onset is shifted to a larger R 0 value. As explained 602 previously, this is due to a minority of super-spreaders that add 603 a degree of robustness to the network (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint high social contact individuals, the affected population is reduced from 75% to 50% for an exponential network and 30% 616 to zero for a power-law network when R 0 = 6. The affected 617 population for a Poisson network is reduced by only 10%. Al-618 though not depicted in this figure, it can be shown that for 619 a 1.5% reduction in directed social distancing, there is a 25% 620 reduction in the affected population for a power-law network 621 when R 0 = 6. For a 2% reduction in directed social distancing, 622 none of the population is affected for R 0 ≤ 6. However, the ex- (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint period T r is 7 days, the inception and release of (uniform) social distancing starts at week 6 and ends at week 14, respectively, 635 and testing and contact tracing begins at week 16. Additionally, 636 80% of individuals are social distanced, the percent isolated per 637 week is 20%, and the percent traced per week is 30% per week.

Figures 9 and 10 depict the susceptible, infectious, and out-639 break case load probabilities for a Poisson and exponential net-640 work, respectively, for the example parameters outlined above.

The pandemic begins to build after week 2 and is arrested start- (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . Figure 10 . Susceptible, infectious, and outbreak case load probabilities as a function of time for an exponential network based on example parameters 26 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint network for Massachusetts and an exponential network for New York. The estimated transmissibility for Massachusetts is 0.089, 685 which translates into a R 0 value of 1.8. Similarly, the estimated 686 transmissibility for New York is 0.185, which translates into a 687 R 0 value of 3.8. Both networks provide reasonable fits to their 688 respective observed differential outbreak case histories ∆C(t) 689 and reinforces the notion that the rapidity of disease spread in 690 New York was much more severe. The instantaneous basic reproduction number R 0 (t) can also 692 be estimated by tracking the differential outbreak cases over 693 time using a particle filter. Figures 14 and 15 illustrate the 694 estimation technique for New York given the collected outbreak 695 28 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint history. The fit to the observed differential outbreak case data 696 ∆C(t) is quite good. In addition, the instantaneous R 0 (t) values 697 are consistent with the average R 0 value using an exponential 698 network model depicted in Fig. 13. However, Fig. 15 is more 699 illustrative because it allows one to examine the trend in R 0 700 over the progression of the pandemic. In this case, there is 701 a downward trend after reaching a peak of approximately 3.5.

There is also an up tick in R 0 at later times probability due to 703 an increase in testing. Note that estimation errors may result 704 in negative R(t) values that are not realistic. Figure 14 . Estimated and predicted change in outbreak cases ∆C(t) based on New York State outbreak data using a particle filter Figure 15 . Estimated instantaneous basic reproduction number R 0 based on New York State outbreak data using a particle filter 705 29 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 4, 2020. . https://doi.org/10.1101/2020.07.31.20166025 doi: medRxiv preprint

In the present work, we developed and illustrated a coarse-707 grain analytic modeling procedure that can be used to assess 

Epidemic spreading in 738 scale-free networks

Spread of epidemic disease on networks

The structure and function of complex 743 networks

Network robustness and fragility: Percolation on random 746 graphs

Networks: An Introduction

Principles 750 of condensed matter physics

Dynamical patterns of epidemic outbreaks in com-754 plex heterogeneous networks

Network theory and sars: pre-760 dicting outbreak diversity

Real time bayesian es-763 timation of the epidemic potential of emerging infectious 764 diseases

Model parameters 767 and outbreak control for sars

Why your friends have more friends than you do

Conse-772 quences of delays and imperfect implementation of isolation 773 in epidemic control

Generalized poisson regression 775 model

Beyond the kalman 778 filter: Particle filters for tracking applications