Optimal periodic closure for minimizing risk in emerging disease outbreaks

Without vaccines and treatments, societies must rely on non-pharmaceutical intervention strategies to control the spread of emerging diseases such as COVID-19. Though complete lockdown is epidemiologically effective, because it eliminates infectious contacts, it comes with significant costs. Several recent studies have suggested that a plausible compromise strategy for minimizing epidemic risk is periodic closure, in which populations oscillate between wide-spread social restrictions and relaxation. However, no underlying theory has been proposed to predict and explain optimal closure periods as a function of epidemiological and social parameters. In this work we develop such an analytical theory for SEIR-like model diseases, showing how characteristic closure periods emerge that minimize the total outbreak, and increase predictably with the reproductive number and incubation periods of a disease– as long as both are within predictable limits. Using our approach we demonstrate a sweet-spot effect in which optimal periodic closure is maximally effective for diseases with similar incubation and recovery periods. Our results compare well to numerical simulations, including in COVID-19 models where infectivity and recovery show significant variation.


Introduction
The COVID19 pandemic, caused by the novel RNA virus SARS-CoV-2 [1], has resulted in devastating health, economic, and social consequences. In the absence of vaccines and treatments, non-pharmaceutical intervention (NPI) strategies have been adopted to varying degrees around the world. Given the nature of the virus transmission, NPI measures have effectively reduced human contacts-both slowing the pandemic, and minimizing the risk of local outbreaks [2,3]. The use of drastic NPI strategies in China reportedly reduced the basic reproductive number, R 0 , to a value smaller than 1, strongly curbing the epidemic within a short period of time [3,4]. On the other hand widespread testing protocols and contact tracing, in e.g., South Korea, significantly controlled spread during the initial phase of the pandemic [5]. In other countries, the implementation of NPI policies has not been as strict [2], with an optimistic reduction in transmission of roughly a half. To complicate the containment of the disease, early reports indicated significant amounts of pre-symptomatic and asymptomatic transmission [6,7]. For instance, recent estimates point to asymptomatic infection accounting for around 20-30% of the total, with a similar percentage for pre-symptomatic infections [8]-together producing a majority. These findings have been supported by other experimental studies [9] and analysis of the existing data [10,11].
As NPI controls such as quarantine, social distancing and testing are enforced, it is important to understand the impact of early release and relaxation of controls on the affected populations [12,13]. Recent studies have attempted to address how societies can vary social contacts optimally in time in order to maintain economic activity while controlling epidemics [14]. For instance, preliminary numerical studies suggest that periodic closure to control outbreak risk, where a population oscillates between 30-50 days of strict lockdown followed by 30-50 days of relaxed social restrictions, may efficiently contain the spread of COVID-19 and minimize economic damage [15]. These studies test interesting hypotheses, but cannot be immediately generalized to new emerging diseases. A basic understanding of why and when such risk minimizing strategies are effective remains unclear, and may benefit from a general analytical approach.
As a first step in this direction we analyze SEIR-like models with tunable periodic contact rates. Our methods reveal the existence of a characteristic optimal period of contact-breaking between individuals that minimizes the risk of observing a large outbreak, and predicts exactly how such an optimal period depends on epidemic and social parameters. In particular, we show that the optimal period for closure increases (or decreases) predictably with R 0 and the incubation period of a disease, and exists as long as R 0 is below a predictable threshold, and when there is not a time-scale separation between incubation and recovery. We demonstrate analytically that periodic closure is maximally effective for containing disease outbreaks when the typical incubation and recovery periods for a disease are similar-in such cases suppressing large outbreaks with R 0 's as large as 4. Our results compare well to numerical simulations and are robust to the inclusion of heterogeneous infection and recovery rates, which are known to be important for modeling COVID-19 dynamics.
To begin, we first consider the canonical SEIR model with a time-dependent infectious contact rate parameter, β(t). Individuals in this model are in one of four possible states: susceptible, exposed, infectious, and recovered. Following the simplest mass-action formulation of the disease dynamics, and assuming negligible background births and deaths, the fraction of susceptible (s), exposed (e), infectious (i), and recovered (r) individuals in a population satisfy the following differential equations in time (t), where dots denote time derivatives: Such equations are valid in in the limit of large, well-mixed populations and constitute a baseline description for the spreading of many diseases [16,17]. Note that α is the rate at which exposed individuals become infectious, while γ is the rate at which infected individuals recover. If β(t) = β 0 = constant, it is straightforward to show that the basic reproductive number for the SEIR model, R 0 , which measures the average number of new infections generated by a single infectious individual in a fully susceptible population, is R 0 = β 0 /γ [17][18][19]. Note in this work when R 0 is written as a constant (no time dependence) it should be taken to mean this value. Typical values for the R 0 of COVID-19 range from 1-4, depending on local population contact rates [4,20].

