key: cord-0638639-dc4n97lt
authors: Fujiwara, Naoya; Onaga, Tomokatsu; Wada, Takayuki; Takeuchi, Shouhei; Seto, Junji; Nakaya, Tomoki; Aihara, Kazuyuki
title: Analytical estimation of maximum fraction of infected individuals with one-shot non-pharmaceutical intervention in a hybrid epidemic model
date: 2021-03-03
journal: nan
DOI: nan
sha: a20518380464cb227ecf7b378a813514e063e546
doc_id: 638639
cord_uid: dc4n97lt

Facing a global epidemic of new infectious diseases such as COVID-19, non-pharmaceutical interventions (NPIs), which reduce transmission rates without medical actions, are being implemented around the world to mitigate spreads. One of the problems in assessing the effects of NPIs is that different NPIs have been implemented at different times based on the situation of each country; therefore, few assumptions can be shared about how the introduction of policies affects the patient population. Mathematical models can contribute to further understanding these phenomena by obtaining analytical solutions as well as numerical simulations. In this study, an NPI was introduced into the SIR model for a conceptual study of infectious diseases under the condition that the transmission rate was reduced to a fixed value only once within a finite time duration, and its effect was analyzed numerically and theoretically. It was analytically shown that the maximum fraction of infected individuals and the final size could be larger if the intervention starts too early. The analytical results also suggested that more individuals may be infected at the peak of the second wave with a stronger intervention. This study provides quantitative relationship between the strength of a one-shot intervention and the reduction in the number of patients with no approximation. This suggests the importance of the strength and time of NPIs, although detailed studies are necessary for the implementation of NPIs in complicated real-world environments as the model used in this study is based on various simplifications.

Because of the global spread of COVID-19, our human society is facing a major public health crisis. The COVID-19 pandemic is caused by an emerging pathogen, SARS-CoV-2, for which there is no immunized population, causing an overshooting increase in the number of infected patients and depleting medical resources in many countries. Medical institutions are facing a difficult situation in which they must control second transmissions while treating critically ill patients, and as the number of patients increases, the medical system becomes swiftly tighter. When the number of patients exceeds the capacity, the quality of medical care deteriorates drastically, and the number of medical devices required for life support reaches its limit. This situation further increases the fatality rate of this infectious disease and causes serious damage to our society.

No effective treatment for COVID-19 has been established yet as of September, 2020, and only a public health approach can function as a control measure for the epidemic. To mitigate the spread of COVID-19, each country is implementing nonpharmaceutical interventions (NPIs) [1] to regulate social activities. NPIs comprise policies such as case isolation, voluntary home quarantine, closure of schools and universities, social distancing, stopping mass gatherings, and border closure. In major European countries, these NPIs were implemented, depending on the epidemic situation, during the first part of spring in 2020, with different timings and intensities [2, 3] . Although the first phase of the epidemic appeared to be suppressed by these mitigation measures, the re-epidemic became clearer in almost all countries because of deregulation after the first epidemic.

Because the control of epidemics by NPIs has caused a situation involving the imposition of strong restrictions on human socioeconomic activities, it would be desirable to study in advance the optimal timing, intensity, and duration of interventions that could bring about more promising results with minimal damage to society. Regarding COVID-19, for the purpose of ex-post verification, the influence of NPIs implemented in each country on the effective reproduction number was estimated. In European countries, the correlation between the decay of the effective reproduction number and the implementation of various NPIs has been verified [3] . The impact of travel limitations in China on the spread of infection has been discussed [4, 5] . However, a reliable estimation of the effects of NPIs is difficult, because the differences in NPI strategies employed by each country are strongly related to various background factors, such as the epidemic situations, social structure, legal systems, and culture [2] .

Mathematical modelling is an important method for estimating the effect of NPIs. In particular, the global pandemic of COVID-19 revealed that a situation in which only NPIs are effective against an emerging infectious disease is possible in societies in the 2020s. Such a situation had been "neglected" as a practical research target, which enhances the importance of theoretical approaches. Recently, compartmental models such as the SIR model [6] [7] [8] [9] have been extended to estimate the effects of NPIs on the number of patients [10] [11] [12] [13] . In other settings, optimal policies where intervention intensity can change continuously over time have been discussed in the context of minimizing objective functions [14] [15] [16] [17] . These solutions are very useful if we can estimate the effect of NPIs on the change in transmission parameters precisely.

The impacts of NPIs should be evaluated as a discontinuous change of the transmission rate in a model to represent the temporal discontinuity of intervention in the real world. The dynamics of the number of infected patients in continuous time in a system that includes discrete parameters, state spaces, and continuous-time dynamics can be modeled as a hybrid dynamical system [18] [19] [20] [21] [22] . This framework has been applied to mathematical models of infectious diseases [23] [24] [25] . In the simplest case, the one-shot intervention model, in which the intervention is implemented only once during the epidemic, can be used to discuss the theoretical dependence of the effects of NPIs on the timing and intensity [26] . As the COVID-19 epidemic continues, the accumulation of theoretical research on the effects brought about by NPIs has become even more significant. Recently, compartmental models with intervention have been studied both numerically [27, 28] and using some analytical methods [29, 30] in line with the COVID-19 epidemic.

In this study, we provide exact solutions of a simple SIR model with one-shot intervention, represented by a single discrete reduction in the transmission parameter during an epidemic. These solutions describe the dependence of the peak number of infected patients on the reproduction number under consideration of the implementation of NPIs and intervention timing. Theoretical and numerical analyses revealed non-trivial relations among the intensity of suppression of pandemics via NPIs, the number of infected individuals at the peaks, and the final size of infection cases. The methods and results shown in this study provide basic theoretical understanding in the context of the evaluation of NPIs.

