Epidemics forecast from sir-modeling, verification and calculated effects of lockdown and lifting of interventions

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 [ 3 – 17 ]. The considered population of $N\gg 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\left(t\right)$) and recovery rates ($\mu \left(t\right)$) which in general vary with time due to non-pharmaceutical interventions taken during the pandemic evolution.

Let $I\left(t\right)=i\left(t\right)/N$,$S\left(t\right)=s\left(t\right)/N$ and $R\left(t\right)=r\left(t\right)/N$ denote the infected, susceptible and recovered/removed fractions of persons involved in the infection at time t , with the sum requirement $I\left(t\right)+S\left(t\right)+R\left(t\right)=1$. In terms of the reduced time $\tau \left(t\right)={{\int }^{\text{}}}_0^{t}d\xi a\left(\xi \right)$, accounting for arbitrary but given time-dependent infection rates, the SIR-model equations are [ 1 , 2 , 18 ]

$\frac{dI}{d\tau }=j-KI,\frac{dS}{d\tau }=-j,\frac{dR}{d\tau }=KI\left(1\right)$

in terms of the time-dependent ratio $K\left(t\right)=\mu \left(t\right)/a\left(t\right)$ of the recovery and infection rates and the medically interesting daily rate of new infections

$\dot{J}\left(t\right)=a\left(t\right)j\left(\tau \right)=\dot{\tau }j\left(\tau \right),\left(2\right)$

where the dot denotes a derivative with respect to t .

For the special and important case of a time-independent ratio $K\left(t\right)=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\in [-\infty ,\infty ]$, and that time $t=\tau =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 $\dot{J}\left(t\right)$ is started. In paper A it has been shown that, for arbitrary but given infection rates $a\left(t\right)$, apart from the peak reduced time ${\tau }_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 ${\tau }_0$ is controlled by k and the value $\epsilon =-\text{ln}S\left(0\right)$, indicating as only initial condition at the observing time the fraction of initially susceptible persons $S\left(0\right)=e^{-\epsilon }$. This suggests to introduce the relative reduced time $\text{Δ}=\tau -{\tau }_0$ with respect to the reduced peak time. In real time t the adopted infection rate $a\left(t\right)$ acts as second parameter, and the peak time $t_m$, where $\dot{J}\left(t\right)$ reaches its maximum must not coincide with the time, where the reduced j reaches its maximum, i. e.,${\tau }_m\equiv \tau \left(t_m\right)\ne {\tau }_0$, in general.

2 Results and Discussion

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

$S\left(\tau \right)=1-J\left(\tau \right),I\left(\tau \right)=\frac{j\left(\tau \right)}{1-J\left(\tau \right)},R\left(\tau \right)=-k\text{ln}[1-J\left(\tau \right)]\left(3\right)$

can be expressed in terms of the cumulative number $J\left(\tau \right)$ and differential daily rate $j\left(\tau \right)=dJ/d\tau$ of new infections. The cumulative number satisfies the nonlinear differential equation

$j\left(\tau \right)=\frac{dJ}{d\tau }=\left(1-J\right)[J+k\text{ln}\left(1-J\right)]\left(4\right)$

Two important values are $J_0\left(k\right)=J\left({\tau }_0\right)$, where j attains its maximum with ${\left(dj/dJ\right)}_{J_0}=0$, and the final cumulative number $J_{\infty }\left(k\right)$ at $\tau =t=\infty$, when the second bracket on the right-hand side of the differential Eq. 4 vanishes, i. e.,$J_{\infty }+k\text{ln}\left(1-J_{\infty }\right)=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 )

$J_0\left(k\right)=\left(3-k\right)\left(1-k\right)\left(1+k+k^2\right)/6,\left(5\right)$$J_{\infty }\left(k\right)=1-\text{exp}[-\left(1-k\right)\left(1+\kappa \right)/k],\left(6\right)$$\kappa \left(k\right)=\left(4-k\right)k/3\left(7\right)$

Without any detailed solution of the SIR-model equations the formal structure of Eqs 3 and 4 then provides the final values $I_{\infty }=j_{\infty }=0$,$R_{\infty }=J_{\infty }$, and $S_{\infty }=1-J_{\infty }$. We list these values together with κ in Table 1 . We emphasize that the final cumulative number $J_{\infty }$, 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

(A)Approximants $J_0$,$J_{\infty }$, and κ (thick green) used in this manuscript, cf. Eqs (5)–(7), compared with the exact functions [ 19 ] (thin black).(B)$J_{\infty }$ vs.$J_0$ and(C)$J_{\infty }$ vs.$j_{\text{max}}$ for the SIR model.(D) k as function of $J_0$ according to Eq. 14 . For $J_0\le \mathrm{0.\; 1}$ this is well approximated by $k\approx 1-J_0$, where $J_0$ 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

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 )

