key: cord-0321029-wl8r1ohd
authors: Khataee, H.; Kibble, J.; Scheuring, I.; Czirok, A.; Neufeld, Z.
title: Transition from growth to decay of an epidemic due to lockdown
date: 2020-12-14
journal: nan
DOI: 10.1101/2020.12.14.20248154
sha: 7e25b01abe990cd88191672958654560b46264f1
doc_id: 321029
cord_uid: wl8r1ohd

We study the transition of an epidemic from growth phase to decay of the active infections in a population when lockdown measures are introduced to reduce the probability of disease transmission. While in the case of uniform lockdown a simple compartmental model would indicate instantaneous transition to decay of the epidemic, this is not the case when partially isolated active clusters remain with the potential to create a series of small outbreaks. We model this using a connected set of stochastic susceptible-infected-removed/recovered (SIR) models representing the locked-down majority population (where the reproduction number is less than one) weakly coupled to a large set of small clusters where the infection may propagate. We find that the presence of such active clusters can lead to slower than expected decay of the epidemic and significantly delayed onset of the decay phase. We study the relative contributions of these changes to the additional total infections caused by the active clusters within the locked-down population. We also demonstrate that limiting the size of the inevitable active clusters can be efficient in reducing their impact on the overall size of the epidemic outbreak.

Quantitative epidemiological models are being used as a tool to inform public health decisions about infectious disease outbreaks, like recently during the COVID-19 pandemic. Simpler mathematical models were useful to introduce concepts such as, the basic reproduction number R 0 [1] , while more complex agent based models were used to estimate hospital occupancy, or the effect of planned social distancing measures [2] [3] [4] . Arguably the simplest and most influential epidemic model operates with the susceptible, infected and removed/recovered subpopulations and hence referred to as the SIR model [5] . The corresponding system of nonlinear ordinary differential equations exhibits an initial exponential spread of an epidemic, followed by saturation as the susceptible proportion of the population decreases. Eventually, as the immune sub-population reaches a critical threshold, often referred to as herd-immunity threshold, the infected subpopulation decays to zero. The simplest SIR model assumes that the epidemic parameters (e.g., the probability of transmission from an infected to a susceptible individual) remain constant during the outbreak, and the whole population is perfectly mixed: any two individuals have the same chance for disease transmission. However, in the COVID-19 pandemic a range of public health measures were introduced with the goal of reducing transmission efficiency. Accordingly, the initial exponential growth phase of the infectious fraction of the population was followed by an exponential decay after the introduction of social distancing measures, primarily reflecting the reduced transmission rate well before a significant fraction of immune population would have accumulated [6] .

Detailed analysis of epidemic and mobility data from various European countries during the first wave of the COVID-19 pandemic has revealed that after the introduction of social distancing measures the transition from growth to decay phase may involve a fairly long, and country-specific delay [6] . The 4-5 days long incubation period of the disease [7, 8] , as well as the average survival time of patients with fatal COVID-19 illness are expected to introduce a certain delay between the alteration of social behavior and its observable effect in epidemic data such as the daily new infections and death rate. However, in most countries the delay between the introduction of control measures to the peak of the epidemic was significantly longer, usually lasting for several weeks. The same study also revealed that the delay was longer in countries where the control measures were less strict [6] .

Within the context of an SIR model a substantial sudden reduction of the infection transmission rate due to social distancing efforts should flip the exponential growth into exponential decay instantaneously. In this paper we demonstrate that the delay τ between the introduction of social distancing measures and the onset of the decaying phase of an epidemic can arise when control measures affect various groups of society highly selectively. For example, under official lockdown orders transmission can be reduced for the majority of the population, while the disease remains able to spread among various distinct clusters -like workers in warehouses, meat processing plants or among prison inmates. To study this effect we propose a generalization of the SIR model to take into account the heterogeneity of the population with respect to epidemic control measures. We consider a dual population structure where most of the population has a strongly reduced rate of transmission, while the transmission rate remains high within distinct clusters of a certain size. By stochastic simulations we explore how the delay τ , the decay rate and the total number of infected people are influenced by the size of the non-restricted sub-population as well as the typical size of the clusters within the non-restricted population.