In this study, we focus on the dependence of the maximum fraction of infected individuals on the timing of the NPIs. Note that this analysis is motivated by the COVID-19 epidemic, but we consider a hypothetical epidemic, whose properties are not necessarily the same as that of COVID-19.

We here introduce the SIR model, where the time evolution of the fraction of susceptible (s(t)), infected (i(t)), and removed (r(t)) individuals is given by the following ordinary differential equation (ODE):

taking into account the one-shot intervention and the second wave after the intervention. For simplicity, the total population is assumed to be unity, that is, s(t) + i(t) + r(t) = 1 holds. The summation of the right-hand sides of Eqs. (1)-(3) vanishes, which guarantees conservation of the total population. Therefore, as seen below, only Eqs. (1) and (2) are numerically integrated to obtain the time evolutions of s(t), i(t), and r(t). In this model, the basic reproduction number and the effective reproduction number at time t are given by R 0 = β/γ and R t = βs(t)/γ, respectively.

Kermack and McKendrick [6] derived the equation for the maximum fraction of infected individuals and showed that the final size is obtained by solving a transcendental equation with a given R 0 in the SIR model without any intervention [7] [8] [9] .

By applying the technique they employed in the derivation, we here provide the relationship between the fractions of susceptible and removed individuals at arbitrary times t 0 and t 1 , that is, s(t 0 ), r(t 0 ), and s(t 1 ), r(t 1 ) as

See Supplementary Section 1, Additional File 1 for details of the derivation. All analytical results presented in this paper are derived based on this equality.

In the present framework, an NPI in an isolated population is represented by a change in the transmission rate β. Let β = β off (> γ) be the transmission rate without the intervention, and it is switched to β on (< β off ) when the intervention starts at t = t on , and restored to β off at t = t off = t on + ∆t. Let the corresponding basic reproduction numbers be R 0,off = β off /γ and R 0,on = β on /γ, respectively. Here, ∆t denotes the duration of the intervention. In this paper, we study the effect of one-shot intervention, one of the simplest implementation schemes, where the intervention is implemented only once. The mathematical methods employed in this study can be applied to more complex cases, for example, where multiple interventions are implemented intermittently. If the fraction of susceptible individuals remains large enough and the herd immunity is not achieved after the intervention, a second wave occurs. It is thus necessary to consider this second wave to evaluate the effect of the intervention. The one-shot intervention setting has been numerically studied by Bootsma and Ferguson [26] , showing that the final size, which depends on the timing t on and t off , can be smaller for weaker intervention with larger R 0,on . Di Lauro et al. [28] numerically studied the dependence of the final size and the peak fraction of infected individuals on the timing of the intervention, where the intervention was assumed to start when the fraction of infected and recovered individuals exceeds a threshold value. In particular, they concluded that the onset timing should be chosen so that the two peaks during and after the intervention are comparable. Morris et al. [29] showed that the maximum fraction of infected individuals with one-shot interventions can approach that achieved by the optimal intervention, which requires an unrealistic intervention, such as R 0,on = 0. Sadeghi et al. [30] also suggested the existence of the optimal timing of the intervention based on discussion using a linearized equation and numerical simulation.

Equations (1) and (2) were numerically integrated using ODEPACK [31] , which consists of nine Fortran solvers for the initial value problem for ODE systems. The time step was set to 0.001. The nonintervention transmission rate and the recovery rate are fixed as β off = 2/7 days −1 and γ = 1/7 day −1 , that is, R 0,off = 2. The initial condition s(0) = 1 − ǫ, i(0) = ǫ, and r(0) = 0, where ǫ = 0.001 is employed in this study. As discussed in detail below, the system behaves qualitatively differently for different values of R 0,on . Here, we primarily focus on the cases for R 0,on = 1.4 ( Fig. 1 ) and R 0,on = 0.7 (Fig. 2 ), corresponding to relatively weaker and stronger intervention intensities, respectively. Note that R 0,on > 1 in the former case implies that the infection may spread even in the presence of the intervention. The time series without the intervention is shown in Figs. 1(A) and 2(A). They are identical, because R 0,on does not affect the dynamics without the intervention. Figures  1(B) -(E) and 2(B)-(E) show the time series with an intervention with a constant intervention duration ∆t = 60 days. The second wave can occur if the herd immunity is not achieved when the intervention ends (Figs. 1(C) and 2(B), (C)). The intervention is implemented for ton ≤ t ≤ t off = ton + ∆t and are represented by grey intervals in these panels. Conditions for the peaks of the fraction of infected individuals are given by s(t) = 1/R 0,off , when the intervention is not implemented, and by s(t) = 1/R 0,on during intervention. The final size of the outbreak for each case is represented by r(∞). Panels (F) and (G) represent the maximum fraction of infected individuals imax normalized by that without the intervention i 0 max , plotted in terms of r(ton) and ∆r = r(t off ) − r(ton), and ton and ∆t, respectively. Symbols (A)-(E) in panels (F) and (G) denote the intervention timings of the time series of the corresponding panels in (A)-(E). Symbols (i)-(iv) in panel (F) denote the timing that the maximum infected fraction is observed, as described in the main text. The boundaries between regions shown in panel (F) are obtained analytically. See the main text for the details.

In addition to the numerically obtained time series, we analytically show that peaks of the fraction of infected individuals can appear at the following four timings by applying Eq. (4): after the intervention, during the intervention, at the onset of the intervention, and before the intervention. Note that there can be two peaks of the fraction of infected individuals if a second wave occurs. In such a case, the peak with a larger fraction gives the global maximum. Let us describe the four cases with respect to the timing of the maximum in detail. • (i) The maximum appears after the intervention (Figs. 1(B), (C), 2(B), and (C)). If the intervention ends before achieving herd immunity, a peak is observed during the second wave after the intervention. The fraction of infected individuals for this peak gives the global maximum if it is higher than the first peak before or during the intervention.