Methods
As a simple model for periodic closure we assume a step function for β(t) with infectious contacts occurring for a period of T days with rate β 0 , followed by no contacts for the same period, . A schematic of β(t) is plotted in the inlet panel of Fig  1(a). In S1 Appendix we show results for smoothly varying β(t) and asymmetric closure, where lockdown and open contacts occur for different amounts of time. It is demonstrated that the results presented in the main text do not qualitatively change under these generalizations. Also in Fig 1(a), we plot an example time-series of the infectious fraction, normalized by the initial fraction of non-susceptibles, for three different closure periods: green (short), blue (intermediate), and red (long). For periods that are not too long or short, the disease remains in a linear spreading regime (as we will show below), and therefore normalizing by the initial conditions gives time series that are initial-condition independent.
Intuitively, since the incubation period, α −1 , is finite, it takes time to build-up infection from small initial values. As a consequence, we expect that it may be possible to allow some disease exposure, before cutting contacts, and the result may be a net reduction in infection at the end of a closure period. For instance, notice that all i(t) decrease over a full closure cycle, 2T, in Fig 1(a). If the closure period is too small, infection can still grow (e.g., as T ! 0, R 0 (t) � hR 0 (t)i t = R 0 /2 which could be above the epidemic threshold), while if the period is too long, a large outbreak will occur before the control is applied. Between these two limits, there should be an optimal T (T min ), that results in a minimum outbreak. To illustrate, in Fig 1(b) we show an example of the final outbreak-size, r(t ! 1) � r f starting from i(t = 0) = 10 −3 , as a function of the closure period for different, equally spaced values of R 0 : the bottom curves correspond to smaller values of R 0 , while the top curves correspond to larger values.
As expected from the above intuitive argument, simulations show an optimal period that minimizes r f . A natural question is, how does T min depend on model parameters? Our approach in the following is to develop theory for T min in the SEIR-model, and then show how such a theory can be easily adapted to predict T min in more complete models, e.g., in COVID-19 models that include heterogeneous infectivity and asymptomatic spread [11,20].