To model the spread of the epidemic in a population we use the standard SIR epidemic model, defined by the following system of ODEs [5, 9] :

where I, S, and R denote the infected, susceptible, and recovered fractions of the population, respectively, i.e., S +I +R = 1. Here, k is the rate of transmission from I to S per infected individual, while γ represents the rate of recovery from I to R, either through acquired immunity or death.

The behavior of the model is determined by a single non-dimensional parameter, the reproduction number R 0 = k/γ, which gives the number of new infections caused by a single infected person within a fully susceptible population. The condition for an epidemic outbreak is R 0 > 1, while for R 0 < 1 the epidemic dies out.

In the initial phase of an epidemic described by the SIR model most of the population is susceptible and the number of infected individuals grows exponentially. We are interested in the epidemic dynamics within a population under social distancing measures that have been implemented sufficiently early (i.e., I, R 1 and S ≈ 1) and are efficient enough to reduce the reproduction number below one. In this regime, Equation (1) simplifies to dI/dt = kI − γI and yields:

where parameters R 01 > 1 and R 02 < 1 characterise the initial growth and post-intervention decay phases of I, respectively.

To model the heterogeneity of the population with respect to epidemic control measures, we divide the total population of size T into sub-populations: (i) the locked-down population fraction of size 0 < L < 1, where the social distancing measures are effective and the infection transmission rate reduces (R 0 < 1), and (ii) the rest of the population (e.g., nursing home inhabitants, prison inmates, essential service providers) are sorted into partially isolated clusters, and within the clusters the infection transmission rate remains unaffected (R 0 > 1). For simplicity, we assume that each cluster has the same size, M , a model parameter that reflects both the organization of social institutions as well as potential efforts to control the epidemic. The number of such clusters is N , and thus the size of the non-socially distancing sub-population is N M = (1 − L)T . Table 1 .

The locked-down sub-population and each cluster has its own SIR dynamics, which are also coupled as depicted in Fig. 1(A) . Specifically, we define the transmission rates as follows. Within the locked-down sub-population (L), we define k L→L = R 02 γ. Thus, in this sub-population, in the absence of external interactions the number of infections would decay exponentially. In contrast, for each cluster j = 1, .., N , we assume that the infection transmission rate within the cluster remains unchanged, i.e., k j→j = R 01 γ. Transmission rates are also reduced between clusters and the locked-down sub-population: k L→j = k j→L = R 02 γ. Finally, we also assume that there is no direct transmission of infection between different clusters.

We can now establish the dynamics for each of the SIR subsystems. For each cluster j, new infections occur at the rate given by:

where I j and S j are the proportions of infected and susceptible in the cluster, the first term represents new infections transmitted within the cluster, and the second term describes infections introduced from the locked-down population. The parameter 0 < β j < 1 represents the probability that an individual in a cluster j interacts with others within her own cluster (we use β j = 0.9). In the second term, the factor M/T captures the probability that an individual of the locked-down sub-population interacts with a member of cluster j.

Likewise, in the locked-down sub-population, the total rate of new infections is given by:

where β L = 1 − β j = 0.1 represents the probability of interaction between individuals of cluster j and the locked-down sub-population.

3 . CC-BY-ND 4.0 International license It 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 this this version posted December 14, 2020. ;  Now, the dynamics of the proportion of infected individuals in the locked-down sub-population (I L ) and in the active clusters (I j ) are determined by the normalised rates of new infections combined with the rate of recovery, Q L /LT − γI L and Q j /M − γI j , respectively. Note, while the SIR dynamics is the same for each cluster j, the temporal dynamics differs as the first in-cluster infection transmission could be initiated at various time points.

In the case of the standard SIR model it is natural to use a deterministic description, since in a large population the fluctuation of individual discrete events can be neglected and the overall dynamics is well approximated by continuous variables following mass action type kinetics. However, in our model, while the total population is still assumed to be large (for concreteness, we consider a city with total population T = 10 6 ), after implementing lockdown strategies, the size of partially isolated clusters (where the infection is still able to propagate) are relatively small; we will assume typically values in the range M = 20 − 150. Thus, a single new infected person in a cluster can increase I j by a much larger extent than a similar infection in the larger population increases I L . Therefore, the discreteness of the population and the stochastic nature of the epidemic dynamics is important, especially at the initiation of epidemic outbreaks within the clusters, triggered by a single random infection event.