• (ii) The maximum appears during the intervention ( Fig. 1(D) ). If the effective reproduction number R t declines and crosses unity during the intervention, there is a peak in this timing. The condition R 0,on > 1 is necessary for the existence of this peak, because the effective reproduction number has to be larger than unity at the onset of the intervention. Therefore, this peak does not appear for R 0,on < 1. This peak is the global maximum if it is larger than the second peak. The peak also appears at this timing in Fig. 1 (C), but the second peak is slightly higher than this peak.

• (iii) The maximum appears at the onset of the intervention at t = t on ( Fig. 2(D) ). If R ton < 1 holds when the intervention starts, the fraction of infected individuals decreases during the intervention, and we observe a local peak at t = t on . There may be another peak if herd immunity is not achieved during the intervention. If the fraction of infected individuals at the second peak is less than at this peak, the timing of the global maximum is given at t on . Note that this region also appears for R 0,on > 1, although the time series is not shown in Figs. 1(B)-(E).

• (iv) The maximum appears before the onset of the intervention (Figs. 1(E) and 2(E)). The fraction of infected individuals reaches its maximum before the intervention. This case implies that the intervention starts too late and fails to mitigate outbreaks in terms of the maximum fraction of infected individuals. The fraction of infected individuals at the peaks can be calculated analytically. Conditions for the peaks are given in terms of the fraction of susceptible individuals as s(t) = 1/R 0,off = γ/β off without the intervention and s(t) = 1/R 0,on = γ/β on with the intervention. The maximum fraction of infected individuals for cases (ii), (iii), and (iv) does not depend on the fraction of removed individuals at the offset of the intervention r(t off ), because these peaks appear before t off . For case (i), the maximum fraction of infected individuals i max is given as

which depends on ∆r := r(t off ) − r(t on ), the difference in the fraction of removed individuals between the onset and offset of the intervention. See Supplementary Section 2, Additional File 1 for the explicit form of i max for all cases and its derivation. The boundaries between regions (i) and (ii), (i) and (iii), (ii) and (iii), and (iii) and (iv) can be obtained analytically with respect to r(t on ) and ∆r (Figs. 1(F) and 2(F)). On the boundaries, the fraction of infected individuals at two peaks is comparable (Figs. 1(C) and 2(C)). See Supplementary Section 3, Additional File 1 for details of the derivation.

The final size also reflects the effect of the intervention. The final size of removed individuals r(∞) with the intervention is obtained by solving the equation

in a self-consistent manner [32] . Specifically, as r(∞) appears on both sides, this equation can be solved numerically or using the Lambert W function, except for some special cases. This equation implies that the final size depends on ∆r = r(t off ) − r(t on ). Another important implication of this equation is that r(t off ) has the upper boundr depending on r(t on ), which is given bỹ

See Supplementary Section 4, Additional File 1 for the details of the derivation. Using this equality, one can show that the final size in the presence of the intervention is always smaller than that without the intervention. For

one can achieve r(∞) ≈ 1 − 1/R 0,off by setting t on properly, which is the smallest prevalence to achieve herd immunity, with an intervention duration ∆t that is large enough [32] . See Supplementary Section 5, Additional File 1 for details of the derivation. Numerical results regarding the final size are summarized in Supplementary Section 6, Additional File 1.

We report the numerical and analytical results, when the reproduction number under the intervention R 0,on is large ( Fig. 1 ) and small (Fig. 2) , showing qualitatively different behaviors. In Figs. 1(F), (G), 2(F), and (G), the dependence of the maximum fraction of infected individuals on the timing of the intervention is plotted. Here, i max is normalized by that in the absence of the intervention i 0 max (Figs. 1(A) and 2(A)). As the maximum infected fraction drops in the presence of the intervention, i max /i 0 max is less than unity and quantifies the effectiveness of the intervention in terms of the maximum fraction of infected individuals. As this ratio decreases, the intervention shows more success in reducing the maximum fraction.

It is difficult to obtain the time series of the SIR model analytically without any approximations, for example, linearization, or the method presented in [6] . Therefore, i max is numerically computed with different onset and offset times for the intervention, t on and t off (Figs. 1(G) and 2(G)). However, i max can be analytically calculated with respect to the fraction of the recovered individuals at the onset, r(t on ), and the offset, r(t off ), of the intervention (Figs. 1(F) and 2(F)). Note that there exists a one-to-one correspondence between {t on , t off } (panels (F)) and {r(t on ), r(t off )} (panels (G)) in Figs. 1 and 2.

Weak intervention (large R 0,on = 1.4) In this case, the reproduction number is larger than unity even in the presence of the intervention. Therefore, the fraction of infected individuals may increase during the intervention period. The peaks of infected individuals can be observed during the intervention, and the maximum infected fraction can appear at any of the four timings (i)-(iv) classified above. Figures 1(B) -(E) show that i max is minimized in the intermediate onset time t on near (C), where the peaks during and after the intervention are comparable. This non-monotonic dependence on t on is clearly visualized in Fig. 1(G) . The peak of infected individuals during the intervention is smaller than that without intervention. This intervention mitigates the second wave.

As shown in Eq. (5) and Fig. 1(F) , the maximum fraction of infected individuals is linear in ∆r if the peak of the second wave is the maximum, that is, case (i). This is verified by the fact that the contours lie horizontally in case (i). To clarify this point, the contours are explicitly shown in Supplementary Section 7, Additional File 1. If the fraction of the infected individuals reaches the maximum during or at the onset of the intervention (cases (ii) and (iii), respectively), the maximum fraction depends only on t on , which has one-to-one correspondence to r(t on ). This is verified in Fig. 1(F) , (G), and Supplementary Section 7, Additional File 1. If the maximum appears before the intervention, that is, case (iv), the maximum is independent of both t on and t off .