$\tau ={\int }_{1-e^{-\epsilon }}^J\frac{dy}{\left(1-y\right)[y+k\text{ln}\left(1-y\right)]},\left(8\right)$

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\left(\text{Δ}\right)$ as a function of the relative reduced time $\text{Δ}=\tau -{\tau }_0$, with the reduced peak time approximated by

${\tau }_0=\frac{1-k}{f_m\left(k\right)}\left[\text{ln}\frac{J_0}{1-J_0}-\text{ln}\left(e^{\epsilon }-1\right)\right],\left(9\right)$

corresponding to $J=J_0$ in Eq. 8 , and where $f_m\left(k\right)=1-k+\text{ln}k$, is reasonably well captured by ( Appendix C )

$J\left(\text{Δ}\right)=\frac{1}{2}\left[1+\text{tanh}{Y}_1\left(\text{Δ}\right)\right]\text{Θ}\left[{\text{Δ}}_s\left(k\right)-\text{Δ}\right]+\left\{1-\frac{1-J_{\infty }}{2}\left[1+\text{coth}{Y}_2\left(\text{Δ}\right)\right]\right\}\text{Θ}\left[\text{Δ}-{\text{Δ}}_s\left(k\right)\right]\left(10\right)$

with the Heaviside step function $\text{Θ}\left(x\right)=1\left(0\right)$ for $x\ge (<)0$. In Eq. 10

$Y_1=\frac{1}{2}\left[\frac{f_m\left(k\right)\left(\text{Δ}-{\text{Δ}}_s\right)}{1-k}+\text{ln}\frac{1-k}{k}\right],Y_2=\frac{1}{2}\left[E_0\left(k\right)\left(\text{Δ}-{\text{Δ}}_s\right)+\text{ln}\frac{k}{\left(1-k\right)\kappa }\right],\left(11\right)$

with

${\text{Δ}}_s=\frac{1-k}{f_m\left(k\right)}\text{ln}\frac{\left(1-k\right)\left(1-J_0\right)}{k{J}_0},E_0\left(k\right)=\left[\frac{k}{\left(1-k\right)\kappa }-1\right]f_m\left(k\right)\left(12\right)$

also tabulated in Table 1 . We note that ${\text{Δ}}_s\left(k\right)$ is always positive. Figure 2 shows the approximation ( Eq. 10 ) for the cumulative number as a function of the relative reduced time $\text{Δ}$ 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

(A)Cumulative number $J\left(\text{Δ}\right)$ for different values of k according to Eq. 10 . The vertical gray lines starting at the $\text{Δ}$-axis indicate the respective values of ${\text{Δ}}_s\left(k\right)$.(B)Same as in(A), divided by the final $J_{\infty }$.

For the corresponding reduced differential rate $j\left(\text{Δ}\right)$ in reduced time we use the right hand side of Eq. 4 with $J=J\left(\text{Δ}\right)$ from Eq. 10 , cf. Figure 3 . Note, that this j is not identical with the one obtained via $j=dJ/d\tau$, because J does not solve the SIR equations exactly. The peak value $j_{\text{max}}$ in the reduced time rate occurs when $J=J_0$ and is thus determined by $j_{\text{max}}=\left(1-J_0\right)\left(1-J_0-k\right)$, also tabulated in Table 1 .

FIGURE 3

Reduced differential rate $j\left(\text{Δ}\right)$ of newly infected fraction corresponding to the cumulative $J\left(\text{Δ}\right)$ 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(\text{Δ}\ll 0)\simeq e^{{\text{Γ}}_1\text{Δ}}$ with ${\text{Γ}}_1\left(k\right)=f_m\left(k\right)/\left(1-k\right)$ well before the peak time, and strictly monoexponentially decreasing well above the peak time $j(\text{Δ}\gg 0)\propto e^{-{\text{Γ}}_2\text{Δ}}$ with the ${\text{Γ}}_2=\left(1-J_{\infty }\right){\text{Γ}}_1/\kappa$. These exponential rates exhibit a noteworthy property and correlation in reduced time:

$\frac{{\text{Γ}}_2}{{\text{Γ}}_1}=\frac{1-J_{\infty }}{\kappa }\left(13\right)$

The SIR parameter k affects several key properties of the differential and cumulative fractions of infected persons. If the maximum $\dot{J}\left(t_m\right)$ 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\left(t_m\right)$ at this time $t_m$. While the maximum of $\dot{J}\left(t\right)$ must not occur exactly at $\tau \left(t_m\right)={\tau }_0$( Appendix F ), we can still use $J_0$ as an approximant for the value of $J\left(t_m\right)$ and the relationship between $J_0$ and k can be inverted to read ( Appendix E )