Optimal control
It is possible to estimate T min by calculating its value in the linearized SEIR model, applicable when the fraction of non-susceptibles is relatively small. When e(t), i(t), r(t), 1 − s(t) � 1, the dynamics of Eqs (1)-(4) are effectively driven by a 2-dimensional system: The first step in calculating T min is to construct eigen-solutions of Eqs (5) and (6) in the form where ν(T) is the largest such eigenvalue; the superscript p denotes the corresponding principal eigenvector. Ignoring the subdominant eigenvalues assumes that after a sufficiently large number of iterations of periodic closure, the dynamics is well aligned with the principle solution no matter what the initial conditions. Unless stated otherwise, simulations are started in this state so that initial-condition effects are minimized. The second step is to calculate the integrated incidence, r(2T) from the solution of Eq (7), by integrating i(t) over a full cycle where [C p (t)] 2 denotes the infectious-component of C p (t). The third step is to calculate the final outbreak size from r(2T). To this end, it is important to realize that as long as ν(T) < 1, the outbreak will decrease geometrically after successive closure cycles, and therefore r f (T) = r(2T) + ν(T)r(2T) + ν(T) 2 r(2T) + . . ., or Finally, we can find the local minimum of r f (T) when ν(T) < 1 by solving This algorithm gives a single fixed-point equation that determines T min . Since our analysis is based on a piecewise 2-dimensional linear system, it is possible to give every quantity in the previous paragraph an exact expression [22] in terms of epidemiological and social parameters. See S1 Appendix for full derivation and exact expressions for Eqs (7)-(10). Following our procedure gives the prediction curves shown in Fig 2(a). The solid red line indicates the solution to Eq (10), and agrees well with simulation-determined minima of r f (T) over a range of R 0 given initial fractions of infectious 10 −6 (circles), 10 −4 (squares), and 10 −2 (diamonds). The simulation-determined minima are computed from r f (T) curves like Fig 1(b). It is important to note that our optimal-control theory assumes the validity of the linearized SEIR model, applicable when the total outbreak size, r f � 1 . In general, the total outbreak size will increase with the initial fraction of infectious and R 0 , and hence, the larger both are, the more simulations will disagree with theory. For example, this explains the better agreement for initial fractions of infectious 10 −6 , as compared to 10 −2 in Fig 2(b).
On the other hand, the solid blue line in Fig 2(a) indicates the threshold closure period, satisfying The closure period T thresh results in the largest eigenvalue of Eqs (5) and (6) equalling unity such that the principal component of exposed and infectious fractions is unchanged after a full closure cycle. If T < T thresh , ν(T) > 1 and a large outbreak occurs, even with closure, as infection grows over a full cycle for any small non-zero C(0). Given this property, T thresh gives a lower bound for the optimal period, T min > T thresh . Note: the red curve is always above the blue curve in Fig 2(a).
Before analyzing Eqs (5)-(10) further, we point out two basic dependencies in the (normalized) optimal period T min � γ. The first is intuitive: as the reproductive number R 0 increases, so does T min � γ. Hence, the faster a disease spreads the longer a population's closure-cycle must be in order to contain it. The second is more interesting. Notice in Fig 2(c) that T min � γ ! 1 as a ! 0, and T min � γ ! 0 as a ! 1. Therefore, recalling a = α/γ, if a disease has a long incubation period, then the optimal closure cycle is similarly long. On the other hand, if a disease has a short incubation period, then the optimal closure cycle is short. In order for periodic closure to be a practical strategy, with a finite T min , our results indicate that a � Oð1Þ, roughly speaking, or that the recovery and incubation periods should be on the same time scale-a condition that generally applies to acute infections [19].
Another observation from our approach that we can make is that periodic closure is not an effective strategy for arbitrarily large R 0 , as one might expect. One way to see this from the analysis is to notice that the optimal period diverges for the linear system at some R max 0 , as T thresh ! T min ! 1 (at fixed a). This transition can be seen in Fig 2(a), as the blue and red curves collide. Above the transition R 0 > R max 0 , no periodic closure can keep a disease from growing over a cycle. In this sense R max 0 ðaÞ gives an upper bound on contact rates between individuals that can be suppressed by periodic-closure as a control strategy. We note that an optimal T min still exists even when our linear approximation no longer applies, e.g., R 0 > R max 0 (in the sense that r(t ! 1) is minimized by some T min ), but the benefit of control becomes smaller and smaller as R 0 is increased, and the optimal period becomes increasingly dependent on initial conditions. In such cases, one must resort to numerical simulations of the full non-linear system, Eqs (1)-(4).
A sharper analytical understanding can be found by making the additional approximation that C(t) � exp[λ 11 γt]v 11 , for t < T and β(t) = β 0 , where ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ða þ 1Þ 2 þ 4aðR 0 À 1Þ Eq (12) is the largest eigenvalue of M(t<T) with eigenvector v 11 . Hence, we ignore the time-decaying part, C(t) dec � exp[−(a + 1 + λ 11 )γt]v 12 , of a general solution where v 12 is the other eigenvector of M(t<T). Our assumption becomes increasingly accurate with increasing T, and Eqs (7)-(11) simplify significantly: nðTÞ � e Tgl 11 ½ fe À Tg þ ð1 À f Þe À Tag �; ð13Þ rð2TÞ � r � e Tgl 11 À 1 l 11 þ e Tgl 11 1 À a and � r is a constant that depends on β 0 , α, γ and initial conditions, but is independent of T. Substituting Eqs (13)-(15) into Eqs (10) and (11) gives a single fixed-point equation for the approximate T min and T thresh each, which can be easily solved. See S1 Appendix for further details. Examples of the approximate solutions are plotted with dotted and dashed lines in Fig  2, and are almost indistinguishable from the complete linear-theory predictions shown with solid lines. Using the simplified expressions, we can now show several interesting features of periodic closure. First, since Eqs (13) and (14) are exact for large T, we can determine R max 0 as a function of a. As T ! 1, Eq (13) has two scaling limits depending on whether a � 1 or a < 1. In the former, the second term on the RHS of Eq (13) becomes negligible. As T ! 1 the solution of ν = 1 is λ 11 ! 1. Solving for R 0 in λ 11 = 1 gives R max 0 . Similarly when a < 1, as T ! 1 the solution of ν = 1 is λ 11 ! a. Putting the two cases together, gives R max 0 ðaÞ, and the phase-diagram for optimal-periodic closure: Eq (16) is plotted in Fig 3. In region I, the optimal period is predicted to be finite, in which case small outbreaks can be contained by optimal closure. In region II, such outbreaks can not be contained. The blue squares plot numerically-determined thresholds for the piecewise linear system Eqs (5)-(7) in the long closure-time limit. We compute each point by: picking a fixed value of a (starting with R 0 = 2), solving for T thresh according to Eqs (5)- (7) and (11), and then repeatedly incrementing R 0 in small steps of 0.001 and solving for T thresh (R 0 ; a) until it is a large number, i.e., T thresh (R 0 , a) � = 500. Note that as long as the system Eqs (1)-(4) is below threshold, we can always start with initial fractions of infectious and exposed that are small enough for the linear system to apply.
There are several important cases to notice in Fig 3. The first is that R max 0 has a peak when a = 1 (α = γ). The implication is that periodic closure has the largest range of effectiveness, as measured by the ability to keep infection from growing over any closure-cycle, for diseases with equal exposure and recovery times. In this symmetric case, periodic closure can prevent large outbreaks as long as R 0 < 4 (compare this to the usual epidemic threshold without closure, R 0 = 1). On the other hand, when there is a time-scale separation between incubation and recovery, a ! 1 or a ! 0, the phase-diagram nicely reproduces the intuitive, time-averaged effective epidemic threshold hR 0 (t)i t = 1, or R max 0 ¼ 2.

