key: cord-1017918-fev4bwag
authors: Köhler-Rieper, Felix; Röhl, Claudius H. F.; De Micheli, Enrico
title: A novel deterministic forecast model for COVID-19 epidemic based on a single ordinary integro-differential equation
date: 2020-07-25
journal: Eur Phys J Plus
DOI: 10.1140/epjp/s13360-020-00608-0
sha: 00e2738de3413c60f4bf1996fc6bc48fff83ea20
doc_id: 1017918
cord_uid: fev4bwag

In this paper, we present a new approach to deterministic modelling of COVID-19 epidemic. Our model dynamics is expressed by a single prognostic variable which satisfies an integro-differential equation. All unknown parameters are described with a single, time-dependent variable R(t). We show that our model has similarities to classic compartmental models, such as SIR, and that the variable R(t) can be interpreted as a generalized effective reproduction number. The advantages of our approach are the simplicity of having only one equation, the numerical stability due to an integral formulation and the reliability since the model is formulated in terms of the most trustable statistical data variable: the number of cumulative diagnosed positive cases of COVID-19. Once this dynamic variable is calculated, other non-dynamic variables, such as the number of heavy cases (hospital beds), the number of intensive-care cases (ICUs) and the fatalities, can be derived from it using a similarly stable, integral approach. The formulation with a single equation allows us to calculate from real data the values of the sample effective reproduction number, which can then be fitted. Extrapolated values of R(t) can be used in the model to make reliable forecasts, though under the assumption that measures for reducing infections are maintained. We have applied our model to more than 15 countries and the ongoing results are available on a web-based platform [1]. In this paper, we focus on the data for two exemplary countries, Italy and Germany, and show that the model is capable of reproducing the course of the epidemic in the past and forecasting its course for a period of four to five weeks with a reasonable numerical stability.

Our primary aim is to set up a deterministic model that can be easily tuned with available data in order to make numerically stable forecasts. We found that existing methods are not well suited to reach this goal.

Empirical top-down modelling, i.e., approaches that start from data and make prognoses, mostly ignores underlying dynamics. The easiest approach is curve fitting of available data. In [2] the number of cumulative diagnosed positive COVID-19 cases P(t) was assumed to be an error function. This is true if the number of daily new cases P (t) ( with respect to t) can be described by a Gaussian distribution. As we will show, a symmetric distribution function P (t) corresponds to an effective reproduction number that converges rapidly to zero. This might be true for China data. In Italy and Germany we observed a final (at the time of writing) value for R eff between 0.6 to 0.8, leading to an asymmetric function P (t) with a long tail. Although the peak date has been predicted well in [2] , the predicted total cases and fatalities differ by more than 30%. Current deterministic models were developed with the aim of simulating possible scenarios and showing the effect of containment and mitigation measurements. They are "bottom-up" in the sense that they are based on the knowledge of typical epidemiological parameters, such as the basic reproduction number R 0 or the time between contacts T c , just to name a few. However, three reasons make it difficult to set up these complex models for forecasting:

-The epidemiological parameters are unknown and change in time; -For most of the compartmental model variables, such as susceptible, exposed, infected or removed individuals, the availability of surveillance data is limited; -Model tuning requires fitting many variables simultaneously-making it difficult to find an optimum.

In [3] the classical SIR model has been applied to Italy, dividing the country into three parts: north, centre and south. The problem they face, in our opinion, is that the official number of infected individuals I contains people who are officially not cured. But in Italy people enter the statistics as cured when they have been tested as negative twice or even three times in a week's distance. Thus, from a dynamic point of view they remain "infected" for too long. The model cannot capture this feature appropriately and, in order to keep track of the statistical data, it has to be re-tuned within days.

The German Robert Koch-Institut (RKI) uses an extended SEIR model to show various scenarios for the course of the COVID-19 epidemic in Germany by applying different seasonality of the epidemic and immunity of the population [4] .