$k=\frac{2\left(1-J_0\right)-1}{1+\text{ln}\left(1-J_0\right)}=1-J_0-\frac{J_0^2}{2}+O\left(J_0^3\right)\left(14\right)$

The dependency of k on $J_0$ is shown in Figure 1C . With the so-obtained value for k at hand, the infection rate $a\left(t_m\right)$ at peak time can be inferred from $a\left(t_m\right)=\dot{J}\left(t_m\right)/j_{\text{max}}\left(k\right)$. 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\left(t\right)$. 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 $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.

2. 2 Modeling in Real Time of Lockdowns

The corresponding daily rate $\dot{J}\left(t\right)$ and cumulative number $J\left(t\right)$ of new infections in real time t for given time-dependent infection rates $a\left(t\right)$ are $J\left(t\right)=J\left(\tau \left(t\right)\right)$ and $\dot{J}\left(t\right)$ given by Eq. 2 . From a medical point of view the daily rate $\dot{J}\left(t\right)$ is most important as it determines also i) the fatality rate [ 23 ]$d\left(t\right)\simeq f\dot{J}\left(t-t_d\right)$ with the fatality percentage $f\simeq \mathrm{0.\; 005}$ in countries with optimal medical services and hospital capacities and the delay time of $t_d\simeq 7$ days, ii) the daily number of new seriously sick persons [ 24 ]$\text{NSSPs}=2f\dot{J}\left(t-t_d\right)$ 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

$a_{\text{LD}}\left(t\right)=\frac{a_0}{2}\left[1+q-\left(1-q\right)\text{tanh}\frac{t-t_a}{t_b}\right]\simeq \{\begin{array}{cc}{a}_0& \text{for}t\ll t_a\\ q{a}_0& \text{for}t\gg t_a,\end{array}\left(15\right)$

which implies

${\tau }_{\text{LD}}\left(t\right)=\frac{a_0}{2}\left[\left(1+q\right)t-\left(1-q\right)t_b\text{ln}\left(\frac{\text{cosh}\frac{t-t_a}{t_b}}{\text{cosh}\frac{t_a}{t_b}}\right)\right]\simeq \{\begin{array}{cc}{a}_{0}t& \text{for}t\ll t_a\\ q{a}_{0}t& \text{for}t\gg t_a\end{array}\left(16\right)$

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\ll t_a$, ii) the final constant infection rate $a_1=q{a}_0$ at late times $t\gg t_a$ described by the quarantine factor $q=a_1/a_0\le 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\simeq 7$–14 days reflecting the typical 1–2 weeks incubation delay. Consequently, the parameter q mainly affects the amplitude $\dot{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\simeq 58$ days ⁻¹ and $t_b\simeq 11$ 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\ll 1$) and rapid ($t_a$ small), 2) drastic (small $q\ll 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 .

FIGURE 4

(A)Infection rate $a\left(t\right)$,(B)reduced time $\tau \left(t\right)$,(C)daily rate of new infections $\dot{J}\left(t\right)$,(D)cumulative fraction $J\left(t\right)$ of infected persons. In each panel we consider four basic types of reductions: 1) drastic (small $q=\mathrm{0.\; 1}$) and rapid ($t_a=20$), 2) drastic (small $q=\mathrm{0.\; 1}$) and late ($t_a=40$), 3) mild ($q=\mathrm{0.\; 5}$) and rapid ($t_a=20$), and 4) mild ($q=\mathrm{0.\; 5}$) and late ($t_a=40$). Remaining parameters due to Germany:$t_b=11$ days,$a_0=57$ days ⁻¹ , and $k=\mathrm{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\left(t\right)$ 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\left(t\right)\simeq a_0$ implying $\tau =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=\mathrm{8.\; 3}\times {10}^7$).

In Germany the first two deaths were reported on March 9 so that $\epsilon =\mathrm{4.\; 8}\times {10}^{-6}$ corresponding to about 400 infected people 7 days earlier, on March 2 ($t=0$). The maximum rate of newly infected fraction,${\dot{J}}_{\text{max}}\simeq 380/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_{\mathrm{m}}=\mathrm{0.\; 009}$. This implies $k\approx 1-J_0\approx \mathrm{0.\; 991}$ according to Eq. 14 (and not far from the value $k=\mathrm{0.\; 989}$ to be determined from the fit shown in Figure 5 ). From the initial exponential increase of daily fatalities in Germany we extract ${\text{Γ}}_1\left(k\right)a_0\simeq \mathrm{0.\; 28}$, corresponding to a doubling time of $\text{ln}\left(2\right)/{\text{Γ}}_1{a}_0\simeq \mathrm{2.\; 3}$ days, as we know ${\text{Γ}}_1\left(k\right)=f_m\left(k\right)/\left(1-k\right)\simeq \mathrm{0.\; 0046}$ already from the above k . The quantity $a_m$ we can estimate from the measured ${\dot{J}}_{\text{max}}$, as ${\dot{J}}_{\text{max}}=a_m{j}_{\text{max}}$ and $j_{\text{max}}\left(k\right)\approx \mathrm{4.\; 2}\times {10}^{-5}$. Using the mentioned value for ${\dot{J}}_{\text{max}}$, we obtain $a_m\approx 22/$ days as a lower bound for $a_0$.

