Skip to main content

ORIGINAL RESEARCH article

Front. Phys., 20 January 2021
Sec. Social Physics
Volume 8 - 2020 | https://doi.org/10.3389/fphy.2020.593421

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

www.frontiersin.orgR. Schlickeiser1,2* www.frontiersin.orgM. Kröger3*
  • 1Institut für Theoretische Physik, Lehrstuhl IV: Weltraum- und Astrophysik, Ruhr-Universität Bochum, Bochum, Germany
  • 2Institut für Theoretische Physik und Astrophysik, Christian-Albrechts-Universität zu Kiel, Kiel, Germany
  • 3Polymer Physics, Department of Materials, ETH Zurich, Zurich, Switzerland

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.

1 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 [317]. The considered population of N1 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)=0tdξa(ξ), accounting for arbitrary but given time-dependent infection rates, the SIR-model equations are [1, 2, 18]

dIdτ=jKI,dSdτ=j,dRdτ=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

J˙(t)=a(t)j(τ)=τ˙j(τ),(2)

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/R0. 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 tm, where J˙(t) reaches its maximum must not coincide with the time, where the reduced j reaches its maximum, i.e., τmτ(tm)τ0, in general.

2 Results and Discussion

According to paper A the three fractions of the SIR-model

S(τ)=1J(τ),I(τ)=j(τ)1J(τ),R(τ)=kln[1J(τ)](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

j(τ)=dJdτ=(1J)[J+kln(1J)](4)

Two important values are J0(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(1J)=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)

J0(k)=(3k)(1k)(1+k+k2)/6,(5)
J(k)=1exp[(1k)(1+κ)/k],(6)
κ(k)=(4k)k/3(7)

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=1J. 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.

FIGURE 1
www.frontiersin.org

FIGURE 1. (A) Approximants J0, J, and κ (thick green) used in this manuscript, cf. Eqs (5)(7), compared with the exact functions [19] (thin black). (B)J vs. J0 and (C)J vs. jmax for the SIR model. (D)k as function of J0 according to Eq. 14. For J00.1 this is well approximated by k1J0, where J0 can be replaced by the cumulative fraction of infected people at the time of the maximum in the daily number of newly infected people.

TABLE 1
www.frontiersin.org

TABLE 1. Exact parameter values depending on the inverse basic reproduction number k.

2.1 New infections

The exact solution of the differential Eq. 4 is given in inverse form by (Appendix A)

τ=1eεJdy(1y)[y+kln(1y)],(8)

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

τ0=1kfm(k)[lnJ01J0ln(eε1)],(9)

corresponding to J=J0 in Eq. 8, and where fm(k)=1k+lnk, is reasonably well captured by (Appendix C)

J(Δ)=12[1+tanhY1(Δ)]Θ[Δs(k)Δ]+{11J2[1+cothY2(Δ)]}Θ[ΔΔs(k)](10)

with the Heaviside step function Θ(x)=1(0) for x(<)0. In Eq. 10

Y1=12[fm(k)(ΔΔs)1k+ln1kk],Y2=12[E0(k)(ΔΔs)+lnk(1k)κ],(11)

with

Δs=1kfm(k)ln(1k)(1J0)kJ0,E0(k)=[k(1k)κ1]fm(k)(12)

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).

FIGURE 2
www.frontiersin.org

FIGURE 2. (A) Cumulative number J(Δ) for different values of k according to Eq. 10. The vertical gray lines starting at the Δ-axis indicate the respective values of Δs(k). (B) Same as in (A), divided by the final J.

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 jmax in the reduced time rate occurs when J=J0 and is thus determined by jmax=(1J0)(1J0k), also tabulated in Table 1.

FIGURE 3
www.frontiersin.org

FIGURE 3. Reduced differential rate j(Δ) of newly infected fraction corresponding to the cumulative J(Δ) shown in Figure 2. (A) linear scale, (B) semilogarithmic scale.