Onset and offset times for the intervention with constant ∆t corresponding to Figs. 1(A)-(E) are plotted in Figs. 1(F) and (G). As suggested by the time series, the maximum infected fraction is smallest in the intermediate intervention onset t on , near (C). It is clear from panel (F) that this point is located close to the boundary between cases (i) and (ii), where peaks during and after the intervention are comparable. In case (ii), the maximum does not depend on ∆r, which implies that a longer intervention does not reduce the maximum. Note that the relationship between (t on , t off ) and (r(t on ), r(t off )) is non-monotonous.

Strong intervention (small R 0,on = 0.7) When the transmission rate is small during the intervention, the maximum infected fraction is minimized in the intermediate starting time of the intervention t on , namely, early implementation of the intervention does not necessarily minimize the infected fraction (Figs. 2(F), (G)). This result is intuitively understood as follows. If the intervention starts too early, the infection does not spread because of the small intervention transmission rate. Therefore, the second wave after the intervention is large, thus the early intervention is not effective. If the timing of the intervention is characterized in terms of r(t on ) and r(t off ) ( Fig. 2(F) ), the maximum of the second wave depends only on ∆r. If the maximum fraction is found during the intervention, its value is independent of t off and r(t off ). The peak does not appear during the intervention, that is, case (ii) does not appear because R t,on < 1 holds.

Maximum infected fractionī max versus reproduction number under the intervention R 0,on For each R 0,on , there exist onset and offset timings of the intervention that minimize the maximum fraction of infected individuals i max (Figs. 1(F), (G) and 2(F), (G)). Let this value beī max (R 0,on ). Figure 3 plots the dependence ofī max on R 0,on . As seen in the figure, it is minimized at a non-trivial intermediate value of R 0,on = R * 0,on ≈ 1.23 > 1. This implies that the maximum fraction of infected individuals is minimized for a weak intervention under the one-shot condition. For R 0,on ≥ R * 0,on , i max is achieved at the boundary between regions (i) and (ii) at t on = 0, where the peaks during and after the intervention are comparable ( Fig. 1(C) ). For R * 0,on ≥ R 0,on ≥ R (1) 0,on , the boundary between regions (i) and (ii) with t on > 0 givesī max , where R (1) 0,on ≈ 1.08 is the parameter value below which region (ii) does not exist. For strong intervention R 0,on ≤ R (1) 0,on ,ī max is found at the boundary between regions (i) and (iii), where the peaks of the onset and after the intervention are comparable ( Fig. 2(C) ), at t on > 0. Conditions forī max are analytically derived in all cases, shown by the solid line in Fig. 3 . The conditions forī max for R 0,on ≥ R 

It is possile to optimize the onset and offset timings of the NPI by minimizing an objective function under certain constraints. In the present framework, this can be formulated as an optimization problem in a hybrid nonlinear dynamical system. The optimal intervention strategy depends on the objective function. Here, we discuss the following two simple scenarios to minimize i max . In general, more complex objective functions can be used for optimization. The following discussion may provide some intuition for considering such cases, taking into account the second wave.

Minimizing i max with a constraint in the intervention duration

Let us start with a case where the intervention duration ∆t is less than a certain value. As i max is monotonically decreasing in ∆t for a fixed t on , we can assume that ∆t is a constant, and t on is varied. This case has been studied in Figs. 1(B) -(E) and 2(B)-(E). This scheme fixes the intervention duration, so it is easier to anticipate the economic impact of the intervention, which depends on the duration of the intervention, than in the next scenario based on the fraction of recovered individuals.

As discussed above, i max is minimized in the intermediate onset time t on , when the first peak during the intervention and the second peak after the intervention are comparable. For larger R 0,on , t on giving the minimum i max converges to zero for large ∆t (Fig. 1(G) ), but converges to non-zero for large ∆t for smaller R 0,on (Fig. 2(G) ). These results are understood as follows: for a strong intervention with small R 0,on , the intervention immediately suppresses the fraction of infected individuals. Therefore, the early onset of the intervention prevents more individuals from achieving immunity during the intervention and eventually increases the maximum fraction of infected individuals.

Another possible constraint is to minimize the intervention duration ∆t, keeping the maximum fraction of infected individuals i max constant. This corresponds to choosing t on along a contour of i max , which prevents the overcapacity of medical support. In Figs. 1(G) and 2(G), ∆t is minimized for an intermediate t on . Namely, a contour that crosses point (C) (∆t = 60) reaches ∆t ≈ 20 if the onset time is later than t on = 19.2 ( Fig. 1(C) ). As we have already discussed, the maximum fraction depends on ∆r in region (i). It takes a shorter time to achieve the same ∆r if the intervention starts later in this case. As the maximum fraction of infected individuals depends on r(t on ) but is independent of ∆r in cases (ii) and (iii), the optimal r(t on ) is determined by the boundary between regions (i) and (ii) or that between (i) and (iii) in Figs. 1(F) and 2(F). Evidently this timing is an intermediate value of t on > 0. It may be easier to set a plan for the intervention with respect to r(t on ) and r(t off ), rather than the timing t on and t off .

Similar optimization problems can be considered using the final size as the objective functions. These problems are discussed in Supplementary Section 9, Additional File 1.