We implemented the time evolution of the interacting system of sub-populations as a series of discrete events in continuous time, by using the Gillespie stochastic simulation algorithm [10] [11] [12] .

The parameter values used in the simulations are summarised in Table 1 . We used R 01 = 3 and R 02 = 0.6 since the estimations of R 01 and R 02 are in the ranges 2 − 3.5 [13] [14] [15] and 0.36 − 0.82 [6] , respectively. Although the exact value of γ for COVID-19 is uncertain [16, 17] , we used γ = 0.1 (day −1 ) as a reasonable estimate for this parameter [18] . Yet, our analysis of the epidemic dynamics is rather general and does not rely on the values of the model parameters. Stochastic trajectories of the total number of infected individuals are shown in Fig. 1(B) . The typical time course of the epidemic (I tot ) was obtained as an average over 10 runs of the algorithm for each set of parameters N and M ; see Fig. 1 (B). 

We first explore how the time course of the total number of infected individuals in the whole population, I tot (t), is affected by the number and size of the clusters which remain active after the lock-down. We find that raising the number N and size M of these clusters can lead to a dramatic increase in the number of infections in comparison to the reference case when the whole population is locked down; dashed green line in Fig. 2 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 this this version posted December 14, 2020. ; Figure 2 : Total number of infected individuals (I tot (t)) versus time for different numbers (N ) and sizes (M ) of active clusters in the not socially distancing sub-population. Circles: I tot (t) calculated using the Gillespie algorithm after averaging over 10 repeated simulations. Solid red line is obtained by fitting I 02 e −αt to the simulation data points highlighted by orange using the Mathematica routine NonlinearModelFit. To perform the fitting with minimal effects of the stochasticity, in most cases, the fitting range is set to start at a time after lockdown time t LD where 100 < I tot(t) < 3000 (see Table S1 for the fit parameters I 02 and α). Dashed green line: I LD e γ(R0 2 −1)(t−t LD ) that passes through I LD = I tot (t LD ), where t LD is the time of lockdown. Black dashed line: t LD . Blue line: I 01 e γ(R0 1 −1)t . Other simulation parameters are in Table 1 .

The significant differences observed in Fig. 2(A-I) in the number of infected individuals with varying the numbers and sizes of the active clusters in the population persuade us to further examine the transition to the decay of the epidemic. We characterise the transition by three properties: the decay rate of the number of infected individuals (α), the delay in the decay phase caused by the presence of the active clusters (τ ), and the additional number of infected individuals (I extra ) due to the slower and delayed decay of the epidemic.

The decay rate is calculated by fitting I tot (t) I 02 e −αt to the simulation data (averaged over 10 repeated realisations) within a selected range where there is a well defined exponential decay; see 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 this this version posted December 14, 2020. ;

Next, the delay τ of the onset of the exponential decay phase, caused by the not locked-down clusters, is calculated as the time-shift relative to decay with the same rate α passing through I LD = I tot (t LD ) ≈ 5000; see Fig. 3(A) . This decay follows I 02 e −α(t+τ ) (illustrated using dotted brown line in Fig. 3(A) ). Therefore, I LD = I 02 e −α(t LD +τ ) results in:

Finally, the number of new additional infected individuals, caused by the presence of not lockeddown clusters, is calculated as:

which corresponds to the area of the yellow shaded region in Fig. 3(A) , multiplied by γ. Note that the second integral, representing the area under the dashed green line in Fig. 3(A) , gives the number of individuals infected during the lockdown (t > t LD ) if a complete lockdown was implemented with no active clusters (N = 0). This can be calculated explicitly and has the value I LD /(1 − R 02 ), which we use as a reference value for the extra infections. Thus the effect of the clusters in the incomplete lockdown N > 0 can be characterised by the relative increase of total infections:

