key: cord-0885722-g9zar987
authors: Schlickeiser, R.; Kroger, M.
title: Epidemics forecast from SIR-modeling, verification and calculated effects of lockdown and lifting of interventions
date: 2020-08-15
journal: nan
DOI: 10.1101/2020.08.12.20173294
sha: 77a8da0a592b02f8d466d02a56e7ebe0e2d5d654
doc_id: 885722
cord_uid: g9zar987

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 and the testing of lockdown and lifting interventions as they hold for arbitrary time dependence of the infection rate. Here we present simple analytical approximations for the rate and cumulative number of new infections.

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] 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, 16 

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 infectionṡ

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 SIRmodel (1) have been recently derived 17 -hereafter referred to * rsch@tp4.rub.de (R.S.), mk@mat.ethz.ch (M.K.) 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 per-sonsJ(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 ε = − ln S(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 , whereJ(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. SIR-model results: According to paper A the three fractions of the SIR-model

(3) 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 equation (4) vanishes, i.e.

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 (Fig. 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 Tab. I. New infections: The exact solution of the differential equation (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 semiquantitatively 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 + ln k, 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 Tab. I. We note that ∆ s (k) is always positive. Figure 2 shows the approximation (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.

For the corresponding reduced differential rate j(∆) in reduced time we use the right hand side of Eq. (4) with J = .

CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted August 15, 2020. . https://doi.org/10.1101/2020.08.12.20173294 doi: medRxiv preprint J(∆) from Eq. (10), cf. Fig. 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 Tab. I. As can be seen in Fig. 3 the rate of new infections (12) is strictly monoexponentially increasing j(∆ 0) e Γ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 maximumJ(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 ofJ(t) must not occur exactly at τ (t m ) = τ 0 , 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)

The dependency of k on J 0 is shown in Fig. 1c . With the so-obtained 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 time-dependencies of the infection rate a(t). Such time-dependencies result at times greater than the observing time t = 0 from non-pharmaceutical interventions (NPIs) taken after the pandemic outbreak 18 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 rigour 19 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 orderr to discriminate this incresase from the onset of a second wave. Modeling in real time of lockdowns: The corresponding daily rateJ(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)) andJ(t) given by Eq. (2) . From a medical point of view the daily rateJ(t) is most important as it determines also (i) the fatality rate 20 

with the fatality percentage f 0.005 in countries with optimal medical services and hospital capacities and the delay time of t d 7 days, (ii) the daily number of new seriously sick persons 21 NSSPs = 2fJ(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 (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. 19, 21 , (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 7-14 days reflecting the typical one-two weeks incubation delay.

Moreover, the initial constant infection rate a 0 characterizes the Covid-19 virus: if we adopt the German values . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 15, 2020. a 0 58 days −1 and t b 11 determined below, with the remaining two parameters q and t a we can represent with the chosen functional form (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 Fig. 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 non-pharmaceutical interventions (NPIs) taken during the pandemic evolution. Only at the beginning of the pandemic wave it is justified to adopt a timeindependent injection rate a(t) a 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 380/f P , 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 0 = 0.009. This implies k ≈ 1 − J 0 = 0.991 according to Eq. (14) . From the initial exponential increase of daily fatalities in Germany we extract Γ 1 (k)a 0 0.28, correspond-ing to a doubling time of ln(2)/Γ 1 a 0 2.3 days, as we know Γ 1 (k) = f m (k)/(1 − k) 0.0046 already from the above k. The quantity a m we can estimate from the measuredJ max , aṡ J max = a m j max and j max (k) ≈ 4.2 × 10 −5 . Using the mentioned value forJ 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 (Fig. 5) . This yields for Germany k 0.989, t a 21 days, q 0.15, a 0 58 days, and t b = 11 days. The obtained parameters allow us to calculate the dimensionless peak time τ m 1353, the dimensionless time τ 0 1390, as well as J m 0.009, J 0 0.011, Γ 1 0.0056 and Γ 2 0.0055.

We note that the value of k = 0.989 implies for Germany that J ∞ (0.989) = 0.022 according to Fig. 1 , so that at the end of the first Covid-19 wave in Germany 4 percent of the population, i.e. 1.83 million persons will be infected. This number corresponds to a final fatality number of D ∞ = 9310 persons in Germany. 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 Γ 2 for k = 0.989, the exponential decay is by a factor q smaller than the exponential rise prior the climax, i.e.,J 3 days. For Germany we thus know that their lockdown was drastic and rapid: the time t a March 23 is early compared to the peak time t m Apr 8 resulting in a significant decrease of the infection rate with the quarantine factor q a m /a 0 = 0.15. In Fig. 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 Fig. 7 , parameters for the remaining countries are tabulated in Appendix I.

The presented equations hold for arbitrary a(t). An example, on how lockdown lifting can be modeled is described in Appendix H. The situation is depicted in Fig. 6 . 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 the tables of Appendix I provide a hint on how likely is a 2nd wave. These values correspond to the population fraction that had been infected already. While this . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 15, 2020. 

Consequently, as the cumulative number of new infections is given (see Eq. (A-37)) by

with the abbreviation ψ = J(0) = 1 − e −ε for the initial value. This inverse relation τ (J) is the general solution of the SIR-model for constant k. It is not in parametrized form.

Taking the derivative of Eq. (A3) with respect to τ we obtain

or the exact SIR relation

The maximum value j max occurs for (dj/dJ) J0 = 0 providing

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 15, 2020. which is of the form (A-G1), and solved in terms of the nonprincipal Lambert function W −1 as

so that

The maximum value is then given by

and this can also be written as

with J 0 from Eq. (A10). According to Eq. (4) the reduced peak time in the dimensionless rate of new infections is then given by

which is the only quantity depending besides on k also on ε via ψ = 1 − e −ε . In order to have one approximation depending only on k we therefore introduce the relative reduced time with respect to the peak reduced time

which is still exact, independent of ε and only determined by the value of k.

Appendix B: Approximating the function f (y)

The function f (y) defined in Eq. (A3) vanishes for y c + k ln(1 − y c ) = 0, or 1 − y c = e −yc/k with the solution

where κ was already stated in the introduction. According to Eq. (A14) the value J ∞ corresponds to ∆ = τ = ∞, so the maximum value of the cumulative number of new infections is J max = J ∞ . Moreover, the function f (y) attains its maximum value f m (k) = f (y = 1 − k) = 1 − k + k ln k at y m = 1 − k.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 15, 2020. .

which is shown in Fig. 8 in comparison with the function f (y).

The agreement is reasonably well with maximum deviations less than 30 percent. 

denoting the relative time corresponding to the value J = 1 − k. We consider each case in turn.

so that the difference of Eqs. (C3) and (C4) yields

or after inversion

Here Eq. (C2) with Eq. (C3) yields

so that

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 15, 2020. . https://doi.org/10.1101/2020.08.12.20173294 doi: medRxiv preprint After straightforward but tedious algebra we obtain

and consequently

Using the identities 2(1 + e −2Y ) −1 = 1 + tanh Y and 2(1 − e −2Y ) −1 = 1 + coth Y we combine the results (C6) and (C11) to the analytical approximation of the SIR-model equations at all reduced times, stated in Eqs. (10) , (11) , and (12) . A comparison with the exact numerical solution of the SIR model is provide in Fig. 10 . The corresponding j(∆) is obtained from J(∆) via Eq. (A5).

Appendix D: SI-limit k = 0

In the limit k = 0 Eq. (A7) provides J 0 (k = 0) = 1/2 so that with lim k→0 f m (k = 0) = 1 the time scale (C3) becomes

With this result

Consequently, the cumulative number (10) Here we prove Eq. (14) . According to paper A the quantity J 0 is given by J 0 = 1 − e −G0 with

where e denotes Euler's number and W −1 the non-principal solution of Lambert's equation z = W e W . Equation (E1) is of the form x = r + c −1 W (ce −cr /β) upon identifying c = 1, r = 1/k, β = −ke/2, and x = −[1+ln(1−J 0 )]. From paper A we thus know that e −cx = β(x − r) holds, or equivalently

This is readily solved for k, and thus proves Eq. (14) .

Appendix F: Time of maximum in the measured differential rateJ(t)

One has J(t) = J(τ (t)) andJ(t) =τ (t)J (τ (t)) = a(t)j(τ (t)) since j = J if we let the prime denote a derivative with respect to τ . The maximum inJ(t) thus fulfills

From part A we know that

and our J 0 = J(τ 0 ) solves 1 = 2J 0 + k ln(1 − J 0 ) + k. That is, j (τ 0 ) = 0. If a does not depend on time, τ 0 = τ (t m ) = a 0 t m , but this is not generally the case. To find t m and τ m ≡ τ (t m ) one has to solve Eq. (F1), or Eq. (F2). Eq. (F2) with (F3) is solved by

The corresponding j is, according to Eq. (4),

The smaller C m , the closer is J m to J 0 .

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted August 15, 2020. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted August 15, 2020. . https://doi.org/10.1101/2020.08.12.20173294 doi: medRxiv preprint country P/M t = 0 ta tm t0 k a0 t b q t2 t 2 dark J∞ D∞P

Proc. Third Berkeley Symp

Stability and persistence analysis on the epidemic model multi-region multi-patches

Impact of nonpharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand

0.931 10.1 15 0.03 1.9 73.7 5.1 13.50% 4763 SWE 9

Appendix G: Fitting the data As discussed in length in paper A we base our analysis of existing data on the reported cumulative number of deaths, D(t), from which we estimate the cumulative number of in-days. From the cumulative value J m = J(t m ) at the time t m of the maximum inJ(t) we estimate k via Eq. (14) upon assuming J m ≈ J 0 . Similarly, a m = a(t m ) is estimated from a m =J(t m )/j max (k). These t m , k, a m are not the final values, but provide starting values which are then used in the minimization of the deviation between measured and modeled J(t). The minimization is performed assuming the time-dependent a(t) parameterized by Eq. (H1) involving pa- 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.

Similarly to the lockdown modeling a later lifting of the NPIs can be modeled by adopting the infection ratewhere t s denotes the stop time of the lockdown still represented by the infection rate (15) , and where a LD is given by Eq. (15) . The infection rate after t s is assumed to bewith q s = a LD (t s )/a 0 the quarantine factor reached at the time t s of lifting. Together with the reduced time (15) we now findandwith τ LD (t) from Eq. (16) . For the four basic types of Fig. 4 we demonstrate in Fig. 6 the effect of incomplete lifting.

We performed the parameter estimation for additional countries with sufficient data on 7th Aug 2020. Tables II-III we report these results.