In the present study, in consideration of the actual COVID-19 situation, we studied the situation in which the reproduction number of an infectious disease is temporarily reduced by implementing an NPI once during the epidemic, using a simple mathematical model. The results provide theoretical implications as to how strong NPIs should be introduced during an epidemic of an emerging infectious disease. If the effective reproduction number during the intervention is too small, the fraction of infected individuals at the peak in the second wave may be higher than the first peak. It was also shown numerically and analytically that the fraction of infected individuals can also increase if the intervention is started too early. The upper limit of medical capacity is an essential practical constraint. In particular, for infectious diseases such as COVID-19, which is too emergent to expect effective treatments, it is more important to avoid exhausting the medical system. This study suggests that it will be necessary to be alert for a larger second wave that may occur after strong intervention in such cases.

We analytically derived the peak fractions of infected individuals in the SIR model with the one-shot intervention. These analytical results suggest that the peak fractions can be smaller with non-trivial intervention timings. The maximum fraction is smallest for R 0,on > 1, that is, the intervention reproduction number is not in the disease-free regime. In the literature, Bootsma and Ferguson [26] showed that the final size can be minimized for non-trivial intervention timing and R 0,on numerically. Di Lauro et al. [28] numerically showed that the peak fraction of infected individuals also depends on the timing of the intervention. We obtained analytical expressions for these quantities in this study. The analytical results for the final size are presented in Supplementary Section 6, Additional File 1. It should be noted that the formulae for the peak fraction depend on the timing of the peak, resulting in various cases compared with the final size. Some analytical results are available for this system. Sadeghi et al. [30] explained these non-trivial effects based on solutions of the linearized equation, which exponentially grows and decays without and with the intervention, respectively. Linearization is one of the simplest approximations and is applicable in this case; in particular, this approximation is useful in discussions regarding timing. Conversely, linearization cannot be used to discuss the important case where R 0,on > 1, as the linearized equation cannot explain the declining number of infected individuals. The proposed method in this study provides a unified framework, including the cases where linearization is not feasible. Morris et al. derived [29] an equation for the peaks of the fraction of infected individuals in terms of s(t on ), s(t off ), i(t on ), and i(t off ). In this work, we further show the peak fractions in terms of the two parameters r(t on ) and r(t off ) for a general initial condition.

Although concrete measures and guidelines for COVID-19 are required, it is emphasized again that the results in this study were derived using a simplified model with many assumptions. First, the results presented in this study are based on the simplified SIR model, which takes into account neither the realistic pathology of COVID-19 nor societal response. This model assumes uniform and random contact within a group and does not consider interactions between different subgroups in the population. Recovered patients are assumed to have complete immunity in this model. These assumptions are not applicable to COVID-19, where there are still many unclear factors regarding heterogeneity in contact networks and the immune response of patients after recovery. If a society develops other public health measures during the intervention, the basic reproduction numbers before and after the NPI may be different. Next, this study is a one-shot intervention model, in which the transmission coefficient returns to the original value after a single intervention. Practically, each re-pandemic requires multiple intermittent interventions [33] , making the intervention process much more complex [34] . Furthermore, the model is based on a deterministic dynamical system of an isolated population. Therefore, important NPI measures such as border control cannot be estimated in the present framework. The deterministic nature of the model assumes that the infectious disease cannot be eradicated, as the number of infected individuals remains non-zero for a finite time. If the population is small enough and no imported cases are assumed, strong and early intervention, which is not necessarily recommended in this study, may eradicate the disease, and a second wave does not occur. Strong intervention would be necessary in other cases, for example, when the number of infected patients approaches the capacity of the medical system. These complex situations may be analyzed in detail by extending the methods presented in this study, which may lead to different conclusions from those reached in this study.

In this study, the effect of NPIs was modeled as a hybrid dynamical system, which may further enable us to approach more refined models in future investigations. The influence of NPIs with respect to political decisions and behavioral changes of people can be expressed more accurately by introducing a hybrid dynamical system. In recent years, dynamical systems theory and control theory have been developed, and phenomena specific to hybrid systems such as Zeno solutions and sliding motions have been discussed [18] [19] [20] [21] [22] . Some studies, such as [35] and [36] , have proposed control of infectious diseases using sliding mode control. The optimal policy under NPIs can be discussed by modeling the effect of economic damages associated with execution and then minimizing the cost function. Discrete changes in parameters such as the transmission coefficient in NPI implementation are strongly linked to the intensity of measures, suppression of economic activity, and changes in human mobility. For example, the correlation between the decrease in human mobility with NPIs and the effective reproduction number for the COVID-19 pandemic has been studied [37, 38] . Such studies can contribute to modeling the costs of NPIs. The mathematical model of the epidemic suppression effect can be constructed using a hybrid dynamical system, taking into account the negative socioeconomic impact of NPIs.

The authors declare that they have no competing interests. Supplementary materials: Analytical estimation of maximum fraction of infected individuals with one-shot non-pharmaceutical intervention in a hybrid epidemic model S1 Method: Derivation of Eq. (4) [6] [7] [8] [9] We want to derive Eq. (4) in the main text from the evolution equation of the SIR model

We integrate Eq. (1) from time t 0 to t 1 :

The integral on the right-hand side can be replaced with the integration of Eq. (3),

Then, we obtain

(S16) which is Eq. (4) in the main text. We apply this equation to various t 0 and t 1 .

As described in the main text, there are four possible timings where the maximum fraction of infected individuals i max is observed. We summarize the analytical forms of the maximum fraction in all cases with respect to the fractions of removed individuals at the onset r(t on ) and the offset r(t off ) of the intervention, that is, we assume that the timing of the onset and offset of the intervention is determined by the fraction of removed individuals.

• (i) The maximum appears after the intervention:

where ∆r = r(t off ) − r(t on ).

• (ii) The maximum appears during the intervention:

(S18)

• (iii) The maximum appears at the onset of the intervention:

(S19)

• (iv) The maximum appears before the onset of the intervention:

