Achieving herd immunity against COVID-19 at the country level by the exit strategy of a phased lift of control

The COVID-19 pandemic has affected the entire world causing substantial numbers of cases and deaths in most countries. Many have implemented nationwide stringent control to avoid overburdening the health care system. This has paralyzed economic and social activities and may continue to do so until the large-scale availability of a vaccine. We propose an alternative exit strategy to develop herd immunity in a predictable and controllable way: a phased lift of control. This means that successive parts of the country (e.g. provinces) stop stringent control, and COVID-19-related IC admissions are distributed over the country as a whole. Importantly, vulnerable individuals need to be shielded until herd immunity has developed in their area. We explore the characteristics and duration of this strategy using a novel individual-based model for geographically stratified transmission of COVID-19 in a country. The model predicts that individuals will have to experience stringent control for about 14 months on average, but this duration may be almost halved by further developments (more IC beds, better treatments). Clearly, implementation of this strategy would have a profound impact on individuals and society, and should therefore be considered carefully by various other disciplines (e.g. health systems, ethics, economics) before actual implementation.


Introduction
This document provides a technical description of the geographically stratified SEIR model used to predict the transmission dynamics of Covid-19 in a country or part of a country. The source code for the model in the form of the R package virsim is available online a under the CC BY-NC-ND 4.0 license. b At the end of the document, Table S1.1 provides a complete overview of the parameter values used in the main analysis. The implications of these parameters (e.g. distributions of local population sizes) are visualized throughout the document as part of the technical description.