COVID-19 model
Now we turn our attention to more complete models that derive from the basic SEIR-model assumptions, but have more disease classes and free parameters which are necessary for accurate predictions. In particular, epidemiological predictions for COVID-19 seem to require an asymptomatic disease state, i.e., a group of people capable of spreading the disease without documented symptoms. Such asymptomatic transmission is thought to be a significant driver for the worldwide distribution of the disease [23,24], since symptomatic individuals can be easily identified for quarantining while asymptomatics cannot (without widespread testing). Many models have been proposed to incorporate the broad spectrum of COVID-19 symptoms, as well as control strategies such as testing-plus-quarantining [11,20]. A common feature of such models is the assumption that exposed individuals enter into one of several possible infectious states according to a prescribed probability distribution (e.g., asymptomatic, mild, severe, tested-and-infectious, etc.) with their own characteristic infection rates and recovery times. Following this general prescription, we define M infectious classes, i m , where m 2 {1, 2, . . .M}, each with its own infectious contact rate β m (t) and recovery γ m rate, and which appear from the exposed state with probabilities p m . The relevant heterogeneous SEIRmodel equations become Taking a common closure cycle for all individuals in the population, β m (t) = β 0,m � mod (floor{[t + T]/T}, 2) [21], we would like to test our method for predicting T min in the more general model Eqs (17) and (18), and demonstrate robustness to heterogeneity. In terms of an algorithm, we could simply repeat our approach for the effective 1 + M dimensional linear system; though, we loose analytical tractability. On the other hand, because T min is well captured by a linear theory, which depends only on R 0 , a, and γ, we might guess that quantitative accuracy can be maintained for higher dimensional models such as Eqs (17) and (18) by swapping in suitable values for these parameters in our SEIR-model formulas above. This is analogous to the epidemic-threshold condition (R 0 = 1) being maintained in such models, as long as the correct value of R 0 is assumed.
The R 0 for Eqs (17) and (18) is easy to derive using standard methods [17,18], Note: the updated R 0 is simply an average over the reproductive numbers for each infectious class. Using this averaging pattern as a starting point, our approach is to substitute the average values of α/γ m and γ m , into Eqs (7)-(10), or Eqs (13)-(15) for approximate solutions. Namely, for the SEIR model we have an equation 0 = F(R 0 , a, γ, T min ), where F is a function that is determined from Eq (10). Our averaging approximation entails solving the same Eq (10) for T min , but with parameters given by Eqs (19)-(21).
We point out that this approximation is not arbitrary since in the limit of heterogeneous infectivity only, γ m = γ8m, one solution of Eqs (17) and (18) is i m (t) = p m i(t), where i(t) is the total fraction of the population infectious. In this case, the linearized system is still effectively 2-dimensional with parameters γ, α/γ, and R 0 , where R 0 is given by Eq (19). For this reason we expect our averaging approximation to be exact in the limit of heterogeneous infectivity only, and a good approximation when the variation in recovery rates is not too large.
Examples are shown in Fig 4, where each panel shows results for an M = 2 model in which asymptomatics are significantly more (a) and less (b) infectious than symptomatics [11]. Symptomatic infectives are denoted with the subscript 1 and asymptomatics with the subscript 2. The optimal closure period is plotted versus the fraction of asymptomatics, p 2 . Within each panel the different colors correspond to no variation in recovery rates (red), moderate variation (blue), and large variation (green). Simulation determined T min are shown with points and predictions from the averaging theory shown with solid lines. The initial conditions for simulations follow the SEIR model convention-parallel to the principal solution of Eq (7), C p (0)-except that the fraction in each infectious class is i m (0) = p m [C p ] 2 . The model parameters were chosen to match similar models [11,20], which were fit to multiple COVID-19 data sources. As expected, the agreement between theory and simulations ranges from excellent to fair depending on the heterogeneity in recovery rates.