Here, we derive Eqs. (S17)-(S20). The conditions for the fraction of infected individuals to be the local peak are given by s(t) = γ/β off = 1/R 0,off without the intervention and s(t) = γ/β on = 1/R 0,on with the intervention. First, let us apply Eq. (4) until the onset time of the intervention, with β = β off , t 1 = t on , and t 0 = 0. Using the initial condition s(0) = 1 − ǫ, i(0) = ǫ, r(0) = 0, we obtain

In the following, we take the limit ǫ → 0. Substituting Eq. (S21) into Eq. (4), one can derive the fraction of susceptible and removed individuals during and after the intervention.

We derive the fraction of infected individuals after the intervention. Let t after p > t off be the time at which the infected fraction reaches the peak after the intervention. The relationship between s(t off ) and r(t off ) is derived by substituting Eq. (S21) with ǫ → 0 into Eq. (4) with t 0 = t on and t 1 = t off as follows:

Then, by applying Eq. (4) setting t 0 = t off and t 1 = t after p , we have

The condition for the peak after the intervention is s(t after p ) = γ/β off . Substituting this condition into Eq. (S23), the fractions of removed and infected individuals are given by:

which is equivalent to Eq. (S17) with R 0,off = β off /γ and R 0,on = β on /γ. Equation (S26) gives i max if the maximum appears after the intervention. Note that the fraction of infected individuals depends only on the onset and offset timings as r(t off ) − r(t on ). The later the onset or sooner the offset is, the higher the peak is.

(ii) Maximum fraction during the intervention If the effective reproduction number at the onset of the intervention is greater than unity, β on s(t on )/γ > 1, and that at the offset is less than unity, β on s(t off )/γ < 1, then there is a peak during the intervention. Let t during p be the time at which the infected fraction reaches the peak, where t on < t during p < t off . We obtain the relationship between susceptible and removed individuals during the intervention by applying Eq. (4) with t 0 = t on and t 1 = t during p and using Eq. (S21) as the relationship between s(t on ) and r(t on ):

The condition for the peak during the intervention is s(t during p ) = γ/β on . By substituting this condition into Eq. (S27), we obtain

Then, the fraction of infected individuals at this peak is

which is equivalent to Eq. (S18). As β off > β on and β off −βon βon > 0, the fraction of infected individuals at this peak linearly increases compared with the fraction of removed individuals at the onset of the intervention r(t on ). Equation (S30) is the condition for the maximum if this peak is higher than another peak.

The conditions for the existence of this peak are β on s(t on ) > γ and β on s(t off ) < γ, which are interpreted as

Therefore, β on > γ is required for the existence of this peak.

(iii) Maximum fraction at the onset of the intervention At the onset of the intervention, the fraction of infected individuals is given as

by substituting Eq. (S21) . This peak appears if the effective reproduction number before the intervention is larger than unity β off s(t on ) > γ and that after the onset of the intervention is less than unity β on s(t on ) < γ. These conditions are summarized in terms of r(t on ) as follows:

by substituting the above conditions into Eq. (S21).

(iv) Maximum fraction before the intervention Let t 1 = t before p be the time at which the peak appears before the onset of the intervention. We apply Eq. (4) with β = β off , t 0 = 0, and t 1 = t before p . Then, we obtain

.

The peak condition s(t before p ) = γ/β off leads to

Substituting this equation into the conservation of total population i(t) = 1 − s(t) − r(t), we obtain

This peak appears for r(t on ) > γ β off log β off γ , meaning a late onset of the intervention.

As shown in Figs. 1(F) and 2(F) in the main text, the timing giving the maximum fraction of infected individuals switches in the (r(t on ), ∆r) plane, where ∆r = r(t off ) − r(t on ). These figures suggest possible transitions between regions (i) and (ii), between (i) and (iii), between (ii) and (iii), and between (iii) and (iv) listed above. At the boundaries, the peaks of the different timings are expected to be equal. Based on this idea, we sketch the derivations of the equations for the boundaries between these timings.

In regions where peaks of infected individuals during and after the intervention coexist, the global maximum is given by

There occurs a transition of the timing giving the global maximum, and the condition at the boundary is i(t after p ) = i(t during p ), that is,

This condition gives the boundary between regions (i) and (ii) in Fig. 1(F) in the main text. For ∆r smaller than this condition, the peak of the second wave is larger than that during the intervention.

Two peaks at the onset of the intervention and after the intervention can coexist. The maximum fraction of infected individuals is

The condition for the boundary i(t on ) = i(t after p ) can be rewritten as

This boundary is shown in Figs. 1(F) and 2(F) in the main text.

Peaks at the onset and during the intervention cannot coexist. As shown in Eqs. (S31) and (S35),

gives the boundary between regions (ii) and (iii). This condition does not depend on ∆r (Fig. 1(F) in the main text).

As shown in Eqs. (S35) and (S37),

gives the condition for the boundary between regions (iii) and (iv). As in the case of the boundary between (ii) and (iii), this condition does not depend on ∆r (Figs. 1(F) and 2(F) in the main text).

As discussed in Eq. (S23), the relationship between s(t 1 ) and r(t 1 ) after the intervention at t 1 > t off is

We obtain the final size equation for t → ∞ with the intervention by substituting this into the condition s(∞) + r(∞) = 1 as

which is Eq. (6) in the main text.

Let us study the features of Eq. (6) . Note that we assume that R 0,on is fixed in this section. First, we show that the final size monotonically decreases with respect to r(t off ) for a fixed r(t on ). Regareing Eq. (6) as an implicit function of r(∞) and ∆r, we obtain

In the final state, the inequality 0 ≤ 1 − r(∞) ≤ γ/β off must hold because of the herd immunity condition. Therefore,