As can be seen in Figure 3 the rate of new infections (Eq. 12) is strictly monoexponentially increasing j(Δ0)eΓ1Δ with Γ1(k)=fm(k)/(1k) well before the peak time, and strictly monoexponentially decreasing well above the peak time j(Δ0)eΓ2Δ with the Γ2=(1J)Γ1/κ. These exponential rates exhibit a noteworthy property and correlation in reduced time:

Γ2Γ1=1Jκ(13)

The SIR parameter k affects several key properties of the differential and cumulative fractions of infected persons. If the maximum J˙(tm) 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(tm) at this time tm. While the maximum of J˙(t) must not occur exactly at τ(tm)=τ0 (Appendix F), we can still use J0 as an approximant for the value of J(tm) and the relationship between J0 and k can be inverted to read (Appendix E)

k=2(1J0)11+ln(1J0)=1J0J022+O(J03)(14)

The dependency of k on J0 is shown in Figure 1C. With the so-obtained value for k at hand, the infection rate a(tm) at peak time can be inferred from a(tm)=J˙(tm)/jmax(k). It provides a lower bound for a0.

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 [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 a0 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.

2.2 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] d(t)fJ˙(ttd) with the fatality percentage f0.005 in countries with optimal medical services and hospital capacities and the delay time of td7 days, ii) the daily number of new seriously sick persons [24] NSSPs=2fJ˙(ttd) needing access to breathing apparati, and iii) the day of maximum rush to hospitals tr=tm+td. 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

aLD(t)=a02[1+q(1q)tanhttatb]{a0for ttaqa0for tta,(15)

which implies

τLD(t)=a02[(1+q)t(1q)tbln(coshttatbcoshtatb)]{a0tfor ttaqa0tfor tta(16)

The time-dependent lockdown infection rate (Eq. 15) is characterized by four parameters: i) the initial constant infection rate a0 at early times tta, ii) the final constant infection rate a1=qa0 at late times tta described by the quarantine factor q=a1/a01, first introduced in Refs. 21 and 24, iii) the time ta of maximum change, and iv) the time tb regularizing the sharpness of the transition. The latter is known to be about tb7–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 tb controls the rapidness of the transition in the fraction of infected persons per day and therefore the widespread.

Moreover, the initial constant infection rate a0 characterizes the Covid-19 virus: if we adopt the German values a058 days−1 and tb11 determined below, with the remaining two parameters q and ta we can represent with the chosen functional form Eq. 15 four basic types of reductions: 1) drastic (small q1) and rapid (ta small), 2) drastic (small q1) and late (ta large), 3) mild (greater q) and rapid (ta small), and 4) mild (greater q) and late (ta large). The four types are exemplified in Figure 4.

FIGURE 4
www.frontiersin.org

FIGURE 4. (A) Infection rate a(t), (B) reduced time τ(t), (C) daily rate of new infections J˙(t), (D) cumulative fraction J(t) of infected persons. In each panel we consider four basic types of reductions: 1) drastic (small q=0.1) and rapid (ta=20), 2) drastic (small q=0.1) and late (ta=40), 3) mild (q=0.5) and rapid (ta=20), and 4) mild (q=0.5) and late (ta=40). Remaining parameters due to Germany: tb=11 days, a0=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.

2.3 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 time-independent injection rate a(t)a0 implying τ=a0t. 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×107).

In Germany the first two deaths were reported on March 9 so that ε=4.8×106 corresponding to about 400 infected people 7 days earlier, on March 2 (t=0). The maximum rate of newly infected fraction, J˙max380/fP, occurred tm=37 days later, consistent with a peak time of fatalities on 16 April 2020. At peak time the cumulative death number was Dm=3820/P corresponding to Jm=Dm/f=200Dm=0.009. This implies k1J00.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)a00.28, corresponding to a doubling time of ln(2)/Γ1a02.3 days, as we know Γ1(k)=fm(k)/(1k)0.0046 already from the above k. The quantity am we can estimate from the measured J˙max, as J˙max=amjmax and jmax(k)4.2×105. Using the mentioned value for J˙max, we obtain am22/days as a lower bound for a0.