Discussion
Fig 4 demonstrates that the optimal closure period for COVID-19 can depend significantly on the amount of asymptomatic spread, particularly if there is a large difference in infection rates compared to symptomatic cases. Since asymptomatic spread is difficult to measure directly, especially in the early stages of an emerging disease outbreak, it may be difficult to estimate the optimal control accurately enough for periodic closure to be an actionable strategy on its own. A possible solution is to deploy effective and widespread testing within a population, early, and capture the fraction of asymptomatic infections. In any case, if basic parameters are The solid lines are theoretical predictions and the points are simulation-determined minima for initial fractions of non-susceptibles 10 −5 . Each series has different recovery times: red (g À 1 1 ¼ 10 � days, g À 1 2 ¼ 10 � days), blue (g À 1 1 ¼ 12 � days, g À 1 2 ¼ 8 � days), and green (g À 1 1 ¼ 14 � days, g À 1 2 ¼ 7 � days). The incubation period is α −1 = 7 � days. (b) Decreased infectivity for asymptomatics. Model parameters are identical to (a) except β 2 = 1.5 � γ 2 .
https://doi.org/10.1371/journal.pone.0244706.g004 known for an emerging disease dynamics, periodic closure is very effective-producing large reductions in the final outbreak size (e.g. , Fig 1(b))-and can be predicted using our methods.
An additional component of population heterogeneity not treated in this work is age dependence, which is known to be particularly important for modelling the COVID-19 pandemic. When considering expanded models that include age compartments, various mixing mechanisms across age groups generate different reproductive rates of infection [25][26][27]. One extreme compartmented grouping is to decompose a population into young, middle aged, and seniors with age-dependent contact rates between groups, age-dependent recovery periods, and some modest age-dependence in incubation periods. Under weak inter-age mixing assumptions, the result is a system of equations similar to Eqs (17) and (18). As demonstrated in Sec.3.2, the emergence of an optimal periodic control depends primarily on R 0 and the mean incubation period, and persists in spite of population heterogeneity. Although our controls are based on mean epidemiological parameters, it is easy to see how such controls may be distributed across age-dependent groups, and/or spatial clusters. Thus, we expect the inclusion of age-dependent effects to quantitatively change the results presented, but leave our methodology and qualitative findings intact.
Finally, we should remark that in addition to the heterogeneity discussed, parameter fluctuations for COVID-19 spread can occur in space and time. In fact, noise in reporting, differences in local policies, and adherence to the various forms of intervention may cause drastic fluctuations in the local spreading parameters. Given these facts, the well-mixed nature of our model may be insufficient to provide accurate optimal-control predictions. In such cases, a meta-population or network framework may be more appropriate. Yet, the approach that we lay out can be naturally generalized to more accurate and heterogeneous contact-network models, particularly since SIR and SEIR model dynamics on random networks can be described by relatively low-dimensional dynamical systems [28][29][30][31], which could be analyzed using the methods described in Sec.2. For small levels of infection, the main contribution from contact heterogeneity is to increase the effective, network R 0 . Once the correct R 0 is assumed, however, we expect the network results to be similar to those presented here, though this is a subject for future study.

Conclusion
In conclusion, a main socio-economic issue with an emerging virus, in the absence of vaccines and treatments, is the enormous damage at all levels of a population. Here we considered a simple approach to model and control an emerging virus outbreak with a finite incubation period. We show that by tuning periodic control of social contact rates, there exists an optimal period that naturally minimizes the outbreak size of the disease, as long as the reproductive number is below a predictable threshold and there is not a time-scale separation between incubation and recovery. Our basic assumption for the existence of such an optimal control rests on early detection of the disease, in which non-susceptible populations are small. Such a basic assumption allows one to analytically predict the optimal period, and provide parameter regions in which an optimal control exists. While in general it has been suggested that periodic closure may help curb the spread of an infectious disease like COVID-19, the implementation of such measures has been, to the best of our knowledge, mostly based on observations of recovery periods and absence of new cases for a given period of time. In this paper, we provide a general formulation that can be utilized to rationally design optimal intervention release protocols. While we start from an SEIR model and expand to heterogeneous models that capture the basic dynamics of COVID-19, our theory can be generally applied to acute infections, with the caveat that recovery and incubation periods should be roughly on the same time scale.
Supporting information S1 Appendix. Optimal control analysis. Supporting calculations and derivation of the outbreak-minimizing periodic control for the SEIR model. Additional simulation and analysis for both smooth and asymmetric control. (PDF)