Epidemics Forecast From SIR-Modeling, Verification and Calculated Effects of Lockdown and Lifting of Interventions

Due to the current COVID-19 epidemic plague hitting the worldwide population it is of utmost medical, economical and societal interest to gain reliable predictions on the temporal evolution of the spreading of the infectious diseases in human populations. Of particular interest are the daily rates and cumulative number of new infections, as they are monitored in infected societies, and the influence of non-pharmaceutical interventions due to different lockdown measures as well as their subsequent lifting on these infections. Estimating quantitatively the influence of a later lifting of the interventions on the resulting increase in the case numbers is important to discriminate this increase from the onset of a second wave. The recently discovered new analytical solutions of Susceptible-Infectious-Recovered (SIR) model allow for such forecast. In particular, it is possible to test lockdown and lifting interventions because the new solutions hold for arbitrary time dependence of the infection rate. Here we present simple analytical approximations for the rate and cumulative number of new infections.


INTRODUCTION
The Susceptible-Infectious-Recovered (SIR) model has been developed nearly hundred years ago [1,2] to understand the time evolution of infectious diseases in human populations. The SIR system is the simplest and most fundamental of the compartmental models and its variations [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. The considered population of N ≫ 1 persons is assigned to the three compartments s (susceptible), i (infectious), or r (recovered/removed). Persons from the population may progress with time between these compartments with given infection (a(t)) and recovery rates (μ(t)) which in general vary with time due to non-pharmaceutical interventions taken during the pandemic evolution.
Let I(t) i(t)/N, S(t) s(t)/N and R(t) r(t)/N denote the infected, susceptible and recovered/removed fractions of persons involved in the infection at time t, with the sum requirement I(t) + S(t) + R(t) 1. In terms of the reduced time τ(t) t 0 dξa(ξ), accounting for arbitrary but given time-dependent infection rates, the SIR-model equations are [1,2,18] dI dτ j − KI, dS dτ −j, dR dτ KI (1) in terms of the time-dependent ratio K(t) μ(t)/a(t) of the recovery and infection rates and the medically interesting daily rate of new infections where the dot denotes a derivative with respect to t.
For the special and important case of a time-independent ratio K(t) k const. new analytical results of the SIR-model (1) have been recently derived [19] hereafter referred to as paper A. The constant k is referred to as the inverse basic reproduction number k 1/R 0 . The new analytical solutions assume that the SIR equations are valid for all times t ∈ [−∞, ∞], and that time t τ 0 refers to the "observing time" when the existence of a pandemic wave in the society is realized and the monitoring of newly infected persons _ J(t) is started. In paper A it has been shown that, for arbitrary but given infection rates a(t), apart from the peak reduced time τ 0 of the rate of new infections, all properties of the pandemic wave as functions of the reduced time are solely controlled by the inverse basic reproduction number k. The dimensionless peak time τ 0 is controlled by k and the value ε −lnS(0), indicating as only initial condition at the observing time the fraction of initially susceptible persons S(0) e −ε . This suggests to introduce the relative reduced time Δ τ − τ 0 with respect to the reduced peak time. In real time t the adopted infection rate a(t) acts as second parameter, and the peak time t m , where _ J(t) reaches its maximum must not coincide with the time, where the reduced j reaches its maximum, i.e., τ m ≡ τ(t m ) ≠ τ 0 , in general.

RESULTS AND DISCUSSION
According to paper A the three fractions of the SIR-model can be expressed in terms of the cumulative number J(τ) and differential daily rate j(τ) dJ/dτ of new infections. The cumulative number satisfies the nonlinear differential equation Two important values are J 0 (k) J(τ 0 ), where j attains its maximum with (dj/dJ) J0 0, and the final cumulative number J ∞ (k) at τ t ∞, when the second bracket on the right-hand side of the differential Eq. 4 vanishes, i.e., J ∞ + kln(1 − J ∞ ) 0. The two transcendental equations can be solved analytically in terms of Lambert's W function, as shown in paper A. In the present manuscript we are going to avoid Lambert's function completely, and instead use the following approximants ( Figure 1A) Without any detailed solution of the SIR-model equations the formal structure of Eqs 3 and 4 then provides the final values I ∞ j ∞ 0, R ∞ J ∞ , and S ∞ 1 − J ∞ . We list these values together with κ in Table 1. We emphasize that the final cumulative number J ∞ , determined solely by the value of k, remains unchanged ( Table 1). With NPIs one can only flatten and distort the epidemics curve (compared to the case of no NPIs taken) but not change the final cumulative number.

New infections
The exact solution of the differential Eq. 4 is given in inverse form by (Appendix A) which can be integrated numerically (subject to numerical precision issues), replaced by the approximant presented in paper A (involving Lambert's function), or semi-quantitatively captured by the simple approximant to be presented next. The solution J(Δ) as a function of the relative reduced time Δ τ − τ 0 , with the reduced peak time approximated by corresponding to J J 0 in Eq. 8, and where f m (k) 1 − k + lnk, is reasonably well captured by (Appendix C) with the Heaviside step function Θ(x) 1(0) for x ≥ (<)0. In Eq. 10 with also tabulated in Table 1. We note that Δ s (k) is always positive. Figure 2 shows the approximation (Eq. 10) for the cumulative number as a function of the relative reduced time Δ for different values of k. For a comparison with the exact variation obtained by the numerical integration of Eq. 8 see Appendix C. The agreement is remarkably well with maximum deviations less than 30 percent. The known limiting case of k 0 is captured exactly by the approximant (Appendix D).
For the corresponding reduced differential rate j(Δ) in reduced time we use the right hand side of Eq. 4 with J J(Δ) from Eq. 10, cf. Figure 3. Note, that this j is not identical with the one obtained via j dJ/dτ, because J does not solve the SIR equations exactly. The peak value j max in the reduced time rate occurs when J J 0 and is thus determined by j max (1 − J 0 )(1 − J 0 − k), also tabulated in Table 1.
As can be seen in Figure 3 the rate of new infections (Eq. 12) is strictly monoexponentially increasing j(Δ ≪ 0)xe Γ1Δ with Γ 1 (k) f m (k)/(1 − k) well before the peak time, and strictly monoexponentially decreasing well above the peak time j(Δ ≫ 0) ∝ e −Γ2Δ with the Γ 2 (1 − J ∞ )Γ 1 /κ. These exponential rates exhibit a noteworthy property and correlation in reduced time: The SIR parameter k affects several key properties of the differential and cumulative fractions of infected persons. If the maximum _ J(t m ) of the measured daily number of newly infected persons has passed already, we find it most convenient to estimate k from the cumulative value J(t m ) at this time t m . While the maximum of _ J(t) must not occur exactly at τ(t m ) τ 0 (Appendix F), we can still use J 0 as an approximant for the value of J(t m ) and the relationship between J 0 and k can be inverted to read (Appendix E)   Frontiers in Physics | www.frontiersin.org January 2021 | Volume 8 | Article 593421 The dependency of k on J 0 is shown in Figure 1C. With the soobtained value for k at hand, the infection rate a(t m ) at peak time can be inferred from a(t m ) _ J(t m )/j max (k). It provides a lower bound for a 0 .
A major advantage of the new analytical solutions in paper A and here is their generality in allowing for arbitrary timedependencies of the infection rate a(t). Such time-dependencies result at times greater than the observing time t 0 from nonpharmaceutical interventions (NPIs) taken after the pandemic outbreak [20] such as case isolation in home, voluntary home quarantine, social distancing, closure of schools and universities and travel restrictions including closure of country borders, applied in different combinations and rigor [21] in many countries. These NPIs lead to a significant reduction of the initial constant infection rate a 0 at later times. It is also important to estimate the influence of a later lifting of the NPIs on the resulting increase in the case numbers in order to discriminate this increase from the onset of a second wave. Especially in the papers by Dehning et al. [17], Flaxman et al. [22] and those reviewed by Estrada [4] the influence of NPIs on the time evolution of the Covid-19 pandemics has been studied using numerical solutions of the SIR-model equations. Our analytical study presented here is superior to these results from numerical simulations as its predictions are particularly robust for the late forecast of the pandemic wave.

Modeling in Real Time of Lockdowns
The corresponding daily rate _ J(t) and cumulative number J(t) of new infections in real time t for given time-dependent infection rates a(t) are J(t) J(τ(t)) and _ J(t) given by Eq. 2. From a medical point of view the daily rate _ J(t) is most important as it determines also i) the fatality rate [23] with the fatality percentage f x0.005 in countries with optimal medical services and hospital capacities and the delay time of t d x7 days, ii) the daily number of new seriously sick persons [24] NSSPs 2f _ J(t − t d ) needing access to breathing apparati, and iii) the day of maximum rush to hospitals t r t m + t d . In countries with poor medical and hospital capacities and/or limited access to them the fatality percentage is significantly higher by a factor h which can be as large as 10.
To calculate the rate and cumulative number in real time according to Eq. 2 we adopt as time-dependent infection rate the integrable function known from shock wave physics which implies The time-dependent lockdown infection rate (Eq. 15) is characterized by four parameters: i) the initial constant infection rate a 0 at early times t ≪ t a , ii) the final constant infection rate a 1 qa 0 at late times t ≫ t a described by the quarantine factor q a 1 /a 0 ≤ 1, first introduced in Refs. 21 and 24, iii) the time t a of maximum change, and iv) the time t b regularizing the sharpness of the transition. The latter is known to be about t b x7-14 days reflecting the typical 1-2 weeks incubation delay. Consequently, the parameter q mainly affects the amplitude _ J shown in the left columns of Figures 5 and 6 (note that we also plotted the case of no NPIs taken (i.e., q 1) for comparison). Alternatively, the transition time t b controls the rapidness of the transition in the fraction of infected persons per day and therefore the widespread.
Moreover, the initial constant infection rate a 0 characterizes the Covid-19 virus: if we adopt the German values a 0 x58 days −1 and t b x11 determined below, with the remaining two parameters q and t a we can represent with the chosen functional form Eq. 15 four basic types of reductions: 1) drastic (small q ≪ 1) and rapid (t a small), 2) drastic (small q ≪ 1) and late (t a large), 3) mild (greater q) and rapid (t a small), and 4) mild (greater q) and late (t a large). The four types are exemplified in Figure 4.

