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

Abstract

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

Let , and denote the infected, susceptible and recovered/removed fractions of persons involved in the infection at time t, with the sum requirement . In terms of the reduced time , accounting for arbitrary but given time-dependent infection rates, the SIR-model equations are [1, 2, 18]in terms of the time-dependent ratio of the recovery and infection rates and the medically interesting daily rate of new infectionswhere the dot denotes a derivative with respect to t.

For the special and important case of a time-independent ratio 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 . The new analytical solutions assume that the SIR equations are valid for all times , and that time refers to the “observing time” when the existence of a pandemic wave in the society is realized and the monitoring of newly infected persons is started. In paper A it has been shown that, for arbitrary but given infection rates , apart from the peak reduced time 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 is controlled by k and the value , indicating as only initial condition at the observing time the fraction of initially susceptible persons . This suggests to introduce the relative reduced time with respect to the reduced peak time. In real time t the adopted infection rate acts as second parameter, and the peak time , where reaches its maximum must not coincide with the time, where the reduced j reaches its maximum, i.e., , in general.

2 Results and Discussion

According to paper A the three fractions of the SIR-modelcan be expressed in terms of the cumulative number and differential daily rate of new infections. The cumulative number satisfies the nonlinear differential equationTwo important values are , where j attains its maximum with , and the final cumulative number at , when the second bracket on the right-hand side of the differential Eq. 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 (Figure 1A)Without any detailed solution of the SIR-model equations the formal structure of Eqs 3 and 4 then provides the final values , , and . We list these values together with κ in Table 1. We emphasize that the final cumulative number , 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

2.1 New infections

The exact solution of the differential Eq. 4 is given in inverse form by (Appendix A)which can be integrated numerically (subject to numerical precision issues), replaced by the approximant presented in paper A (involving Lambert’s function), or semi-quantitatively captured by the simple approximant to be presented next. The solution as a function of the relative reduced time , with the reduced peak time approximated bycorresponding to in Eq. 8, and where , is reasonably well captured by (Appendix C)with the Heaviside step function for . In Eq. 10withalso tabulated in Table 1. We note that 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 is captured exactly by the approximant (Appendix D).

FIGURE 2

For the corresponding reduced differential rate in reduced time we use the right hand side of Eq. 4 with from Eq. 10, cf. Figure 3. Note, that this j is not identical with the one obtained via , because J does not solve the SIR equations exactly. The peak value in the reduced time rate occurs when and is thus determined by , also tabulated in Table 1.

FIGURE 3

As can be seen in Figure 3 the rate of new infections (Eq. 12) is strictly monoexponentially increasing with well before the peak time, and strictly monoexponentially decreasing well above the peak time with the . These exponential rates exhibit a noteworthy property and correlation in reduced time:The SIR parameter k affects several key properties of the differential and cumulative fractions of infected persons. If the maximum of the measured daily number of newly infected persons has passed already, we find it most convenient to estimate k from the cumulative value at this time . While the maximum of must not occur exactly at (Appendix F), we can still use as an approximant for the value of and the relationship between and k can be inverted to read (Appendix E)The dependency of k on is shown in Figure 1C. With the so-obtained value for k at hand, the infection rate at peak time can be inferred from . It provides a lower bound for .

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 . Such time-dependencies result at times greater than the observing time 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 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 and cumulative number of new infections in real time t for given time-dependent infection rates are and given by Eq. 2. From a medical point of view the daily rate is most important as it determines also i) the fatality rate [23] with the fatality percentage in countries with optimal medical services and hospital capacities and the delay time of  days, ii) the daily number of new seriously sick persons [24] needing access to breathing apparati, and iii) the day of maximum rush to hospitals . 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 physicswhich impliesThe time-dependent lockdown infection rate (Eq. 15) is characterized by four parameters: i) the initial constant infection rate at early times , ii) the final constant infection rate at late times described by the quarantine factor , first introduced in Refs. 21 and 24, iii) the time of maximum change, and iv) the time regularizing the sharpness of the transition. The latter is known to be about –14 days reflecting the typical 1–2 weeks incubation delay. Consequently, the parameter q mainly affects the amplitude shown in the left columns of Figures 5 and 6 (note that we also plotted the case of no NPIs taken (i.e., ) for comparison). Alternatively, the transition time controls the rapidness of the transition in the fraction of infected persons per day and therefore the widespread.

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

FIGURE 4

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

In Germany the first two deaths were reported on March 9 so that corresponding to about 400 infected people 7 days earlier, on March 2 (). The maximum rate of newly infected fraction, , occurred  days later, consistent with a peak time of fatalities on 16 April 2020. At peak time the cumulative death number was corresponding to . This implies according to Eq. 14 (and not far from the value to be determined from the fit shown in Figure 5). From the initial exponential increase of daily fatalities in Germany we extract , corresponding to a doubling time of  days, as we know already from the above k. The quantity we can estimate from the measured , as and . Using the mentioned value for , we obtain days as a lower bound for .

FIGURE 5

With these parameter values the entire following temporal evolution of the pandemic wave in Germany can be predicted as function of 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 ,  days, ,  days, and  days. The obtained parameters allow us to calculate the dimensionless peak time , the dimensionless time , as well as , , and .

We note that the value of implies for Germany that 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 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 is the implied flat exponential decay after the peak. Because for , the exponential decay is by a factor q smaller than the exponential rise prior the climax, i.e.,Equation 17 yields a decay half-live of  days  days to be compared with the initial doubling time of  days. For Germany we thus know that their lockdown was drastic and rapid: the time March 23 is early compared to the peak time April 8 resulting in a significant decrease of the infection rate with the quarantine factor . 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 and 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 for all other countries investigated corresponding to basic reproduction numbers . These values are significantly smaller than the estimates of in the mainstream literature on Covid-19 [4, 22]. Part of these significant differences may be explained by the different definitions of .