FIGURE 5

(A)Measured data $\dot{J}\left(t\right)$ of new daily infected fraction (black circles) for Germany (DEU) compared with the model $\dot{J}\left(t\right)=a\left(t\right)j\left(\tau \left(t\right)\right)$ 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\left(t\right)$(black circles) together with the model prediction (green), and the reference $q=1$ case. Also depicted are $J_0$ and the value at peak time,$J_m$.(C)The infection rate $a\left(t\right)$ corresponding to the curves shown in(A)and(B). Model parameters mentioned in the figure; a dropped from an initial value of $\simeq 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 $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 $k\simeq \mathrm{0.\; 989}$,$t_a\simeq 21$ days,$q\simeq \mathrm{0.\; 15}$,$a_0\simeq 58$ days, and $t_b=11$ days. The obtained parameters allow us to calculate the dimensionless peak time ${\tau }_m\simeq 1353$, the dimensionless time ${\tau }_0\simeq 1390$, as well as $J_m\simeq \mathrm{0.\; 009}$,$J_0\simeq \mathrm{0.\; 011}$,${\text{Γ}}_1\simeq \mathrm{0.\; 0056}$ and ${\text{Γ}}_2\simeq \mathrm{0.\; 0055}$.

We note that the value of $k=\mathrm{0.\; 989}$ implies for Germany that $J_{\infty }\left(\mathrm{0.\; 989}\right)=\mathrm{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_{\infty }=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=\mathrm{0.\; 15}$ is the implied flat exponential decay after the peak. Because ${\text{Γ}}_1\simeq {\text{Γ}}_2$ for $k=\mathrm{0.\; 989}$, the exponential decay is by a factor q smaller than the exponential rise prior the climax, i. e.,

$\dot{J}(t\gg t_m)\propto e^{-{\text{Γ}}_2{a}_{0}qt}=e^{-\frac{{\text{Γ}}_1\left(1-J_{\infty }\right)a_{0}qt}{\kappa }}\simeq e^{-t/\mathrm{21.\; 8}\text{days}}\left(17\right)$

Equation 17 yields a decay half-live of $\text{ln}\left(2\right)\times \mathrm{21.\; 8}$ days $\simeq 15$ days to be compared with the initial doubling time of $\simeq \mathrm{2.\; 3}$ days. For Germany we thus know that their lockdown was drastic and rapid: the time $t_a\simeq$ March 23 is early compared to the peak time $t_m\simeq$ April 8 resulting in a significant decrease of the infection rate with the quarantine factor $q\simeq a_m/a_0=\mathrm{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_{\infty }$ and $D_{\infty }$ 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>\mathrm{0.\; 9}$ for all other countries investigated corresponding to basic reproduction numbers $R_0=1/k<\mathrm{1.\; 11}$. These values are significantly smaller than the estimates of $R_0\in [\mathrm{2.\; 4},\mathrm{5.\; 6}]$ in the mainstream literature on Covid-19 [ 4 , 22 ]. Part of these significant differences may be explained by the different definitions of $R_0$.

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

Model parameters and model implications. The columns are as follows: country (${\alpha }_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 ${\tau }_0$ of $j\left(\tau \right)$, 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^{\prime }$ characterizing the late exponential decrease, estimated unreported number of infections per reported number, estimated final fraction $J_{\infty }$ of infected population, final number of estimated fatalities $D_{\infty }P=J_{\infty }Pf$. We use $f=\mathrm{0.\; 005}$ as the probability to decease from a Covid-19 infection (reported plus unreported).

TABLE 3

Continuation of Table 2 .

While the inverse basic reproduction number $k=1/R_0=\mu \left(t\right)/a\left(t\right)$ 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\left(t\right)$. As discussed in detail in Sect. 4 of Ref. 25 $R\left(t\right)$ has to be calculated from the convolution

$R\left(t\right)=\frac{c\left(t\right)}{{{\int }^{\text{}}}_0^{\infty }dsW\left(s\right)c\left(t-s\right)}\left(18\right)$

of the number of daily cases $c\left(t\right)$ with the serial interval distribution $W\left(t\right)$ 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\left(t\right)$. As $R_0$ is identical to the value $R\left(t_0\right)$ 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 ].

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 $\dot{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 $\dot{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\left(t\right)$, 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\left(t\right)$ 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\left(t\right)$ 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_{\infty }$ 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\text{\%}$ 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

Same as Figure 4 but with incomplete lifting ($\eta =\mathrm{0.\; 8}$) hundred days after breakout ($t_s=100$ days).(A)infection rate $a\left(t\right)$,(B)reduced time $\tau \left(t\right)$,(C)daily rate of new infections $\dot{J}\left(t\right)$,(D)cumulative fraction $J\left(t\right)$ of infected persons. In each panel we consider four basic types of reductions: 1) drastic (small $q=\mathrm{0.\; 1}$) and rapid ($t_a=20$), 2) drastic (small $q=\mathrm{0.\; 1}$) and late ($t_a=40$), 3) mild ($q=\mathrm{0.\; 5}$) and rapid ($t_a=20$), and 4) mild ($q=\mathrm{0.\; 5}$) and late ($t_a=40$). Remaining parameters due to Germany:$t_b=11$ days,$a_0=57$ days ⁻¹ , and $k=\mathrm{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

$\tau ={{\int }^{\text{}}}_{\epsilon }^G\frac{dx}{1-e^{-x}-kx}\left(\mathrm{A1}\right)$

and substitute

$y=1-e^{-x},x=-\text{ln}\left(1-y\right),\frac{dx}{dy}=\frac{1}{1-y}\left(\mathrm{A2}\right)$

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

$\tau ={{\int }^{\text{}}}_{\psi }^J\frac{dy}{\left(1-y\right)f\left(y\right)},f\left(y\right)=y+k\text{ln}\left(1-y\right)\left(\mathrm{A3}\right)$

with the abbreviation $\psi =J\left(0\right)=1-e^{-\epsilon }$ for the initial value. This inverse relation $\tau \left(J\right)$ 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=\frac{j}{\left(1-J\right)[J+k\text{ln}\left(1-J\right)]}\left(\mathrm{A4}\right)$

or the exact SIR relation

$j=\left(1-J\right)[J+k\text{ln}\left(1-J\right)]\left(\mathrm{A5}\right)$

Equation A5 provides

$\frac{dj}{dJ}=1-k-2J-k\text{ln}\left(1-J\right)\left(\mathrm{A6}\right)$

The maximum value $j_{\text{max}}$ occurs for ${\left(dj/dJ\right)}_{J_0}=0$ providing

$1-J_0=\frac{k\text{ln}\left(1-J_0\right)+k+1}{2}\left(\mathrm{A7}\right)$

Setting $1-J_0=e^{-X}$ yields

$e^{-X}=-\frac{k}{2}\left(X-\frac{k+1}{k}\right),\left(\mathrm{A8}\right)$

which is of the form of Eq. G1 from part A, and solved in terms of the non-principal Lambert function $W_{-1}$ as

$X=\frac{k+1}{k}+W_{-1}\left({\alpha }_0\right),{\alpha }_0=-\frac{2}{ke}{e}^{-1/k},\left(\mathrm{A9}\right)$

so that

$J_0=1-e^{-\frac{1+k}{k}-W_{-1}\left({\alpha }_0\right)}=1+\frac{k}{2}{W}_{-1}\left({\alpha }_0\right)\left(\mathrm{A10}\right)$

The maximum value is then given by

$j_{\text{max}}=j\left(J_0\right)=\frac{k^2}{4}\left\{{[1+W_{-1}\left({\alpha }_0\right)]}^2-1\right\}\left(\mathrm{A11}\right)$

and this can also be written as $j_{\text{max}}=\left(1-J_0\right)\left(1-J_0-k\right)$ with $J_0$ from Eq. A10 . According to Eq. 8 the reduced peak time in the dimensionless rate of new infections is then given by

${\tau }_0={{\int }^{\text{}}}_{\psi }^{J_0}\frac{dy}{\left(1-y\right)f\left(y\right)},\left(\mathrm{A12}\right)$

which is the only quantity depending besides on k also on ε via $\psi =1-e^{-\epsilon }$. In order to have our approximation depending only on k we therefore introduce the relative reduced time with respect to the peak reduced time

$\text{Δ}=\tau -{\tau }_0={{\int }^{\text{}}}_{J_0}^J\frac{dy}{\left(1-y\right)f\left(y\right)}\left(\mathrm{A13}\right)$

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

Appendix b: approximating the function ${f}\left(y\right)$

The function $f\left(y\right)$ defined in Eq. A3 vanishes for $y_c+k\text{ln}\left(1-y_c\right)=0$, or $1-y_c=e^{-y_c/k}$ with the solution

$y_c\left(k\right)=J_{\infty }\left(k\right)=\left(1-k\right)\left(1+\kappa \right)\left(\mathrm{B1}\right)$

where κ was already stated in the introduction. According to Eq. A13 the value $J_{\infty }$ corresponds to $\text{Δ}=\tau =\infty$, so the maximum value of the cumulative number of new infections is $J_{\text{max}}=J_{\infty }$.

Moreover, the function $f\left(y\right)$ attains its maximum value $f_m\left(k\right)=f(y=1-k)=1-k+k\text{ln}k$ at $y_m=1-k$. As approximation we use

$f\left(y\right)\simeq \{\begin{array}{cc}{f}_1\left(y\right)& \text{for}y\le 1-k\\ f_2\left(y\right)& \text{for}\left(1-k\right)\le y\le J_{\infty }\end{array}=f_m\{\begin{array}{ l}\frac{y}{1-k}& \text{for}y\le 1-k\\ \frac{J_{\infty }-y}{J_{\infty }-\left(1-k\right)}& \text{for}1-k\le y\le J_{\infty }\end{array}\left(\mathrm{B2}\right)$

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

FIGURE B1

Comparison of the approximation (B2) with the exact curve for $f\left(y\right)$ 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 $J_0\left(k\right)$ is always smaller than $1-k$. In order to calculate the integral in Eq. A13 with the approximation Eq. B1 we then have to investigate two cases: 1) For $J_0<1-k$ and $J<1-k$ only the function $f_1$ contributes and

${\text{Δ}}_{(J<1-k,J_0<1-k)}={{\int }^{\text{}}}_{J_0}^J\frac{dy}{\left(1-y\right)f_1\left(y\right)}\left(\mathrm{C1}\right)$

2) For $J\ge 1-k>J_0$ both functions $f_1$ and $f_2$ contribute and

${\text{Δ}}_{(J\ge 1-k>J_0)}\simeq {{\int }^{\text{}}}_{J_0}^{1-k}\frac{dy}{\left(1-y\right)f_1\left(y\right)}+{{\int }^{\text{}}}_{1-k}^J\frac{dy}{\left(1-y\right)f_2\left(y\right)}={\text{Δ}}_s+{{\int }^{\text{}}}_{1-k}^J\frac{dy}{\left(1-y\right)f_2\left(y\right)}\left(\mathrm{C2}\right)$

with

${\text{Δ}}_s={{\int }^{\text{}}}_{J_0}^{1-k}\frac{dy}{\left(1-y\right)f_1\left(y\right)}=\frac{1-k}{f_m\left(k\right)}{{\int }^{\text{}}}_{J_0}^{1-k}\frac{dy}{\left(1-y\right)y}=\frac{1-k}{f_m\left(k\right)}\text{ln}\frac{\frac{1}{J_0}-1}{\frac{1}{1-k}-1}=\frac{1-k}{f_m\left(k\right)}\text{ln}\frac{\left(1-J_0\right)\left(1-k\right)}{k{J}_0}\left(\mathrm{C3}\right)$

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

FIGURE C1

The ratio $\left(1-k\right)/J_0\left(k\right)$ as a function of k .

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

Here Eq. C1 provides

$\frac{f_m\text{Δ}}{1-k}={{\int }^{\text{}}}_{J_0}^J\frac{dy}{\left(1-y\right)y}=\text{ln}\frac{\frac{1}{J_0}-1}{\frac{1}{J}-1},\left(\mathrm{C4}\right)$

so that the difference of Eqs C3 and C4 yields

$\frac{f_m\left(\text{Δ}-{\text{Δ}}_s\right)}{1-k}=\text{ln}\frac{kJ}{\left(1-k\right)\left(1-J\right)},\left(\mathrm{C5}\right)$

or after inversion

$J\left(\tau \right)={\left[1+\frac{k}{1-k}{e}^{-\frac{f_m\left(\text{Δ}-{\text{Δ}}_s\right)}{1-k}}\right]}^{-1}\left(\mathrm{C6}\right)$

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

Here Eq. C2 with Eq. C3 yields

$f_m\text{Δ}\simeq \left(1-k\right){{\int }^{\text{}}}_{J_0}^{1-k}\frac{dy}{\left(1-y\right)y}+[J_{\infty }-\left(1-k\right)]{{\int }^{\text{}}}_{1-k}^J\frac{dy}{\left(1-y\right)\left(J_{\infty }-y\right)}\left(\mathrm{C7}\right)$$=f_m{\text{Δ}}_s+\left(J_{\infty }-1+k\right){{\int }^{\text{}}}_{1-k}^J\frac{dy}{\left(1-y\right)\left(J_{\infty }-y\right)},$

so that

$\frac{f_m\left(\text{Δ}-{\text{Δ}}_s\right)}{J_{\infty }-\left(1-k\right)}={{\int }^{\text{}}}_{1-k}^J\frac{dy}{\left(1-y\right)\left(J_{\infty }-y\right)}=-\frac{1}{1-J_{\infty }}\text{ln}\frac{1-\frac{1-J_{\infty }}{1-J}}{1-\frac{1-J_{\infty }}{k}}\left(\mathrm{C8}\right)$$=-\frac{1}{1-J_{\infty }}\text{ln}\frac{k\left(J_{\infty }-J\right)}{\left(J_{\infty }-1+k\right)\left(1-J\right)}$

After straightforward but tedious algebra we obtain

$\frac{J_{\infty }-J}{1-J}=\frac{J_{\infty }-\left(1-k\right)}{k}{e}^{-E},\left(\mathrm{C9}\right)$$E\left(\text{Δ}\right)=\frac{1-J_{\infty }}{J_{\infty }-\left(1-k\right)}{f}_m\left(\text{Δ}-{\text{Δ}}_s\right)\left(\mathrm{C10}\right)$

and consequently

$J\left(\text{Δ}\right)=1+\frac{J_{\infty }-1}{1-\frac{\left(1-k\right)\kappa }{k}{e}^{-E\left(\text{Δ}\right)}}\left(\mathrm{C11}\right)$

Using the identities $2{\left(1+e^{-2Y}\right)}^{-1}=1+\text{tanh}Y$ and $2{\left(1-e^{-2Y}\right)}^{-1}=1+\text{coth}Y$ we combine the results Eqs C6 and C11 to the analytical approximation of the SIR-model equations at all reduced times, stated in Eqs 10 – 12 above. A comparison with the exact numerical solution of the SIR model is provided in Figure C2 . The corresponding $j\left(\text{Δ}\right)$ is obtained from $J\left(\text{Δ}\right)$ via Eqs A5 .

FIGURE C2

Comparison for(A)$j\left(\text{Δ}\right)$ and(B)$J\left(\text{Δ}\right)$ 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=0$ Eq. A7 provides $J_0(k=0)=1/2$ so that with ${\text{lim}}_{k\to 0}{f}_m(k=0)=1$ the time scale ( Eq. C3 ) becomes

${\text{Δ}}_s(k=0)=\underset{k\to 0}{\text{lim}}\text{ln}\frac{1-k}{k}=\infty \left(\mathrm{D1}\right)$

With this result

$Y_1(k=0)=\frac{\text{Δ}}{2}-\frac{1}{2}\underset{k\to 0}{\text{lim}}\left[{\text{Δ}}_s+\text{ln}\frac{k}{1-k}\right]=\frac{\text{Δ}}{2}-\frac{1}{2}\underset{k\to 0}{\text{lim}}\left[\text{ln}\frac{1-k}{k}+\text{ln}\frac{k}{1-k}\right]=\frac{\text{Δ}}{2}\left(\mathrm{D2}\right)$

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

$J(\text{Δ},k=0)=\frac{1}{2}\left[1+\text{tanh}\frac{\text{Δ}}{2}\right],j(\text{Δ},k=0)=\frac{1}{4{\text{cosh}}^2\left(\text{Δ}/2\right)}\left(\mathrm{D3}\right)$

Appendix e: relationship between $J_0$ and K

Here we prove Eq. 14 . According to paper A the quantity $J_0$ is given by $J_0=1-e^{-G_0}$ with

$G_0=\frac{1+k}{k}+W_{-1}\left(\frac{-2e^{-1/k}}{ke}\right)=-\text{ln}\left(1-J_0\right)\left(\mathrm{E1}\right)$

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\left(c{e}^{-cr}/\beta \right)$ upon identifying $c=1$,$r=1/k$,$\beta =-ke/2$, and $x=-[1+\text{ln}\left(1-J_0\right)]$. From paper A we thus know that $e^{-cx}=\beta \left(x-r\right)$ holds, or equivalently

$\left(1-J_0\right)e=-\frac{ke}{2}\left[-1-\text{ln}\left(1-J_0\right)-\frac{1}{k}\right]=\frac{ke}{2}\left[1+\text{ln}\left(1-J_0\right)\right]+\frac{e}{2}\left(\mathrm{E2}\right)$

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

Appendix f: time of maximum in the measured differential rate $\dot{J}\left(t\right)$

One has $J\left(t\right)=J\left(\tau \left(t\right)\right)$ and $\dot{J}\left(t\right)=\dot{\tau }\left(t\right)J\text{′}\left(\tau \left(t\right)\right)=a\left(t\right)j\left(\tau \left(t\right)\right)$ since $j=J\text{′}$ if we let the prime denote a derivative with respect to τ. The maximum in $\dot{J}\left(t\right)$ thus fulfills

$0=\ddot{J}\left(t_m\right)=\dot{a}\left(t_m\right)j\left(\tau \left(t_m\right)\right)+a^2\left(t_m\right)j\text{′}\left(\tau \left(t_m\right)\right)\left(\mathrm{E1}\right)$

or equivalently,

$0={\left[\frac{d\text{ln}j}{d\tau }+\frac{\dot{a}}{a^2}\right]}_{t=t_m}\left(\mathrm{E2}\right)$

From part A we know that

$\frac{d\text{ln}j}{d\tau }=1-2J-k\text{ln}\left(1-J\right)-k\left(\mathrm{E3}\right)$

and our $J_0=J\left({\tau }_0\right)$ solves $1=2J_0+k\text{ln}\left(1-J_0\right)+k$. That is,$j\text{‘}\left({\tau }_0\right)=0$. If a does not depend on time,${\tau }_0=\tau \left(t_m\right)=a_0{t}_m$, but this is not generally the case. To find $t_m$ and ${\tau }_m\equiv \tau \left(t_m\right)$ one has to solve Eq. E1 , or Eq. E2 . Equation E2 with Eq. E3 is solved by

$J_m=J\left({\tau }_m\right)=1+\frac{k}{2}{W}_{-1}\left({\alpha }_m\right)\left(\mathrm{E4}\right)$

with

${\alpha }_m=-\frac{2e^{-\left(1+C_m\right)/k}}{ek},C_m=-\frac{\dot{a}\left(t_m\right)}{a^2\left(t_m\right)}\left(\mathrm{E5}\right)$

The corresponding j is, according to Eq. 4 ,

$j\left({\tau }_m\right)=\left(1-J_m\right)[J_m+k\text{ln}\left(1-J_m\right)]\left(\mathrm{E6}\right)$

The smaller $C_m$, the closer is $J_m$ to $J_0$.

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\left(t\right)$, from which we estimate the cumulative number of infections $J\left(t\right)=D\left(t-t_d\right)/f=200D\left(t-t_d\right)$ with $t_d=7$ days. From the cumulative value $J_m=J\left(t_m\right)$ at the time $t_m$ of the maximum in $\dot{J}\left(t\right)$ we estimate k via Eq. 14 upon assuming $J_m\approx J_0$. Similarly,$a_m=a\left(t_m\right)$ is estimated from $a_m=\dot{J}\left(t_m\right)/j_{\text{max}}\left(k\right)$. 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\left(t\right)$. The minimization is performed assuming the time-dependent $a\left(t\right)$ parameterized by Eq. H1 involving parameters $t_a>0$,$t_b>0$,$a_0>0$,$q\in \left[\mathrm{0,\; 1}\right]$,$q<\eta \in \left[\mathrm{0,\; 1}\right]$ and $t_s>t_m$. While $\tau \left(t\right)$ is given by the integrated $a\left(t\right)$, we use three strategies to model $J\left(t\right)$: i) the numerical solution of the SIR model, ii) the approximant $G\left(\tau \right)$ and $J\left(\tau \right)=1-e^{-G\left(\tau \right)}$ developed in part A, and iii) the approximant $J\left(\text{Δ}\right)$ given by Eq. 10 with $\text{Δ}=\tau -{\tau }_m$ and ${\tau }_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 $\text{Δ}$, 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\left(t\right)=a_{\text{LD}}\left(t\right)\text{Θ}\left(t_s-t\right)+a_{\text{stop}}\left(t\right)\text{Θ}\left(t-t_s\right)\left(\mathrm{H1}\right)$