Verification and Forecast
In countries where the peak of the first Covid-19 wave has already passed such as e.g. Germany, Switzerland, Austria, Spain, France and Italy, we may use the monitored fatality rates and peak times to check on the validity of the SIR model with the determined free parameters. However, later monitored data are influenced by a time varying infection rate a(t) resulting from nonpharmaceutical interventions (NPIs) taken during the pandemic evolution. Only at the beginning of the pandemic wave it is justified to adopt a time-independent injection rate a(t)xa 0 implying τ a 0 t. Alternatively, also useful for other countries which still face the climax of the pandemic wave, it is possible to determine the free parameters from the monitored cases in the early phase of the pandemic wave. We illustrate our parameter estimation using the monitored data from Germany with a total population of 83 million persons (P 8.3 × 10 7 ). In Germany the first two deaths were reported on March 9 so that ε 4.8 × 10 − 6 corresponding to about 400 infected people 7 days earlier, on March 2 (t 0). The maximum rate of newly infected fraction, _ J max x380/fP, occurred t m 37 days later, consistent with a peak time of fatalities on 16 April 2020. At peak time the cumulative death number was D m 3820/P corresponding to J m D m /f 200D m 0.009. This implies k ≈ 1 − J 0 ≈ 0.991 according to Eq. 14 (and not far from the value k 0.989 to be determined from the fit shown in Figure 5). From the initial exponential increase of daily fatalities in Germany we extract Γ 1 (k)a 0 x0.28, corresponding to a doubling time of ln(2)/Γ 1 a 0 x2.3 days, as we know Γ 1 (k) f m (k)/(1 − k)x0.0046 already from the above k. The quantity a m we can estimate from the measured _ J max , as _ J max a m j max and j max (k) ≈ 4.2 × 10 − 5 . Using the mentioned value for _ J max , we obtain a m ≈ 22/days as a lower bound for a 0 .
With these parameter values the entire following temporal evolution of the pandemic wave in Germany can be predicted as function of t b and q. To obtain all parameters consistently, we fitted the available data to our model without constraining any of the parameters ( Figure 5). This yields for Germany kx0.989, t a x21 days, q x0.15, a 0 x58 days, and t b 11 days. The obtained parameters allow us to calculate the dimensionless peak time τ m x1353, the dimensionless time τ 0 x1390, as well as J m x0.009, J 0 x0.011, Γ 1 x0.0056 and Γ 2 x0.0055.
We note that the value of k 0.989 implies for Germany that J ∞ (0.989) 0.022 according to Figure 1, so that at the end of the first Covid-19 wave in Germany 2.2% of the population, i.e., 1.83 million persons will be infected. This number corresponds to a final fatality number of D ∞ 9146 persons in Germany.  In each panel we consider four basic types of reductions: 1) drastic (small q 0.1) and rapid (t a 20), 2) drastic (small q 0.1) and late (t a 40), 3) mild (q 0.5) and rapid (t a 20), and 4) mild (q 0.5) and late (t a 40). Remaining parameters due to Germany: t b 11 days, a 0 57 days −1 , and k 0.989. Both curves in (C) for the late cases reach the same value at the maximum, which is plausible as the curve remains unaffected at the time of the maximum. The rapid cases tend to lower the maximum amplitude already at the time of the maximum, and thus tend to decrease it compared with the late cases.
Frontiers in Physics | www.frontiersin.org January 2021 | Volume 8 | Article 593421 Of course, these numbers are only valid estimates if no efficient vaccination against Covid-19 will be available. An important consequence of the small quarantine factor q 0.15 is the implied flat exponential decay after the peak. Because Γ 1 xΓ 2 for k 0.989, the exponential decay is by a factor q smaller than the exponential rise prior the climax, i.e.,  Figure 5 we calculate the resulting daily new infection rate as a function of the time t for the parameters for Germany, and compare with the measured data. In the meantime, the strict lockdown interventions have been lifted in Germany: this does not effect the total numbers J ∞ and D ∞ but it should reduce the half-live decay time further. We also performed this parameter estimation for other countries with sufficient data. For some of them data is visualized in Figure 6, parameters for the remaining countries are tabulated in Tables 2 and 3. Most importantly, with the  [4,22]. Part of these significant differences may be explained by the different definitions of R 0 . While the inverse basic reproduction number k 1/R 0 μ(t)/a(t) in the SIR-model is clearly defined as the ratio of the recovery to infection rate, there are alternative definitions of the basic reproduction number R 0 using the effective reproduction factor R(t). As discussed in detail in Sect. 4 of Ref. 25 R(t) has to be calculated from the convolution of the number of daily cases c(t) with the serial interval distribution W(t) describing the probability for the time lag between a person's infection and the subsequent transmission of the virus to a second person. As different choices of the serial interval distribution are used in the literature this leads to differences in the calculated associated effective reproduction factors R(t). As R 0 is identical to the value R(t 0 ) at the starting time of the outbreak it is not clear in the moment that this R 0 will be identical to the 1/k of the SIR model [26].