Model description
We describe the dynamics of Covid-19 transmission using a standard SEIR model for a closed population of size , ignoring births and deaths. In terms of ordinary differential equations this can be described as: Here, is the force of infection, is the average contact rate in the population (including the average probability that transmission occurs during an average contact), is one over the average incubation time, and is one over the average duration of infectiousness, assuming exponentially distributed sojourn times. To relax assumptions about exponentially distributed durations and capture various other heterogeneities (more details below), we implemented the SEIR model in an individual-based framework in discrete time (one-day time steps). We assume that the durations of a http://www.gitlab.com/luccoffeng/virsim b https://creativecommons.org/licenses/by-nc-nd/4.0/ compartments and each follow a Weibull distribution with mean 0 and 1 and shape 0 and 1 ( Figure S1.1). Infection events (transitions from to ) are assumed to follow an exponential distribution, with the probability of an individual being infected on day defined as: Pr 5,7→0 = 1 − exp (−Δ • 5 ) Figure S1.1. Weibull distributions for duration of compartments and . Duration of compartment was assumed to follow a Weibull distribution with mean 5.5 and shape 20. Duration of compartment was assumed to follow a Weibull distribution with mean 10 and shape 0.8. The red dashed line indicates the arithmetic mean durations; dashed black lines indicate the symmetric 95%-confidence intervals.
To capture differential mixing of individuals in the same community (e.g. a town, ward or village) and administrative unit (e.g. a province), we distribute all individuals over superclusters that each consist of E individuals. Within each supercluster, we further subdivide the E individuals over E clusters of GE individuals. We allow for variation in cluster size GE by distributing the population over all ∑ E I clusters using a multinomial distribution with cluster-specific probability weights drawn from a log-normal distribution with mean 0 and standard deviation ( Figure S1.2). We assume that each supercluster contains the same number of clusters (i.e. unrelated to individual cluster sizes), so E = , for all = 0, … , .
To capture heterogeneity in contact rates of individuals and potential assortative mixing of individuals with similar transmission-related behavior, we assign each individual a life-long weight OGE which represents the individual's contact rate relative to the population average contact rate . We allow relative contact rates to vary between individuals according to a gamma distribution with equal shape and rate such that average relative contact rate is one ( Figure S1.3). Relative contact rates capture inter-individual variation resulting from both contact frequency and the probability of transmission per contact. Figure S1.2. Variation in cluster population size. The red dashed line indicates the arithmetic mean cluster population size (note that the horizontal axis is logarithmic); dashed black lines indicate the symmetric 95%-confidence interval. The histogram is based on ten thousand simulated cluster sizes with an expected size of one thousand population and = . . Figure S1.3. Variation in individual relative contact rate. Relative contact rates are assumed to follow a gamma distribution with equal shape and rate ( = . ) and mean one (red dashed line). Dashed black lines indicate the symmetric 95%confidence interval.
Assortative mixing (i.e. differential mixing of individuals with more similar contact rates) is captured by drawing and assigning individual relative contact rates OGE as follows. First, after assigning individuals to clusters, for each individual we draw a random value O from the unit normal distribution (0,1). Likewise, for each cluster we draw a random value G~( 0,1). For each individual in each cluster we then add up OG = O + • G , where is a parameter between 0 and 1 representing the level of assortative mixing. We then determine the rank of each individual in the entire population based on their value OG . Next, we draw values of OGE~Γ ( , ) and order these (ascending or descending order does not matter). Then finally, we assign each individual the n th of the ordered values of OGE , where n is the individual's rank in terms of OG ( Figure S1.4). If = 0, there is no assortative mixing and cluster-level average contact rates are all 1 (± Monte Carlo sampling variation). If = 1, there is maximum assortative mixing such that individuals' relative contact rates OGE within a cluster are extremely similar and the cluster-level average contact rates follow a gamma distribution Γ( , ). To account for the impact of heterogeneity in individual contact rates, geographical mixing patterns, and potential assortative mixing on transmission, we define the force of infection OGE acting on a susceptible individual in cluster in supercluster as: Here, It OGE u is an indicator function that returns 1 if individual in cluster in supercluster is in the compartment (i.e. infectious), and 0 otherwise. Interventions aimed at reducing contact rates (e.g. social distancing) are assumed to be implemented at the level of superclusters; their effect OE on contact rates may vary between individuals, with its square following a Beta distribution with mean and size (i.e. the sum of the distribution's shape parameters). In case the reduction of contact rates is the same across and within all superclusters ( E = and → ∞, such that OE = √ ), the quantity OE v = represents the reduction in the overall contact rate (i.e. the quantity reported in the main manuscript). However, to capture the effect of implementing control in only part of the population (i.e. when E ≠ ), we define E at the supercluster level. Changes in E over time can be specified per supercluster, allowing the simulation of a geographically heterogeneous intervention. In case E changes over time, individual reductions OE in contact rates are assumed to change proportionally to • E , and to be stable over time otherwise (reflecting the individual's inclination to adhere to control). Although not described above, the model also includes a mechanism that allows the user to specify which (random) fraction of the population will take up the intervention at each time point (default value 100%). For individuals who do not take up the intervention, we assume OE = 1.
Differential mixing of populations in clusters and superclusters is captured by mixing weights for population-level transmission ( ), supercluster-level transmission ( SC ), and cluster-level transmission (1 − SC − ). Isolation of a supercluster for control of transmission is simulated by multiplying a supercluster's contribution and exposure to population-wide transmission by (range 0-1).
In absence of variation in individual contact rates (i.e. OGE = 1) and in case of homogeneous mixing of the population ( = 1, and SC = 0) and uniform implementation of (and adherence to) control measures (i.e. OE = √ ), the above equation for OGE can be reduced to the original formulation of the force of infection in a simple SEIR model:  0.1717, such that the initial exponential growth of the epidemic is equal to one predicted by a homogeneous model with = 0.25, 0 = 5.5, and 1 = 10 ( " = 2.5) Shape and rate of gamma distribution for variation in individual relative contact rates (relative to overall contact rate ).
3.4, such that the 2.5 th and 97.5 th percentiles of the distribution of relative contact rates differ by a factor 10. Level of assortative mixing (range 0-1). 0.45, such that the 2.5 th and 97.5 th percentiles of the distribution of average relative contact rate in each cluster differ by a factor 3.8.