FIGURE 5
www.frontiersin.org

FIGURE 5. (A) Measured data J˙(t) of new daily infected fraction (black circles) for Germany (DEU) compared with the model J˙(t)=a(t)j(τ(t)) outlined here (green). Shown for comparison is the case where no NPIs are imposed (q=1, black dot-dashed). (B) The measured cumulative fraction J(t) (black circles) together with the model prediction (green), and the reference q=1 case. Also depicted are J0 and the value at peak time, Jm. (C) The infection rate a(t) corresponding to the curves shown in (A) and (B). Model parameters mentioned in the figure; a dropped from an initial value of 58/days down to 7.8/days during the 2nd half of march. This realized case can be directly compared with the four hypothetic cases shown in Figure 4. For details on how to obtain the parameters see Appendix G.

With these parameter values the entire following temporal evolution of the pandemic wave in Germany can be predicted as function of tb 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 k0.989, ta21 days, q0.15, a058 days, and tb=11 days. The obtained parameters allow us to calculate the dimensionless peak time τm1353, the dimensionless time τ01390, as well as Jm0.009, J00.011, Γ10.0056 and Γ20.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. 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˙(ttm)eΓ2a0qt=eΓ1(1J)a0qtκet/21.8days(17)

Equation 17 yields a decay half-live of ln(2)×21.8 days 15 days to be compared with the initial doubling time of 2.3 days. For Germany we thus know that their lockdown was drastic and rapid: the time ta March 23 is early compared to the peak time tm April 8 resulting in a significant decrease of the infection rate with the quarantine factor qam/a0=0.15. In 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 exception of the six countries ARM, DOM, IRN, PAN, PER, SMR we found values of k>0.9 for all other countries investigated corresponding to basic reproduction numbers R0=1/k<1.11. These values are significantly smaller than the estimates of R0[2.4,5.6] in the mainstream literature on Covid-19 [4, 22]. Part of these significant differences may be explained by the different definitions of R0.

FIGURE 6
www.frontiersin.org

FIGURE 6. Same as Figure 5 for other countries: (A) Italy (ITA), (B) France (FRA), (C) Sweden (SWE), (D) Iran (IRN), (E) Great Britain (GBR).

TABLE 2
www.frontiersin.org

TABLE 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 ta, estimated time t0 corresponding to the reduced peak time τ0 of j(τ), fitted SIR parameter k, fitted initial infection rate a0, fitted parameter tb, fitted quarantine factor q, estimated doubling time t2 characterizing the early exponential increase, estimated decay half life t2 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 DP=JPf. We use f=0.005 as the probability to decease from a Covid-19 infection (reported plus unreported).

TABLE 3
www.frontiersin.org

TABLE 3. Continuation of Table 2.

While the inverse basic reproduction number k=1/R0=μ(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 R0 using the effective reproduction factor R(t). As discussed in detail in Sect. 4 of Ref. 25R(t) has to be calculated from the convolution

R(t)=c(t)0dsW(s)c(ts)(18)

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 R0 is identical to the value R(t0) at the starting time of the outbreak it is not clear in the moment that this R0 will be identical to the 1/k of the SIR model [26].

2.4 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 time-dependent 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 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 a0. 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.

FIGURE 7
www.frontiersin.org

FIGURE 7. Same as Figure 4 but with incomplete lifting (η=0.8) hundred days after breakout (ts=100 days). (A) infection rate a(t), (B) reduced time τ(t), (C) daily rate of new infections J˙(t), (D) cumulative fraction J(t) of infected persons. In each panel we consider four basic types of reductions: 1) drastic (small q=0.1) and rapid (ta=20), 2) drastic (small q=0.1) and late (ta=40), 3) mild (q=0.5) and rapid (ta=20), and 4) mild (q=0.5) and late (ta=40). Remaining parameters due to Germany: tb=11 days, a0=57 days−1, and k=0.989.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://pomber.github.io/covid19/timeseries.json.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Appendix a: non-parametric solution of the sir model