Summary and Conclusion
In this work we derived for the first time an analytical approximation for the solution for the SIR-model equations with an accuracy better than 30 percent. The explicit approximation refers to the fraction of newly infected persons per day _ J as a function of the relative reduced time with respect to the reduced peak time. This closed form of the analytical solution only depends on a single parameter k, the ratio of infection to recovery rates. We assume that this ratio is independent of time. As _ J can be directly compared with the monitored death and infection rates in different countries, we see no advantage in using the more complicated SEIR-model which currently does not allow for a closed analytical solution. An analytic solution of the SIR model with an accuracy better than 5% is available as well from our yet unpublished work where we did not consider timedependent a(t), but it has the disadvantage that it involves Lambert's function.
For the first time in the history of SIR-research (these equations have been discovered 93 years ago!) we thus have derived an analytical solution which can be applied successfully to all accumulated data of virus diseases in the world. Being of analytic form it is superior to all existing numerical simulations and results in the literature. We also discovered for the first time how to extract the value of k from the monitored data which is highly nontrivial. We applied this new method to the data taken 2 | Model parameters and model implications. The columns are as follows: country (α 3 code), population P in millions (M), outbreak time defining t 0, fitted time t a , estimated time t 0 corresponding to the reduced peak time τ 0 of j(τ), fitted SIR parameter k, fitted initial infection rate a 0 , fitted parameter t b , fitted quarantine factor q, estimated doubling time t 2 characterizing the early exponential increase, estimated decay half life t ′ 2 characterizing the late exponential decrease, estimated unreported number of infections per reported number, estimated final fraction J ∞ of infected population, final number of estimated fatalities D ∞ P J ∞ Pf . We use f 0.005 as the probability to decease from a Covid-19 infection (reported plus unreported).   In each panel we consider four basic types of reductions: 1) drastic (small q 0.1) and rapid (t a 20), 2) drastic (small q 0.1) and late (t a 40), 3) mild (q 0.5) and rapid (t a 20), and 4) mild (q 0.5) and late (t a 40). Remaining parameters due to Germany: t b 11 days, a 0 57 days −1 , and k 0.989.
for the Covid-19 pandemic waves in many countries. Our work includes an estimate on the effects of non-pharmaceutical interventions in these countries. This is possible as our analytical solution holds for arbitrary but given time dependencies a(t) of the infection rates.
An example, on how lockdown lifting can be modeled is described in Appendix H. The situation is depicted in Figure 7. The lifting will increase a(t) from its present value up to a value that might be close to the initial a 0 . While the dynamics is altered, the final values remain unaffected by the dynamics, except, if the first pandemic wave is followed by a 2nd one. The values for J ∞ provided in Tables 2 and 3 provide a hint on how likely is a 2nd wave. These values correspond to the population fraction that had been infected already. While this fraction is extremely large in Peru (53%), it is still below 1% in several of the larger countries. The tables also report the unreported number of infections per reported number (column "dark"), estimated from the number of fatalities, reported infections, and the death probability f.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://pomber.github.io/covid19/timeseries. json.
approximant G(τ) and J(τ) 1 − e −G(τ) developed in part A, and iii) the approximant J(Δ) given by Eq. 10 with Δ τ − τ m and τ m specified by Eq. 9. Because the numerical solution (i) is extremely well approximated by (ii), and (ii) and (iii) compared to (i) not prone to numerical instabilities at small and large Δ, we present results only for method (iii), as they can be readily reproduced by a reader without Lambert's function at hand.