where $t_s$ denotes the stop time of the lockdown still represented by the infection rate Eq. 15 , and where $a_{\text{LD}}$ is given by Eq. 15 . The infection rate after $t_s$ is assumed to be

$a_{\text{stop}}\left(t\right)=a_0\left[q_s+\left(\eta -q_s\right)\text{tanh}\frac{t-t_s}{t_b}\right]\simeq \{\begin{array}{ l}{q}_s{a}_0& \text{for}t=t_s\\ \eta a_0& \text{for}t\gg t_s,\end{array}\left(\mathrm{H2}\right)$

with $q_s=a_{\text{LD}}\left(t_s\right)/a_0$ the quarantine factor reached at the time $t_s$ of lifting. Together with the reduced time given by Eq. 16 we now find

$\tau \left(t\right)={\tau }_{\text{LD}}\left(t\right)\text{Θ}\left(t_s-t\right)+{\tau }_{\text{stop}}\left(t\right)\text{Θ}\left(t-t_s\right)\left(\mathrm{H3}\right)$

and

${\tau }_{\text{stop}}\left(t\right)={\tau }_{\text{LD}}\left(t_s\right)+q_s\left(t-t_s\right)a_0+\left(\eta -q_s\right)a_0{t}_b\text{ln}\left[\text{cosh}\left(\frac{t-t_s}{t_b}\right)\right]\left(\mathrm{H4}\right)$

with ${\tau }_{\text{LD}}\left(t\right)$ 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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