We start from the Eq. 19 from part A

τ=εGdx1exkx(A1)

and substitute

y=1ex,x=ln(1y),dxdy=11y(A2)

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

τ=ψJdy(1y)f(y),f(y)=y+kln(1y)(A3)

with the abbreviation ψ=J(0)=1eε 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.

Appendix a.1: Maximum of j

Taking the derivative of Eq. 37 from part A with respect to τ we obtain

1=j(1J)[J+kln(1J)](A4)

or the exact SIR relation

j=(1J)[J+kln(1J)](A5)

Equation A5 provides

djdJ=1k2Jkln(1J)(A6)

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

1J0=kln(1J0)+k+12(A7)

Setting 1J0=eX yields

eX=k2(Xk+1k),(A8)

which is of the form of Eq. G1 from part A, and solved in terms of the non-principal Lambert function W1 as

X=k+1k+W1(α0),α0=2kee1/k,(A9)

so that

J0=1e1+kkW1(α0)=1+k2W1(α0)(A10)

The maximum value is then given by

jmax=j(J0)=k24{[1+W1(α0)]21}(A11)

and this can also be written as jmax=(1J0)(1J0k) with J0 from Eq. A10. According to Eq. 8 the reduced peak time in the dimensionless rate of new infections is then given by

τ0=ψJ0dy(1y)f(y),(A12)

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

Δ=ττ0=J0Jdy(1y)f(y)(A13)

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 yc+kln(1yc)=0, or 1yc=eyc/k with the solution

yc(k)=J(k)=(1k)(1+κ)(B1)

where κ was already stated in the introduction. According to Eq. A13 the value J corresponds to Δ=τ=, so the maximum value of the cumulative number of new infections is Jmax=J.

Moreover, the function f(y) attains its maximum value fm(k)=f(y=1k)=1k+klnk at ym=1k. As approximation we use

f(y){f1(y)for y1kf2(y)for (1k)yJ=fm{y1kfor y1kJyJ(1k)for 1kyJ(B2)

which is shown in Figure B1 in comparison with the function f(y). The agreement is reasonably well with maximum deviations less than 30%.

FIGURE B1
www.frontiersin.org

FIGURE B1. Comparison of the approximation (B2) with the exact curve for f(y) for different values of k. Vertical solid lines mark the position of the maximum of the function.

Appendix c: approximations for J (τ)

Figure C1 demonstrates that J0(k) is always smaller than 1k. In order to calculate the integral in Eq. A13 with the approximation Eq. B1 we then have to investigate two cases: 1) For J0<1k and J<1k only the function f1 contributes and

Δ(J<1k,J0<1k)=J0Jdy(1y)f1(y)(C1)

2) For J1k>J0 both functions f1 and f2 contribute and

Δ(J1k>J0)J01kdy(1y)f1(y)+1kJdy(1y)f2(y)=Δs+1kJdy(1y)f2(y)(C2)

with

Δs=J01kdy(1y)f1(y)=1kfm(k)J01kdy(1y)y=1kfm(k)ln1J0111k1=1kfm(k)ln(1J0)(1k)kJ0(C3)

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

FIGURE C1
www.frontiersin.org

FIGURE C1. The ratio (1k)/J0(k) as a function of k.

Appendix C.1 Case (1): J ≤ 1 − k, J0 < 1 − k

Here Eq. C1 provides

fmΔ1k=J0Jdy(1y)y=ln1J011J1,(C4)

so that the difference of Eqs C3 and C4 yields

fm(ΔΔs)1k=lnkJ(1k)(1J),(C5)

or after inversion

J(τ)=[1+k1kefm(ΔΔs)1k]1(C6)

Appendix C.2 Case (2): J ≥ 1 − k > J0

Here Eq. C2 with Eq. C3 yields

fmΔ(1k)J01kdy(1y)y+[J(1k)]1kJdy(1y)(Jy)(C7)
=fmΔs+(J1+k)1kJdy(1y)(Jy),

so that