Another Italian team has set up a model with eight prognostic variables, SIDARTHE [5] , taking also into account asymptomatic cases and detection issues. Again, these efforts allow precise simulation of scenarios but are difficult to be set up with real data to make forecasts. The comparison with real data looks good but is restricted to the initial period of the epidemic when the case numbers grew simply exponential.

In [6] statistical parameters are obtained to feed parametric models, though not explicitly specified. Ensemble calculations using various data sources and different models allow for evaluating the statistical spread of the obtained forecasts-a procedure which is widely used in meteorological forecasting. The overall approach seems successful but remains complex.

Estimates of the disease transmissibility obtained through the evaluation of the timedependent reproduction number R t have been proposed by various authors (see, e.g., [7] [8] [9] . Wallinga and Teunis [7] proposed a statistical approach to compute an effective reproduction number which requires as input only the number of daily cases and the distribution of times intervals between the appearance of symptoms in primary cases and the onset of symptoms of secondary cases. The main drawback of this method is that in order to obtain estimates of R at time t, incidence data from times later than t are required.

To check the efficacy of restrictive measures adopted to contrast infectious disease, Bayesian estimation of the reproduction number R t along with Markov chain Monte Carlo and Monte Carlo sampling are employed in [8] to infer the temporal pattern of R t up to the last observation.

Cori et al. define in [9] the instantaneous reproduction number R t as the ratio of the new infections at time t to the total infectiousness of infected individuals at the same time. In this way, R t represents the average number of secondary cases that each infected individual would infect if the conditions remained as they were at time t. It is interesting to note that infectiousness of each individual is modulated by a weighting infectivity function which mainly depends on individual biological factors such as pathogen shedding (see also next Eq. (4)). Forecasting is then based on Bayesian statistical inference which leads to a simple expression of the posterior distribution of R t , assuming a prior gamma distribution for R t [9] .

If data-based forecasts are the primary scope, it seemed reasonable to us developing a hybrid approach: a simple dynamic model that can be easily tuned with available data. This goal is obtained with our approach based on a single prognostic variable, which is found to satisfy an ordinary integro-differential equation.

To our knowledge, there are only a few approaches that are equally simple and effective. In [10] a delay model is presented with a single prognostic equation that has even an analytic solution. Arguments and results are comparable to ours, though our integral formulation is more general and more robust when extracting parameters from available data to feed the prognostic model. Delay models [11, 12] can also be described with several variables and many parameters, which again makes them difficult to set up as forecasting model.

The paper is organized as follows. In Sect. 2 we derive the model and show that it can be interpreted as a generalisation of classical compartmental models, such as SIR. Section 3 is devoted to the analysis of real data from the COVID-19 epidemic in Italy and Germany. A summary of how the model is capable of handling data from other countries is given in Sect. 4. We conclude the presentation of our model with some remarks on stability and numerical robustness in Sect. 5. Finally, some conclusions are drawn in Sect. 6.

Many compartmental models, such as SIR, use deterministic equations for susceptible S, removed R and currently infected individuals I -all these variables being difficult to obtain from real data for various reasons. In our opinion, the most reliable statistic variable is the number of cumulative diagnosed positive cases. We choose this quantity as our model variable and denote it by P. We are aware of the fact that the diagnosed cases are only a part of all positive cases but we assume that:

(i) They are a statistically relevant part of the population; (ii) The fraction of diagnosed to all cases does not change through time and therefore, (iii) The dynamics applied to the visible part of the epidemic is representative for the entire epidemic.

Concerning point (ii), at the early stages of the epidemic the fraction of diagnosed cases obviously increased, also in response to the rapid increasing of the number of tests, reaching then an approximate stationary value. Variations of this stationary value are, however, negligible since it soon appeared that a large percentage (up to 80%, depending on the area) of people testing positive for COVID-19 may be asymptomatic. Therefore, symptom-based screenings, which are the most frequent epidemiologic investigations adopted in many countries, are likely to miss a lot of them (see [13] and references therein). Though the diagnosed cases are only a fraction of the total cases, what is important is that the proportion of asymptomatic cases be nearly constant over time (in this regard, see also the interesting simulation study reported in Web Appendix 8 of [9] ). One of the main objectives of an epidemiological analysis is to give estimates of the reproduction number R after restrictive measures (e.g., confinement, social distancing) have been adopted to limit epidemic spreading. With such constraints, assuming constant environment and exponential increase in new case counts appears unjustified [14] and a more empirical data-based approach is more appropriate to follow the temporal evolution of the reproducing number. We derive our model in a discrete version, using discrete daily values as they are given by various data sources. Successively, in Sect. 2.2 we illustrate its continuous version.

We refer to P n as the number of cumulative diagnosed positive cases on day n and to ΔP n as the number of newly infected COVID-19 cases on day n. We denote by R n the ratio between the new cases ΔP n on day n and the weighted sum of new cases on the previous N r days:

{g i } being a set of N r fixed weights with the property N r −1 i=0 g i = 1, where N r is the average number of days until an infectious person is removed from the infection process and ΔP n is a suitable average of ΔP n .

Let us first illustrate the main idea supporting our hybrid approach. The numbers R n can be calculated easily from the existing epidemic data. Then, a regression curve R(t) can be fitted to the set of numbers { R j } n j=N r , thus providing us with a law, purely based on data, on how this epidemic variable evolves. The data we observed show that, prior that restrictive measures have been adopted, R n has a nearly constant initial value corresponding to the basic reproduction number R 0 . Once contact behaviour changes (due to media information, measurements, quarantine, and so on) from a certain time T Q on, R n is no longer constant but manifests an evident decay towards a final asymptotic value. Therefore, we are prompted to choose the following model to describe the data R n :

α, R 0 and R ∞ being parameters to be determined from data R n . Note the difference between the discrete values R n , calculated from data by Eq. (1), which show sample fluctuations and the numbers R n ≡ R(n), which are the values of the regression curve evaluated at the sample day t ≡ n. Also note that phases of different severity in mitigation measures lead to small intermediate plateaus in the time course for R n . This behaviour of the data was also noted in [15] , especially for what concerns fatalities. Since the steps in R n are strongly smeared out, we simply model data by one single step with an asymptotic decay. Obviously, in order to gain higher accuracy, R(t) could be described with a more complicated, piece-wise defined function but at the price of introducing additional fitting parameters. The time T Q can be read rather easily from the course of the data and set approximately as the time when R n starts decreasing. Moreover, we will see in the next section that in the special case of constant weights g i = 1/N r , the quantity R(t) is actually seen to have the meaning of an effective reproduction number.

Resolved for ΔP n , Eq. (1) turns into a prognostic model for future infection cases ΔP n but provides also a model-based, smoothed curve for representing present and past data, i.e.:

(3)

Introducing weights g i is sensible because clinical data show that the probability to infect others is not equally distributed over time. The incubation time is known to be between 1 and 14 days, with an average of 5 days [4] . The infectiousness begins probably before symptoms manifest and is maximal at the beginning of the disease. All these characteristics can be captured with a suitable choice of the summation weights g i . The weights g i are samples of an infectivity probability (see, e.g., Fig. 2 ) over the time period of N r days in the past (i.e., t ≤ 0) that we have assumed to be Gamma-distributed with shape p and rate b. The corresponding probability density function then reads:

b and p being parameters to be fixed on the basis of the biological/clinical data of the disease. Note the definition of g(t) for negative values of t in order to provide probability values for the "past time" (see next Eqs. (5) and (6)). A Gamma distribution with suitable parameters describes well what we know about the temporal distribution of infectiousness of the disease: no or low infectiousness in the first few days, a rapid slope towards a maximum followed by a slow decay. Assuming that the infectiousness of a single individual is Gamma-distributed, the infectiousness of the sum of all individuals is again Gamma-distributed and therefore also the average infectiousness used in our deterministic model has this characteristic distribution. The Gamma distribution in the context of epidemic modelling was introduced in [16] for stochastic epidemic models and later applied in [17] for the derivation of quasi-stationary distributions of the SIS and SEIS model. The values of the parameters N r , b and p are given by the clinical observations and are typical of the disease. In Sect. 3, where we discuss the epidemiologic data, we used N r = 14, b = 0.75 and p = 4, which corresponds to the peak of infectivity after 4 days. Therefore, we have introduced a prognostic model with six degrees of freedom: three degree of freedom set by the clinical information, i.e., N r representing the infection or removal time in days, p and b for describing the infectivity probability; the other three degrees of freedom R 0 , α and R ∞ , obtained by fitting the regression curve R(t) in (2) to the data R n , represent the time dependance of the effective reproduction number upon the social restrictive measures. The logical steps associated with the proposed model are summarized in Fig. 1 .

The fitting parameters entering the model can be seen either as pure tuning numbers and also can be interpreted in epidemiological terms. In fact, the model as a whole can be compared to standard compartmental models, such as SIR, as we will show in the next section.

In order to compare our discrete model to classical continuous deterministic models, we write Eq. (3) in the following continuous form: where P(t) is the number of total COVID-19 cases, prime stands for derivative, R(t) represents, as we will show, an effective reproduction number, T r is the continuous generalization of N r given in Eq.

(1) and represents the time during which infected individuals take part in the infection process, and g(t) is a weighting function representing the infectivity probability:

Notice that g(t) is defined for negative values of t, emphasizing thus the fact that the averaging process described by the integral in (5) works over past times or, roughly speaking, that an individual found infected at time t (secondary case) has been actually infected some time before. Finally, recall that the numbers g i , which appear in (1) and (3), are samples of probability distribution g(t). It is worth observing that Eq. (5), read as the equation ruling R(t):

can be obtained as a particular case of the renewal equation for the birth process, where P (t) represents the observed birth rates and the time-varying infection rate R(t)g(−τ ) refers to the rate of production by a mother at the (positive) age τ [14, 18] .

The SIR model The SIR model, initiated by Kermack and McKendrick in 1927 to describe the plague spread mechanisms in Mumbai [19] , is a classic compartmental epidemic model that works with three prognostic time-dependent variables: the susceptible individuals S(t), the infected individuals I (t) and the people removed from the infection process r (t). There are transitions from S to I to r , which lead to the following system of ODE's [20] :

where N = S + I +r is the total number of population, β = 1/T c is the contact frequency, T c being the average time between contacts and 1/γ = T r is the mean time between infection and removal. We skip here the discussion of the SIR model (for further details, the interested reader is referred to [21] ) and ask: How does our model compare to SIR model? For this purpose, let us rewrite the SIR model in terms of our prognostic variable P(t) = I (t) + r (t).

If we add Eqs. (9) and (10), we obtain

which, in terms of P(t), becomes

Introducing relative susceptible number s = S/N , Eq. (12) can be written in the following form:

Now, let us rewrite our model (see Eq. (5)) by splitting the integral into two parts and replacing g(t) with constant weights g 0 = 1/T r . We have:

which, after integration, yields Comparing the SIR model in the form (13) with our model (15), we can identify

Recalling the standard definitions, γ = 1/T r , β = 1/T c and introducing the basic reproduction number R 0 = T r /T c = β/γ , we can state the following epidemiological interpretation of the parameters:

, the quantity R(t) assumes indeed the meaning of a time-dependent effective reproduction number:

(ii) The positive cases P(t − T r ) correspond to the removed individuals r (t); thus, T r can be consequently interpreted as the time until removal from the infection process.

In general, g(t) is not constant and it is typical of the disease. In this case, we can extend what discussed in point (i) above and give to R(t) the meaning of generalized effective reproduction number associated with the probability distribution g(t).

Delay models Delay models [22, 23] follow essentially the same strategy as SIR, the main difference being that the removal process is not modelled by a separate variable but with a time shift in the function describing the number of cumulative cases. In fact, Eq. (15) is identical to the functional retarded differential equation (11) in [10] , with the averaging weights g i set to the constant value g i = 1/T r .

The SEIR models The SEIR models (see [20] for a comprehensive review of the fundamental dynamics of these models) introduce a further group of people, the exposed E, that is, people who are infected but not yet infectious. This effect is accounted for in our model by excluding the first days into the integral in Eq. (5) or, equivalently, by using null weights, g i = 0, for those days. Though, from the clinical observation of the COVID-19 epidemic it seems that the probability that people are infectious already from the first days after infection [4] .

What makes the difference?

1. First of all, our approach does not explicitly model S(t) with a coupled prognostic equation. This sounds reasonable to us because the assumption that susceptible individuals are removed only by the infection process is wrong for the current COVID-19 epidemic. In fact, severe quarantine measures, including lockdowns and social distancing, have been implemented in almost all countries. 2. Other compartmental models take such measures into account by introducing, e.g., direct transfers from S to r compartments. However, in our opinion, this makes these models complicated and hides the fact that political measures and their effects are extremely difficult to model. Our approach is a very practical one: we do not model S(t). We focus on R(t), which we have seen being related to the product of the basic reproduction number R 0 and the time-dependent relative susceptible number s(t). We extract R(t) from real data using our model assumptions and apply a curve fitting procedure to allow for extrapolation. Therefore, our approach can be called "hybrid": a mixture of curve fitting and modelling. 3. By using the number of cumulative diagnosed positive cases P(t) as prognostic variable, we automatically have the numbers of new infections as ΔP n . In our opinion, this is the best variable to describe how the epidemic evolves. In SIR models, this number is not automatically obtained since I represents the "currently infected people" and ΔI is a net difference mixing the "new positive cases", that is, the transfer from S to I , with the "removed cases", i.e., the transfer from I to r . 4. By introducing the weights g i we can model the incubation time, namely, the time between "getting infected" and "being infectious", as well as the time before detection. 5. We define the removed people r (t), which counts the individuals that no longer take part in the infection process, as all positive cases P(t − T r ) at a certain previous time T r . Again, we have no prognostic equation for r (t). Therefore, we do not have to model and elaborate how these individuals are removed from the infection process. We simply assume they are removed after a certain time T r because, prior to curing or deceasing, people are isolated in hospitals or, in the case of weak symptoms, are put into quarantine. 6. Parameters or variables that are not known, such as T c , R 0 , S(0) and N , are subsumed into a single function R(t), which is obtained from real data, without having to speculate on how it comes about. This makes the model simple-and setting the weights constanteven simpler so that it can be set up quickly to produce satisfying prognoses. 7. From the comparison with the compartmental SIR models, we see that our epidemiological variable R(t) assume the meaning of generalized effective reproduction number associated with the infectivity probability g(t), which is typical of the disease.

With these assumptions we have been able to describe the infection process with a single prognostic variable P(t) in an integro-differential equation. From our perspective, the computation of deceased and cured people is a secondary process, which does not influence the dynamics of the epidemic. Nonetheless, they are important numbers to know and, however, they can be simply obtained from P(t), as it will be shown in the next section.

The analysis of the "number of fatalities", "number of cured" and "number of active cases" follows mutatis mutandis the analysis presented in Sect. 2.

The number of fatalities Let V n be the total number of deceased individuals on day n from the beginning of the epidemic. We assume that the casualties on the nth day, i.e. ΔV n , is related to the weighted sum of new cases over the last N V days, yielding the definition of the following empirical ratio:

where N V is the maximum number of days after which the people decease and the weights h i allow for taking into account a probability distribution of deceasing. The numbers μ n can be interpreted as case fatality ratios [24] . Similarly to what we have done in Sect. 2, the values μ n can be computed from existing data, and then fitted to a model function μ(t) to extrapolate future values: μ n ≡ μ(t n ), where t n is a day in the future. The corresponding discrete prognostic equation reads

Note that from the analysis of currently available data, it results that the peak of fatalities lacks about 7-8 days behind the peak of the daily new infections. This has been accounted for in the weights h i , having set a maximum value at 7-8 days prior to the current day t n .

The number of cured Let C n denote the number of currently cured individual and ΔC n the number of new cured individuals at day n. The empirical discrete curing ratio ν n can be defined analogously (see Eqs. 1 and 18):

where, as in the previous cases, w i denote suitable weights. Even in this case, the discrete curing ratio can be fitted with a suitable model function ν(t) to obtain predictions: ν n ≡ ν(t n ), which can be used to predict future number of cured people through the relation

Active cases The number of currently infected individuals I (t), also known as active cases, is the difference between all cases P(t) and the deceased and cured cases:

Note that calculating the number of removed cases in the usual way is not valid for our model, that is, we have:

since our removal process is obtained by cutting off the corresponding integral after a removal time T r , taking thus into account not only the usual removal processes due to curing and deceasing, but including even other processes such as quarantine or isolation. However, we have to guarantee by means of the suitable choice of our tuning parameters that, in the long term,

since, in this limit, I (t) = 0. 

On the left of Fig. 2 we see the Γ −distributed weights g i we used to obtain the effective reproduction numbers R n (see Eq.

(1) and the figure caption for numerical details) shown on the right of the same figure. The integration time is N r = 14 days, i.e. only individuals registered positive within this time period actually take part in the (model) infection process. The effective reproduction numbers R n are shown for the time period from March, 3 to April, 14. The data-based R n are shown with blue dots and the model-fitting curve with the orange line. The data shows two short plateaus. The first one at R 0 ≈ 3.3 represents the The original time series data show significant weekly fluctuations; therefore, we show only the 7-day average. The model is capable of reproducing the exponential growth in the beginning of the epidemic, as well as the peak and the slow decay of the curve afterwards. The deviation of the curve ΔP(t) from the data Δ P n is the consequence of the deviation between R n and fitting curve R(t). The cumulative number of diagnosed positive cases P(t) is shown on the right of Fig. 3 . The model curve P(t) follows rather accurately the data P n . Note that the Daily new fatalities ΔV (t) (left) and cumulative number of fatalities V (t) (right) from February 26 to July 1, 2020 in Italy. Only data before the vertical black line, i.e., before April 13, have been used to tune the model last model-tuning was made on April 13, 2020. Two weeks later the relative deviation of the cumulative number of cases is about 2%. Approximately 2.5 months later, at the moment of revision, using the tuning of April, 13 we had a deviation of less than 10%.

On the left of Fig. 4 we see the Gaussian weights h i (see Eq. (18)) used to obtain the model case fatality rate μ n shown on the right of the same figure. The integration time is N r = 18 days, i.e., only individuals diagnosed positive within the last 18 days are considered in the model calculation of fatalities. Note that setting the Gaussian time-shift, that is the location of the Gaussian peak, at t shift = 6 puts the maximum weight on patients that have been diagnosed positive 6 days before. This is sensible because the peak of daily fatalities occurs 6 days after the peak of daily new infections. The number of daily ΔV (t) and cumulative fatalities V (t) is shown in Fig. 5 . The course of fatalities is well represented 

On the left of Fig. 6 we see the Γ −distributed weights g i we used to obtain the effective reproduction numbers R n shown on the right panel. Again, the integration time is N r = 11 days.

The effective reproduction number R(t) is shown for the time period from March, 3 to April, 14. The most restrictive measure in Germany was the school closing on March, 14. Note that some days before the R n show a short plateau at R 0 ≈ 4.0, which can be interpreted as the initial basic reproduction number R 0 of COVID-19 in Germany. After three weeks, R n tends to the asymptotic value of about R ∞ ≈ 0.6. We modelled this behaviour from day March 3, 2020 according to Eq. (1), the parameters being: R 0 = 4.10, α = 0.12, R ∞ = 0.64. If we compare this to the analysis of Italy, we note that 1. The initial effective reproduction number R 0 was higher in Germany; 2. Its final value R ∞ is higher in Italy; 3. The rate parameter α is in both countries approximately the same.

A tentative interpretation can be the following: The initial reproduction number in Germany was higher because at the beginning of the epidemic the disease spread mainly among young people coming from skiing resorts in the Austrian Alps. If we assume that social contacts among young and sporty people are more frequent, this could be an explanation.

Although the measures taken by Italian politicians were more restrictive than in Germany, the effective reproduction number R ∞ at mid April was 30% higher in Italy, leading to a much slower decay of daily new positive cases ΔP(t).

The COVID-19 cases in Germany The number of daily positive diagnosed cases is shown on the left of Fig. 7 . Original data are shown as a 7-day averages (blue dots). The model function (orange line) is capable of reproducing the course of the epidemic correctly but the peak value at March, 30 is underestimated. On the right panel the model function for Fatalities in Germany On the left of Fig. 8 we see the Gaussian weights we used to obtain the case fatality rates μ n shown on the right panel. The integration time is N r = 18 days. Note a time shift of the Gaussian by t shift = 10 days putting a maximum weight on patients that have been diagnosed positive 10 days earlier, which corresponds to the average time in hospital before deceasing given by the Robert-Koch Institut [4] .

The number of daily and cumulative fatalities for Germany is shown in Fig. 9 . The course of fatalities was not well represented due to a sudden rise in case fatality rate about April, 15. This effect is also visible on the right panel of Fig. 8 , where there is a sudden increase Other variables Once the function describing the cumulative positive diagnosed cases P(t) is known, other variables can be derived, such as the number of: -Cured individuals, -Hospital beds needed, -Intensive care units (ICU's) needed, -Active and closed cases, and so forth. Forecasting these variable is made with an analogous weighted-integralapproach as for the fatalities, see Sect. 2.3. In Fig. 10 we show an example for the number of cumulative cured and active cases for Germany. Note that the same diagram for Italy (not shown) was difficult to obtain because the data of cured cases appeared not to be reliable. In fact, cured individuals were registered with a very long time delay (as to date of writing, cured cases made up only 65% of all closed cases).

We have applied our model to the COVID-19 data made available for other countries by the European Centre for Disease Prevention and Control [27] . The corresponding graphs are available on a web platform [1]. Here, we simply summarise the fitting parameters in Table 1 . The date of fitting was April 28, 2020.

The following facts can be observed:

-The lowest final reproduction number is 0.5.

-Some countries, such as Brazil and Sweden, still had R(t) > 1.

-The highest slope can be seen with South Korea, that has quickly introduced severe measures. 

Concerning the numerical error propagation for our model, we expect a high stability due to its integral formulation. This is important when applying the model to noisy time series data because the parameters that govern the overall dynamic are not obtained from a single data point but from a weighted sum of data points, see Eq.

(1). We cannot make a theoretical stability analysis and therefore we limit ourselves to show empirically how P(t) reacts to small changes in the curve fitting parameters. To this purpose, we individually disturb each of the parameters R 0 , α and R ∞ associated with the fitting curve R(t) in Eq. (2) by 5% and observe the resulting change δ P in P(t) after 2, 4, 6 and 8 weeks. The results are given in Table 2 . Table 2 Model reaction to a 5%-perturbation of the fitting parameters δ P(2 weeks) δ P (4 weeks) δ P (6 weeks) δ P (8 weeks) Table 3 Model reaction to a 10%-perturbation of the fitting parameters δ P(2 weeks) δ P (4 weeks) δ P (6 weeks) δ P (8 weeks)

Model deviations remain reasonably limited for mid-term forecasts up to two months. If we double the perturbation to 10% we observe changes in the cumulative number of diagnosed positive cases P(t) as shown in Table 3 .

It can be noted that doubling the perturbation roughly leads to the doubling of the deviation of the function. Thus, the model can be seen as numerically robust. The same results are expected for secondary variables, such as fatalities because they depend linearly on P(t).

The main advantage of our approach is the simplicity of its formulation, the precision with which the real course of significant variables can be reproduced and the effectiveness to make mid-term forecasts.

We were able to show that our model has similarities to classic compartmental models, in particular:

-The variable R(t) can be interpreted as an effective reproduction number.

-The limited integration interval in the deterministic equation models the removal process.

-The Gamma−distributed weights account for infectiousness within a latency period, corresponding to the exposed state of the SEIR model.

Since our model consists of only one deterministic equation it is simpler compared to most approaches but is nevertheless able to capture the time course of the epidemic. In addition, the integral formulation leads to a good numerical robustness. The model contains six parameters: three parameters related to infectiousness, N r , b and p, can be set from biological/clinical data typical of the disease; three free parameters, α, R 0 and R ∞ that can be obtained by the fitting procedure of the set of sample dataR n . In the case of extinction of the epidemic outbreak, the parameter R ∞ is expected asymptotically to vanish. A non-null value of R ∞ indicates that the epidemics remain "latent", with a relative small number of daily positive diagnosed cases for a very long time. This is indeed what seems to happen in various countries where the primary outbreak has been contained by adopting (even severe) social restrictive measures. Hence, deviations from a reliable forecasted value of R ∞ could be interpreted as a warning signal of oncoming second epidemic wave [28] . We applied our model to many countries with a special focus on the data of Italy and Germany. After extracting the three fitting parameters we have been able to model the course of the epidemic in both countries rather well. We found it interesting how the parameters R 0 and R ∞ differ between the two countries-inviting us to interpret them in terms of effectiveness of measures, social organisation (in Italy elderly vulnerable people are more likely found to live with the younger part of the family) as well as the organisation of the health system.

We set up the hypothesis that three parameters suffice to model the epidemic from the outbreak, over the period of social distancing measures until the end-under the assumption that the measures remain effective with respect to infections till the end, i.e., zero new infections. It remains to be shown that this hypothesis remains valid for longer periods of time, especially when mitigation measures are loosened.

We hope that our approach facilitates forecasting of ongoing epidemics for long term periods, providing early warnings of further epidemic waves. Updating the model and retuning is very easy and we collect the results on the website [1].

A Mathematical prediction of the time evolution of the Covid-19 pandemic in some countries of the European Union using Monte Carlo simulations

Modellierung von Beispielszenarien der SARS-CoV-2-Epidemie 2020 in Deutschland

A SIDARTHE Model of COVID-19 Epidemic in Italy

Modeling projections for COVID-19 pandemic by combining epidemiological, statistical, and neural network approaches

Different epidemic curves for severe acute respiratory syndrome reveal similar impacts of control measures

Estimating in real time the efficacy of measures to control emerging communicable diseases

caucheme A new framework and software to estimate time-varying reproduction numbers during epidemics

Solvable delay model for epidemic spreading: the case of Covid-19 in Italy

Global stability analysis for a generalized delayed SIR model with vaccination and treatment

Dynamic models for Coronavirus Disease 2019 and data analysis

COVID-19: What proportion are asymptomatic?

How generation intervals shape the relationship between growth rates and reproductive numbers

Estimating the number of infections and the impact of nonpharmaceutical interventions on COVID-19 in 11 European countries

Stochastic epidemics in dynamic populations: quasi-stationarity and extinction

An application of queuing theory to SIS and SEIS epidemic models

On the integral equation of renewal theory

A contribution to the mathematical theory of epidemics

The mathematics of infectious diseases

Dynamic behaviors of a modified SIR model in epidemic diseases using nonlinear incidence and recovery rates

Some epidemiological models with delays

Time Delays in Epidemic Models

Epidemiological Research: Terms and Concepts

The end of social confinement and COVID-19 re-emergence risk

Acknowledgements We wish to thank Prof. Klaus Dietz for a very useful and illuminating discussion.