FIGURE 6

While the inverse basic reproduction number 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 using the effective reproduction factor . As discussed in detail in Sect. 4 of Ref. 25 has to be calculated from the convolutionof the number of daily cases with the serial interval distribution 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 . As is identical to the value at the starting time of the outbreak it is not clear in the moment that this will be identical to the 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 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 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 , 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 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 from its present value up to a value that might be close to the initial . 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 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 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

References

• 1.

  KermackWOMcKendrickAG. A contribution to the mathematical theory of epidemics. Proc R Soc A (1927). 115:700–21.

  • Google Scholar

• 2.

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

  • Google Scholar

• 3.

  HethcodeHW. The mathematics of infectious diseases. SIAM Rev (2000) 42:599–653. 10.1137/s0036144500371907

  • CrossRef
  • Google Scholar

• 4.

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

  • CrossRef
  • Google Scholar

• 5.

  O’ReganSMKellyTCKorobeinikovAO’CallaghanMJAPokrovskiiAV. Lyapunov functions for SIR and SIRS epidemic models. Appl Math Lett (2010) 23:446–8. 10.1016/j.aml.2009.11.014

  • CrossRef
  • Google Scholar

• 6.

  SatsumaJWilloxRRamaniAGrammaticosBCarsteaAS. Extending the SIR epidemic model. Phys Stat Mech Appl (2004) 336:369–75. 10.1016/j.physa.2003.12.035

  • CrossRef
  • Google Scholar

• 7.

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

  • CrossRef
  • Google Scholar

• 8.

  CadoniMGaetaG. Size and timescale of epidemics in the sir framework. Phys D (2020) 411:132626. 10.1016/j.physd.2020.132626

  • CrossRef
  • Google Scholar

• 9.

  ChekrounAKuniyaT. 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. 10.1016/j.jde.2020.04.046

  • CrossRef
  • Google Scholar

• 10.

  ImronCHariyantoYunusMSurjantoSDDewiNAC. Stability and persistence analysis on the epidemic model multi-region multi-patches. J Phys Conf Ser (2019). 1218: 012035. 10.1088/1742-6596/1218/1/012035

  • CrossRef
  • Google Scholar

• 11.

  KarajiPTNyamoradiN. Analysis of a fractional SIR model with general incidence function. Appl Math Lett (2020). 108:106499. 10.1016/j.aml.2020.106499

  • CrossRef
  • Google Scholar

• 12.

  MohamadouYHalidouAKapenPT. 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. 10.1007/s10489-020-01770-9

  • CrossRef
  • Google Scholar

• 13.

  SamantaSSahooBDasB. Dynamics of an epidemic system with prey herd behavior and alternative resource to predator. J Phys Math Theor (2019). 52:425601. 10.1088/1751-8121/ab264d

  • CrossRef
  • Google Scholar

• 14.

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

  • CrossRef
  • Google Scholar

• 15.

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

  • CrossRef
  • Google Scholar

• 16.

  TianCZhangQZhangL. Global stability in a networked SIR epidemic model. Appl Math Lett (2020). 107:106444. 10.1016/j.aml.2020.106444

  • CrossRef
  • Google Scholar

• 17.

  DehningJZierenbergJSpitznerFPWibralMNetoJPWilcekMet alInferring change points in the spread of COVID-19 reveals the effectiveness of interventions. Science (2020). 369:eabb9789. 10.1126/science.abb9789

  • CrossRef
  • Google Scholar

• 18.

  BarmparisGDTsironisGP. Estimating the infection horizon of covid-19 in eight countries with a data-driven approach. Chaos Solitons Fractals (2020) 135:109842. 10.1016/j.chaos.2020.109842

  • CrossRef
  • Google Scholar

• 19.

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

  • CrossRef
  • Google Scholar

• 20.

  FergusonNLaydonDNedjati-GilaniG. Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand. London: Imperial College London (2020) 10.25561/77482

  • CrossRef
  • Google Scholar

• 21.

  SchlickeiserRKrögerM. Dark numbers and herd immunity of the first covid-19 wave and future social interventions. Epidem Int J (2020). 4:000152. 10.23880/eij-16000152

  • CrossRef
  • Google Scholar

• 22.

  FlaxmanSMishraSMishraSGandyAUnwinHJTMellanTAet alEstimating the effects of non-pharmaceutical interventions on Covid-19 in Europe. Nature (2020) 584:257–261. 10.1038/s41586-020-2405-7

  • CrossRef
  • Google Scholar

• 23.

  SchüttlerJSchlickeiserRSchlickeiserFKrögerM. Covid-19 predictions using a gauss model, based on data from april 2. Physics (2020). 2:197–212. 10.3390/physics2020013

  • CrossRef
  • Google Scholar

• 24.

  SchlickeiserRSchlickeiserF. 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. 10.3390/physics2020010

  • CrossRef
  • Google Scholar

• 25.

  KrögerMSchlickeiserR. Gaussian doubling times and reproduction factors of the Covid-19 pandemic disease. Front Phys (2020). 8:276. 10.3389/fphy.2020.00276

  • CrossRef
  • Google Scholar

• 26.

  SchlickeiserRKrögerM. 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. 10.3929/ethz-b-000456421

  • CrossRef
  • Google Scholar