fm(ΔΔs)J(1k)=1kJdy(1y)(Jy)=11Jln11J1J11Jk(C8)
=11Jlnk(JJ)(J1+k)(1J)

After straightforward but tedious algebra we obtain

JJ1J=J(1k)keE,(C9)
E(Δ)=1JJ(1k)fm(ΔΔs)(C10)

and consequently

J(Δ)=1+J11(1k)κkeE(Δ)(C11)

Using the identities 2(1+e2Y)1=1+tanhY and 2(1e2Y)1=1+cothY we combine the results Eqs C6 and C11 to the analytical approximation of the SIR-model equations at all reduced times, stated in Eqs 1012 above. A comparison with the exact numerical solution of the SIR model is provided in Figure C2. The corresponding j(Δ) is obtained from J(Δ) via Eqs A5.

FIGURE C2
www.frontiersin.org

FIGURE C2. Comparison for (A)j(Δ) and (B)J(Δ) between exact solution of the SIR model (green) and the approximant used here (black) for various k. Our approximant [19] in terms of Lambert’s function is shown as well, but coincides with the exact solution (green) in this representation.

Appendix d: si-limit k=0

In the limit k=0Eq. A7 provides J0(k=0)=1/2 so that with limk0fm(k=0)=1 the time scale (Eq. C3) becomes

Δs(k=0)=limk0ln1kk=(D1)

With this result

Y1(k=0)=Δ212limk0[Δs+lnk1k]=Δ212limk0[ln1kk+lnk1k]=Δ2(D2)

Consequently, the cumulative number Eq. 10 and the rate Eq. 4 in this case for all times correctly reduce to

J(Δ,k=0)=12[1+tanhΔ2],j(Δ,k=0)=14cosh2(Δ/2)(D3)

Appendix e: relationship between J0 and K

Here we prove Eq. 14. According to paper A the quantity J0 is given by J0=1eG0 with

G0=1+kk+W1(2e1/kke)=ln(1J0)(E1)

where e denotes Euler’s number and W1 the non-principal solution of Lambert’s equation z=WeW. Equation E1 is of the form x=r+c1W(cecr/β) upon identifying c=1, r=1/k, β=ke/2, and x=[1+ln(1J0)]. From paper A we thus know that ecx=β(xr) holds, or equivalently

(1J0)e=ke2[1ln(1J0)1k]=ke2[1+ln(1J0)]+e2(E2)

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

Appendix f: time of maximum in the measured differential rate J˙(t)

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

0=J¨(tm)=a˙(tm)j(τ(tm))+a2(tm)j(τ(tm))(E1)

or equivalently,

0=[dlnjdτ+a˙a2]t=tm(E2)

From part A we know that

dlnjdτ=12Jkln(1J)k(E3)

and our J0=J(τ0) solves 1=2J0+kln(1J0)+k. That is, j'(τ0)=0. If a does not depend on time, τ0=τ(tm)=a0tm, but this is not generally the case. To find tm and τmτ(tm) one has to solve Eq. E1, or Eq. E2. Equation E2 with Eq. E3 is solved by

Jm=J(τm)=1+k2W1(αm)(E4)

with

αm=2e(1+Cm)/kek,Cm=a˙(tm)a2(tm)(E5)

The corresponding j is, according to Eq. 4,

j(τm)=(1Jm)[Jm+kln(1Jm)](E6)

The smaller Cm, the closer is Jm to J0.

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 infections J(t)=D(ttd)/f=200D(ttd) with td=7 days. From the cumulative value Jm=J(tm) at the time tm of the maximum in J˙(t) we estimate k via Eq. 14 upon assuming JmJ0. Similarly, am=a(tm) is estimated from am=J˙(tm)/jmax(k). These tm, k, am 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 parameters ta>0, tb>0, a0>0, q[0,1], q<η[0,1] and ts>tm. While τ(t) is given by the integrated a(t), we use three strategies to model J(t): i) the numerical solution of the SIR model, ii) the approximant G(τ) and J(τ)=1eG(τ) 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.