holds. As ∆r monotonically increases with respect to t off for a fixed t on , the final size is smaller if the intervention is longer. Equivalently, r(∞) monotonically increases with respect to r(t off ) for a fixed r(t on ).

Note that there exists the upper bound in r(t off ). Let this upper bound bẽ r(r(t on )). It satisfies Eq. (7) in the main text:

which is obtained forr = r(∞) in Eq. (6) . This upper bound is observed in Figs. 1(F) and 2(F) in the main text.

Finally, let us discuss r(t on ) and r(t off ) which minimizes the final size with fixed R 0,on . The monotonicity of r(∞) with respect to r(t off ) implies that the lower bound of the final size r(∞) with fixed r(t on ) is obtained if r(t off ) is replaced withr(r(t off )). Therefore, we focus on the final size equation in which r(t off ) is replaced withr,

To minimize r(∞) in this equation, we should maximizer(r(t on )) − r(t on ) with respect to r(t on ). This condition is given by

This can be satisfied for

By substituting these equations into Eq. (S54), we obtain the final size equation: For sufficiently large ∆t and small enough R 0,on , the final size reaches close to 1 − 1 R 0,off = 0.5, which is the lower bound of the final size necessary to achieve herd immunity, for r(t on ) ≈ 0.295 (Fig. 2(D) in the main text), which corresponds to Weak intervention (large R 0,on = 1.4) The final size can be minimized with t on = 0, the early implementation of the intervention, in this case (Fig. S4 ). The timing of the intervention that minimizes the final size is different from that which minimizes the maximum fraction of infected individuals. In general, the optimal timing of the intervention depends on the objective function to be minimized. The theoretical prediction in Eq. (6) (Eq. (6)), that the final size depends on ∆r, is numerically verified (Fig. S4(A) ). For a constant intervention duration ∆t (Fig. S4(B) ), the final size is minimized with intermediate t on (Fig. 1(D) in the main text).

Strong intervention (small R 0,on = 0.7) The simulation results (Fig. 2(B) -(E) in the main text) suggest that the final size is minimized with an intermediate starting time for small R 0,on . Indeed, as shown in Fig. S5 , the final size can be minimized at t on = 0. The minimized final size is closer to r(∞) = 1 − 1 R 0,off = 0.5, which is the minimum fraction of removed individuals necessary to achieve herd immunity, compared with the case of R 0,on = 1.4 (Fig. 2(D) in the main text). appears during or at the onset of the intervention (cases (ii) and (iii)), i max depends on r(t on ) and does not depend on ∆r. Therefore, the contours in the (r(t on ), ∆r) plane are the horizontal and vertical lines for the former and latter cases, respectively. To visualize them explicitly, Figs. S6(A) and S7(A) plots the five contours for i max . As a theoretical prediction, they consist of horizontal and vertical lines. Equation (6) theoretically predicts that the final size r(∞) with the intervention depends on ∆r only, which implies that the contours are horizontal lines (Figs. S6(B) and S7(B)).

S8 Appendix: Dependence ofī max on R 0,on We want to minimize the maximum fraction of infected individualsī max for a fixed R 0,on by varying t on and t off . Note that i max is decreasing function in ∆r (Eq. (S19)) for case (i), and i max does not depend on ∆r in regions (ii) and (iii). Therefore, for a fixed r(t on ),ī max is achieved either in region (ii) or (iii). Here,ī max can be analytically or semi-analytically calculated by evaluating i max on the boundaries between (i) and (ii) or between (i) and (iii). We show below that the dependence of i max on R 0,on changes, in which the boundary givesī max , leading to its non-trivial dependence on R 0,on .

Weak intervention (R 0,on ≥ R * 0,on ≈ 1.23) For large R 0,on ,ī max is achieved at the boundary between regions (i) and (ii), where the peaks of the infected individuals during and after the intervention are equal ( Fig. 1(F) in the main text) . To discuss the parameter regions where this boundary exists, it is necessary to take into account the upper bound of removed individuals r(r(t on )). If r(t off ) = r(t on ) + ∆r at the boundary between regions (i) and (ii) given by Eq. (S41) is larger thanr, this boundary is not observed at t on . In particular, it should be verified whetherr(0) is less than ∆r in Eq. (S41) for r(t on ) = 0. Let R * 0,on be the value of R 0,on , below which the boundary between (i) and (ii) does not exist. For R 0,on ≤ R * 0,on , the boundary crossesr at a non-zero intervention onset time r(t on ), because the boundary is decreasing in r(t on ) (Eq. (S41)) andr is increasing in r(t on ) for small r(t on ).

For R 0,on ≥ R * 0,on , where there is a boundary between regions (i) and (ii) at t on = 0 ( Fig. 1(F) in the main text), ∆r can be the largest at this boundary at t on = 0. Equation (S17) suggests that a larger ∆r yields a smaller i max , andī max is given by substituting r(t on ) = 0 into Eq. (S30) as follows:

This equation is used as the theoretical curve for R 0,on ≥ R * 0,on in Fig. 3 in the main text. This equation implies thatī max is increasing in R 0,on for R 0,on ≥ R * 0,on . See the next subsection for the value of R * 0,on .

Medium-intervention (R * 0,on ≥ R 0,on ≥ R

(1) 0,on ≈ 1.08) If R 0,on is smaller than R * 0,on , region (ii) does not give the global maximum of the fraction of infected individuals for t on = 0. As a result, a larger ∆r, which indicates a smaller i max , is realized for t on = 0. The upper bound of r(t off ) at the boundary between regions (i) and (ii) equalsr(r(t on )). Substituting the condition for the boundary between (i) and (ii) (Eq. (S41)) in the limit r(t off ) →r into the condition forr (Eq. (7)), we obtaiñ