Results in Fig. 3(B) show that as the number and size of the active clusters increases the number of infected individuals decays slower than the value expected for a uniformly locked-down SIR model: α 0 = γ(1 − R 02 ) = 0.04 (day −1 ). For the range of values considered in the simulations we find that when the locked-down fraction decreases to around L = 0.8 the characteristic decay time of the epidemic, α −1 , increases from the reference value of 25 days (for L = 1) to around 50 days, i.e., α ≈ 0.02 (day −1 ).

Figure 3(C) shows that the delay of the decay phase, τ , can increase substantially in the range of tens to hundred days when the locked-down fraction is reduced to 80-90%. Also we find again that for a certain locked-down fraction, the delay can be much longer when the active population is composed of larger clusters. The relative increase of the total infections during the lockdown due to the active clusters (µ) shows a similar behavior (see Fig. 3 (D)) indicating up to more than 10-fold increase when the active proportion is increased to around 10% of the population.

6 . CC-BY-ND 4.0 International license It 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 this this version posted December 14, 2020. ; Table 1 .

To further illustrate the differential impact of the number and size of the active clusters on the extra infections during the decay phase of the epidemic, we generate a heat-map representing I extra (M, N ); see Fig. 4(A, B) . Using a logarithmic scale for both variables N and M , the linear iso-contours in Fig. 4(B) suggest to approximate the dependence on these variables by a power law of the form:

We find that this functional form leads to a good fit for the scaling of the data from the stochastic simulations ( Fig. 4(C, D) ) resulting in a ≈ 1 and b ≈ 1.8. This confirms our earlier prediction that increasing the size of the clusters (M ) has a stronger impact on the total infections than increasing the number of clusters (N ) by the same proportion.

7 . CC-BY-ND 4.0 International license It 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 this this version posted December 14, 2020. ; We further examine how the rate and delay of the decay phase contribute to the extra infections in the population. This relationship can be quantitatively expressed as:

where the two terms represent the contributions due to slower and later decay, respectively. The approximation is based on replacing the transition region to a constant segment I(t) = I LD until the onset of the decay phase (see S1(A)). Although this formula slightly underestimates the total infections in some cases, it provides a very simple representation of the two components.

We find a good agreement between the analytical and numerical calculations of the extra number of infected individuals; see Fig. S1 (B), and the increase of infected individuals appears to be mainly due to the delay (τ ) in the decay phase, rather than the rate (α) of the decay; compare Fig.  S1(C, D) . These results provide an explanation for our earlier prediction, that dividing the active sub-population into smaller clusters reduces the total number of infected individuals, mainly due to faster transition from growth phase to decay phase when the size of the outbreaks in infected clusters is smaller.

8 . CC-BY-ND 4.0 International license It 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 this this version posted December 14, 2020. ; https://doi.org/10.1101/2020.12.14.20248154 doi: medRxiv preprint

We studied the effect of a large number of small active clusters within a locked-down majority population on the transition to the decay phase of an epidemic. We have shown that the active clusters can lead to delayed onset of the decay of the number of infected in the population and a slower decay rate of the epidemic. In our representative example simulations, considering a COVID-19-like epidemic in a city of one million inhabitants, we find that the delay can be in the range of several weeks to a few months while the characteristic time of the exponential decay, assumed to be around 3 weeks for the locked-down population, may double. This could lead to up to tenfold increase in the total number of infections occurring from the start of the lockdown relative to the case of a complete uniform lockdown. We found that the size of the clusters has a much stronger impact on the total extra infections in the incomplete lockdown than the number of these clusters, and the extra infections can be approximated with a double power law expression of the form: I extra ∼ N · M 1.8 . We also identified that the dominant factor leading to the extra infections is the delayed transition to the decay phase which has an approximately 5-fold larger contribution than the slower than expected decay.