Appendix h: modeling of lockdown lifting

Similarly to the lockdown modeling a later lifting of the NPIs can be modeled by adopting the infection rate

a(t)=aLD(t)Θ(tst)+astop(t)Θ(tts)(H1)

where ts denotes the stop time of the lockdown still represented by the infection rate Eq. 15, and where aLD is given by Eq. 15. The infection rate after ts is assumed to be

astop(t)=a0[qs+(ηqs)tanhttstb]{qsa0for t=tsηa0for tts,(H2)

with qs=aLD(ts)/a0 the quarantine factor reached at the time ts of lifting. Together with the reduced time given by Eq. 16 we now find

τ(t)=τLD(t)Θ(tst)+τstop(t)Θ(tts)(H3)

and

τstop(t)=τLD(ts)+qs(tts)a0+(ηqs)a0tbln[cosh(ttstb)](H4)

with τLD(t) from Eq. 16. For the four basic types of Figure 4 we demonstrate in Figure 7 the effect of incomplete lifting.

References

1. Kermack WO, McKendrick AG. A contribution to the mathematical theory of epidemics. Proc R Soc A (1927). 115:700–21.

2. Kendall DG. Deterministic and stochastic epidemics in closed populations. In: Conference third Berkeley symposium on mathematical statistics and probability; 1955 Jul–Aug. Vol. 4. Berkeley, CA: University of California Press (1956). p. 149–65.

3. Hethcode HW. The mathematics of infectious diseases. SIAM Rev (2000) 42:599–653. doi:10.1137/s0036144500371907

4. Estrada E. Covid-19 and sars-cov-2. modeling the present, looking at the future. Phys Rep (2020) 869:1–51. doi:10.1016/j.physrep.2020.07.005

PubMed Abstract | CrossRef Full Text

5. O’Regan SM, Kelly TC, Korobeinikov A, O’Callaghan MJA, Pokrovskii AV. Lyapunov functions for SIR and SIRS epidemic models. Appl Math Lett (2010) 23:446–8. doi:10.1016/j.aml.2009.11.014

6. Satsuma J, Willox R, Ramani A, Grammaticos B, Carstea AS. Extending the SIR epidemic model. Phys Stat Mech Appl (2004) 336:369–75. doi:10.1016/j.physa.2003.12.035

CrossRef Full Text

7. Cadoni M. How to reduce epidemic peaks keeping under control the time-span of the epidemic. Chaos Solitons Fractals (2020) 138:109940. doi:10.1016/j.chaos.2020.109940

CrossRef Full Text

8. Cadoni M, Gaeta G. Size and timescale of epidemics in the sir framework. Phys D (2020) 411:132626. doi:10.1016/j.physd.2020.132626

PubMed Abstract | CrossRef Full Text

9. Chekroun A, Kuniya T. Global threshold dynamics of an infection age-structured SIR epidemic model with diffusion under the Dirichlet boundary condition. J Diff Equ (2020). 269:117–48. doi:10.1016/j.jde.2020.04.046

CrossRef Full Text

10. Imron C, Hariyanto , Yunus M, Surjanto SD, Dewi NAC. Stability and persistence analysis on the epidemic model multi-region multi-patches. J Phys Conf Ser (2019). 1218: 012035. doi:10.1088/1742-6596/1218/1/012035

11. Karaji PT, Nyamoradi N. Analysis of a fractional SIR model with general incidence function. Appl Math Lett (2020). 108:106499. doi:10.1016/j.aml.2020.106499

CrossRef Full Text

12. Mohamadou Y, Halidou A, Kapen PT. A review of mathematical modeling, artificial intelligence and datasets used in the study, prediction and management of Covid-19. Appl Intell (2020). 50, 3913–25. doi:10.1007/s10489-020-01770-9

CrossRef Full Text

13. Samanta S, Sahoo B, Das B. Dynamics of an epidemic system with prey herd behavior and alternative resource to predator. J Phys Math Theor (2019). 52:425601. doi:10.1088/1751-8121/ab264d

14. Sene N. SIR epidemic model with mittag-leffler fractional derivative. Chaos Solitons Fractals (2020). 137:109833. doi:10.1016/j.chaos.2020.109833

CrossRef Full Text

15. Simon M. SIR epidemics with stochastic infectious periods. Stoch Process their Appl (2020). 130:4252–4274. doi:10.1016/j.spa.2019.12.003

CrossRef Full Text

16. Tian C, Zhang Q, Zhang L. Global stability in a networked SIR epidemic model. Appl Math Lett (2020). 107:106444. doi:10.1016/j.aml.2020.106444

CrossRef Full Text

17. Dehning J, Zierenberg J, Spitzner FP, Wibral M, Neto JP, Wilcek M, et al. Inferring change points in the spread of COVID-19 reveals the effectiveness of interventions. Science (2020). 369:eabb9789. doi:10.1126/science.abb9789

CrossRef Full Text

18. Barmparis GD, Tsironis GP. Estimating the infection horizon of covid-19 in eight countries with a data-driven approach. Chaos Solitons Fractals (2020) 135:109842. doi:10.1016/j.chaos.2020.109842

CrossRef Full Text

19. Kröger M, Schlickeiser R. Analytical solution of the SIR-model for the temporal evolution of epidemics. part a: time-independent reproduction factor. J Phys A (2020). 53:505601. doi:10.20944/preprints202007.0416.v1

CrossRef Full Text

20. Ferguson N, Laydon D, Nedjati-Gilani G. Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand. London: Imperial College London (2020) doi:10.25561/77482

CrossRef Full Text

21. Schlickeiser R, Kröger M. Dark numbers and herd immunity of the first covid-19 wave and future social interventions. Epidem Int J (2020). 4:000152. doi:10.23880/eij-16000152

CrossRef Full Text

22. Flaxman S, Mishra S, Mishra S, Gandy A, Unwin HJT, Mellan TA, et al. Estimating the effects of non-pharmaceutical interventions on Covid-19 in Europe. Nature (2020) 584:257–261. doi:10.1038/s41586-020-2405-7

PubMed Abstract | CrossRef Full Text

23. Schüttler J, Schlickeiser R, Schlickeiser F, Kröger M. Covid-19 predictions using a gauss model, based on data from april 2. Physics (2020). 2:197–212. doi:10.3390/physics2020013

CrossRef Full Text

24. Schlickeiser R, Schlickeiser F. A Gaussian model for the time development of the sars-cov-2 corona pandemic disease. predictions for Germany made on 30 March 2020. Physics (2020) 2:164–70. doi:10.3390/physics2020010

CrossRef Full Text

25. Kröger M, Schlickeiser R. Gaussian doubling times and reproduction factors of the Covid-19 pandemic disease. Front Phys (2020). 8:276. doi:10.3389/fphy.2020.00276

CrossRef Full Text

26. Schlickeiser R, Kröger M. First consistent determination of the basic reproduction number for the first Covid-19 wave in 71 countries from the SIR-epidemics model with a constant ratio of recovery to infection rate. Global J Front Res F (2020). 20:37–43. doi:10.3929/ethz-b-000456421

CrossRef Full Text

Keywords: coronavirus (2019-nCoV), statistical analysis, pandemic spreading, time-dependent infection rate, parameter estimation

Citation: Schlickeiser R and Kröger M (2021) Epidemics Forecast From SIR-Modeling, Verification and Calculated Effects of Lockdown and Lifting of Interventions. Front. Phys. 8:593421. doi: 10.3389/fphy.2020.593421

Received: 10 August 2020; Accepted: 23 October 2020;
Published: 20 January 2021.

Edited by:

Sen Pei, Columbia University, United States

Reviewed by:

Ming Tang, East China Normal University, China
Jiannan Wang, Beihang University, China

Copyright © 2021 Schlickeiser and Kröger.. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: R. Schlickeiser, rsch@tp4.ruhr-uni-bochum.de; M. Kröger, mk@mat.ethz.ch

Download