where A := 1 + 1 R 0,off − R 0,on R 0,off log (R 0,on ) − R 0,on log (R 0,off ) .

Substituting this into Eq. (S41) again, r(t on ) at the boundary can be explicitly given as

Substituting these into Eq. (S17) or (S18), we obtain

which is the theoretical curve for R * 0,on ≥ R 0,on ≥ R

0,on in Fig. 3 in the main text. Note thatī max is a decreasing function in R 0,on . Therefore,ī max is minimized at R 0,on = R * 0,on . The condition for R * 0,on is obtained by equating Eqs. (S60) and (S64). It is difficult to calculate R * 0,on explicitly, and it was numerically obtained as R * 0,on ≈ 1.23 in the present parameter setting.

Strong intervention (R 0,on ≤ R

(1) 0,on ) For small R 0,on ,ī max is determined differently. It is understood that a peak in region (ii) does not give i max for R 0,on ≤ R (1) 0,on , where R (1) 0,on ≈ 1.08. In this case, the boundary between regions (i) and (iii) satisfying r(t off ) =r(r(t on )) gives the condition forī max (Fig. 2(F) in the main text) . This condition cannot be solved explicitly with respect toī max . Therefore, we give the parametric equations forī max and R 0,on in terms of r(t on ). The parametric equation forī max is i max = 1 − e − β off γ r(ton) − r(t on ).

(S65)

Next, we derive the parametric equation for R 0,on . To this end, let us rewrite (S44) as ∆r = B(r(t on )) R 0,off − R 0,on ,

where B(r(t on )) = R 0,off [e − β off γ r(ton) + r(t on )] − 1 + log (R 0,off ) .

Substituting this into Eq. (7), we obtain a transcendental equation forr:

which can be solved as

where W −1 (·) denotes the branch of the Lambert W function which gives the real value other than the principal branch. Finally, we obtain the parametric equation for R 0,on as R 0,on = R 0,off − B(r(t on )) 1 + 1 R 0,off W −1 (−R 0,off exp [B(r(t on )) − R 0,off ]) − r(t on )

.

Equations (S65) and (S70) are used to plot the theoretical curves for R 0,on ≤ R

(1) 0,on in Fig. 3 in the main text.

S9 Appendix: Optimization problems with respect to final size Two scenarios for the optimization problems in the hybrid nonlinear dynamical systems are presented in the main text. Here, we discuss the problems using the final size as the objective functions.

Targeting the final size to be minimized, we discuss two constraints based on Figs. S4 and S5. As in the case for the minimization problem of the maximum fraction of infected individuals, the final size is minimized at an intermediate t on with the constraint of a constant ∆t (Figs. 1(D) and 2(D) in the main text, and Figs. S4 (right) and S5 (right)). Note that this t on is different from the one which minimizes i max . This result is intuitively understood in the right panels of these figures. The contour of the final size depends on ∆r. For a constant intervention duration, ∆r is largest with intermediate r(t on ), which maximizes the final size with this constraint.

Minimizing intervention duration with a constraint on the final size It is possible to consider this rule to minimize the intervention duration along the contour of the final size in Figs. S4 and S5. This corresponds to intervening while keeping ∆r constant. There exists an intermediate t on that minimizes the final size in this setting as well. 

Non-pharmaceutical public health interventions for pandemic influenza: an evaluation of the evidence base

Countries test tactics in 'war' against covid-19

Estimating the number of infections and the impact of non-pharmaceutical interventions on covid-19 in 11 european countries

The effect of travel restrictions on the spread of the 2019 novel coronavirus

The effect of human mobility and control measures on the covid-19 epidemic in china

A contribution to the mathematical theory of epidemics

Infectious Diseases of Humans: Dynamics and Control

The basic reproduction number of infectious diseases: computation and estimation using compartmental epidemic models

Age-structured Population Dynamics in Demography and Epidemiology

Modelling the influence of human behaviour on the spread of infectious diseases: a review

Determinants of the spatiotemporal dynamics of the 2009 h1n1 pandemic in europe: Implications for real-time modelling

Inferring the causes of the three waves of the 1918 influenza pandemic in england and wales

Appropriate models for the management of infectious diseases

An optimal control theory approach to non-pharmaceutical interventions

Merging economics and epidemiology to improve the prediction and management of infectious disease

Optimal vaccination policies for an sir model with limited resources

Optimal control of an sir epidemic through finite-time non-pharmaceutical intervention

Differential Equations with Discontinuous Righthand Sides: Control Systems. Kluwer

Handbook of Hybrid Systems Control: Theory, Tools, Applications

Theory of hybrid dynamical systems and its applications to biological and medical systems

Discontinuity-induced bifurcations of piecewise smooth dynamical systems

Stability analysis and controller synthesis for hybrid dynamical systems

Pulse mass measles vaccination across age cohorts

Global dynamics of a piece-wise epidemic model with switching vaccination strategy

Global dynamics of sir model with switched transmission rate

The effect of public health measures on the 1918 influenza pandemic in us cities

How will country-based mitigation measures influence the course of the covid-19 epidemic?

The timing of one-shot interventions for epidemic control

Optimal, near-optimal, and robust epidemic control

Universal features of epidemic models under social distancing guidelines

Odepack, a systematized collection of ode solvers. Scientific computing

Effects of infection control policies on the final size

Projecting the transmission dynamics of sars-cov-2 through the postpandemic period

Sliding mode control of outbreaks of emerging infectious diseases

Analysis and control of seir epedemic model via sliding mode control

Association between mobility patterns and covid-19 transmission in the usa: a mathematical modelling study. The Lancet Infectious Diseases

Non-compulsory measures sufficiently reduced human mobility in japan during the covid-19 epidemic