The above results are based on a series of simplifying assumptions which allow us to keep a fairly simple general mathematical framework and efficient computer simulations while capturing the important impact of small minority groups in disease transmission, which is completely missed in standard ODE based compartmental epidemic models. This is particularly relevant in the context of the current COVID-19 pandemic, where unexpected strong over-representation of certain affected subgroups in the population become evident (e.g., very high proportion of nursing home deaths, or frequent outbreaks in meat factories, livestock plants, and prisons) [19] [20] [21] [22] . Our model of course gives a very crude representation of these subgroups, considering a single average size of the clusters and assuming that the reproductive ratio remains the same as in the initial uncontrolled epidemic. Also, for simplicity, we assumed that the outbreaks occurring within the clusters are only stopped by the development of natural herd immunity within the cluster. In reality these outbreaks may be detected before they reach herd immunity and can be controlled through quarantine measures. This could be implemented quite easily in the model and would likely be equivalent to considering a smaller effective cluster size in the current model.

We also note, that the disproportionate impact of small sub-populations in epidemic models have also been observed in various other contexts. For example, in social contact network based epidemic models, it was shown that a very small fraction of highly connected individuals can have a large effect on the spread of infections in the population [23] . In the case of COVID-19, it was suggested that certain super-spreading events may play an essential role in the epidemic dynamics [24] . These examples also demonstrate potential limitations of mean field type models restricted to representing the average properties and behavior of the interacting components. 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 this this version posted December 14, 2020. ; https://doi.org/10.1101/2020.12.14.20248154 doi: medRxiv preprint . CC-BY-ND 4.0 International license It 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 this this version posted December 14, 2020. ; https://doi.org/10.1101/2020.12.14.20248154 doi: medRxiv preprint

H. Khataee, J. Kibble, I. Scheuring, A. Czirok, Z. Neufeld Table S1 : Fit parameters I 02 and α are approximated by fitting I 02 e −αt to the calculations of the average total number of infected individuals (highlighted by orange in Fig. 2) . Parameter values correspond to mean ± standard error (SE). 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 this this version posted December 14, 2020. ; https://doi.org/10.1101/2020.12.14.20248154 doi: medRxiv preprint Figure S1 : Contributions to the extra infected individuals, I extra due to slower vs delayed decay of the epidemic. (A) A sketch illustrating areas associated with I extra , for example with N = 2000 and M = 70. Area A 1 (blue shaded region) represents extra infected due to the slower rate of decay of the infections; first term in Equation (9) . Area A 2 (green shaded region) represents an estimate of the extra infected caused by the delay in the decay phase; second term in Equation (9) . 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 this this version posted December 14, 2020. ; https://doi.org/10.1101/2020.12.14.20248154 doi: medRxiv preprint

Infectious diseases of humans: dynamics and control

Strategies for containing an emerging influenza pandemic in Southeast Asia

Report 9: impact of non-pharmaceutical interventions (npis) to reduce covid-19 mortality and healthcare demand

Individual-based simulation model for covid-19 transmission in daegu, korea

A contribution to the mathematical theory of epidemics

Effects of social distancing on the spreading of covid-19 inferred from mobile phone data

of Novel Coronavirus-Infected Pneumonia

The incubation period of coronavirus disease 2019 (covid-19) from publicly reported confirmed cases: estimation and application

Population biology of infectious diseases: Part I

A general method for numerically simulating the stochastic time evolution of coupled chemical reactions

Stochastic modelling for quantitative description of heterogeneous biological systems

Inverse gillespie for inferring stochastic reaction mechanisms from intermittent samples

Review of the clinical characteristics of coronavirus disease 2019 (covid-19)

Early dynamics of transmission and control of COVID-19: a mathematical modelling study

Unique epidemiological and clinical features of the emerging 2019 novel coronavirus pneumonia (covid-19) implicate special control measures

Epidemiology and transmission of COVID-19 in 391 cases and 1286 of their close contacts in Shenzhen, China: a retrospective cohort study

Temporal dynamics in viral shedding and transmissibility of COVID-19

The challenges of modeling and forecasting the spread of covid-19

Coronavirus in the us: latest map and case count

Sars-cov-2 outbreak investigation in a german meat processing plant

Livestock plants and covid-19 transmission

Evolution and effects of covid-19 outbreaks in care homes: a population analysis in 189 care homes in one geographical region of the uk

Epidemic spreading in scale-free networks

Sars-cov-2 (covid-19) superspreader events