Methodology for modelling the new COVID-19 pandemic spread and implementation to European countries.

After the outbreak of the new COVID-19 disease, the mitigation stage has been reached in most of the countries in the world. During this stage, a more accurate data analysis of the daily reported cases and other parameters became possible for the European countries and has been performed in this work. Based on a proposed parametrization model appropriate for implementation to an epidemic in a large population, we focused on the disease spread and we studied the obtained curves, as well as, investigating probable correlations between the country's characteristics and the parameters of the parametrization. We have also developed a methodology for coupling our model to the SIR-based models determining the basic and the effective reproductive number referring to the parameter space. The obtained results and conclusions could be useful in the case of a recurrence of this insidious disease in the future.

Introduction

The new COVID-19 disease has been pandemic in the world during the first months of 2020. In most of the European countries the curve of daily new cases has reached a mitigation stage in May–June 2020 and it was time to turn our thoughts not only to its subsequent, painful and serious implications of this pandemic (Anastassopoulou et al., 2020; Peng et al., 2020; Li, 2020; Liu et al., 2020) but also to analyse mathematically the epidemiological data, that is, the new daily reported cases. In this work we performed this analysis in different European countries included UK and Switzerland for the time period from the beginning of pandemic spread since February 2020 until the 10 of May 2020, by using an appropriate mathematical function with three or more free parameters depending on the epidemic curve complexity. Together with the data parametrization (or fitting) we developed a methodology to determine the basic reproductive number which is a fundamental spread figure in the SIR epidemic models and its subsequent extensions. Additionally, we investigated the existent of correlation between the determined parameters of the parametrization and some general characteristics of the countries, like the financial level and the population density. The obtained results of this correlation study could give us information for preparing more effective defensive strategies and measures during a possible future return of the pandemic disease of COVID-19. The outline of the paper is as follows. In Section 2 we present a theoretical methodology for parametrizing an epidemic, whereas in Section 3 we explain how to couple the present parametrization model with the basic SIR model. In Section 4 we give results relevant to the end-to-end epidemics growth and in Section 5 we discuss the conclusions.

Theoretical methodology

Epidemic model

Our methodology is based on the parametrization of the growth of the COVID-19 disease that we used in our recent work (Maltezos, 2021) and also in (Vazquez, 2005; Ziff and Ziff, 2020), that we call “semi-gaussian of n-degree”. It was used for fitting the growth of the disease in various indicative countries and it belongs to the model category of “Regression Techniques” for epidemic surveillance. The basic-single term expression of this parametrization model iswhere the function c(t) applied in an epidemic spread represents the rate of the infected individuals as the new daily reported cases (DRC) and coincides with the function I(t) in the SIR model, described in Section 3. Also A is constant while n and τ are model parameters. The more analytical approach, in the general case from the mathematical point of view, comes from the fundamental study of the epidemic growth and includes a number of terms in a form of double summation related to the inverse Laplace Transform of a rational function given in (Kermack and McKendrick, 1927), referring to the “Earlier stages of an epidemic in a large population”. In this hypothesis, the number of unaffected individuals may be considered to be constant, while any alternation is assumed small compared to the total number of exposed individuals. This function, which could be called “Large Population Epidemic Semi-Gaussian model” (LPE-SG), is the followingwhere A are arbitrary amplitudes, n are the degrees of the model (assumed fractional in a general case) and τ are the time constants representing in our case a mean infection time respectively. Also, M and N are the finite number of terms to be included. It is easy to prove that the “peaking time” of the function of each term depends on the product of n by τ , that is, t p = n τ .

In practice, the number of the required terms should be determined according to the shape of the data and the desired achievable accuracy. After investigation of the fitting performance we concluded that, at most, two terms of the above double sum are adequate for our purpose. Also, the cross terms, with indices ij and ji, contribute even less to the accuracy of the model. In particular, a) the degree of the model can “cover” any early or late smaller outbreaks of the daily cases, while b) the mean infection time is a characteristic inherent parameter of the disease under study and thus should be essentially constant. For these reasons, the expression with one term was adequate in most of the cases, however, the use of 6 free parameters allow better flexibility for the fitting. Therefore, we can write

For the fitting procedure, we have used two alternative tactics, based on either the daily model or the cumulative integral of it. The decision depends on the “goodness of the fit” in each case which is based on the criterion of minimizing the χ 2/dof, as we have done in our previous work. (Maltezos, 2021). The starting date (at t = 0) is one day before the first reported case (or cases). The cumulative parametrization model and the fitting model take the following forms respectively where the symbol Γ represents the gamma function and the Γc the upper incomplete gamma function at time t.

Implementation strategy

The above generalized mathematical model has the advantage of providing more flexibility when the raw data include a composite structure or superposition of more than one growth curve which could be coexisting. This could possibly happen due to restriction measures applied during the evolution of an epidemic. Regardless of the number of terms used, the obtained parameters must be well understood, relative to their role in the problem. Let us consider a country where there is an outbreak of a disease due to a small number of infected individuals. In this stage we assume that the country has been isolated to a relatively high degree, given that a complete isolation is practically infeasible. At this point, the disease starts with a transmission rate which depends on the dynamics of the spread in each city and village, while other inherent properties of the disease affect its dynamics (e.g. the immune reaction, the incubation and recovering time).

At this point, we must also clarify the issue of the “size” of the epidemic. The SIR-based models assume that the size, N, that is the total number of individuals potentially exposed to the disease is unchanged during its evolution, a fact which cannot be exactly true. On the other hand, a fraction of the size concerns individuals who are in quarantine for different reasons (due to tracking or for precautionary reasons). Therefore, the size cannot be absolutely constant and the forecasting at the first stage (during the initial spread of the epidemic) would be very uncertain. In the second stage (around the turning point), the situation is clearer by means of more accurate parametrization, although high fluctuations could still be present. In the third stage (which is called mitigation or suppress), an overall parametrization is feasible and the obtained values of the free parameters can be used as initial ones in the parametrization of a future return of the epidemic. In any of the three stages regardless of the level of uncertainties, the parametrization specifies the associated parameters according to the epidemic model used. It is known that the “basic reproductive number” symbolized by R 0 is a very important parameter of the spreading of the epidemic. In the SIR-based models it is determined during the first moments of the epidemic (mathematically at t = 0) and is related to the associated parameters. Moreover, as it is described in Section 3, this parameter does not depend on the size N of the epidemic.

By using the present parametrization model we assume that the size of the epidemic is, not only unknown, but also much smaller than the population of the country or city under study, that is, it constitutes an unbiased sample of the potentially exposed generic population. Once the epidemic is almost completely eliminated, the size could be also estimated “a posteriori” by the help of a SIR-based models. However, in this case, the parameters of the spread, as well as, the reproductive number are already determined by the methodology given in Section 3. We consider that this more generic approach facilitates the fitting process and improves the accuracy because of the existence of an analytical mathematically optimal solution.

Coupling with the SIR-based models

Short description of SIR model

The classical model for studying the spread of an epidemic, SIR, belongs to the Mathematical or State-Space category of models along with a large number of other types which are briefly described in (Siettos and Russo, 2013), while some recent advanced models are given in (Shah et al., 2020; Abdo et al., 2020; Ahmad et al., 2020; Yousaf et al., 2020; Arfan et al., 2021; Rahman et al., 2020). Our model belongs to the continuum deterministic SIR models in the special case applied at the earlier stage of an epidemic, assuming that the population is much greater than the infected number of people (Kermack and McKendrick, 1927). This model can be applied also when the epidemic is in each later stage where the total number of infected individuals is an unbiased sample of the population. Under these assumptions, this model can be related to the classical SIR model, or even to its extensions (SEIR and SIRD), by means of correlating their parameters. Below, we give a brief description of the basic SIR epidemic model.

Let us describe briefly the three state equations of the SIR model:

The function S = S(t) represents the number of susceptible individuals, I = I(t) the number of infected individuals and R = R(t) the number of recovered individuals, all referred to per unit time, usually measured in [days]. The constant factor a is the transmission rate, the constant factor β is the recovering rate and N is the size of the system (the total number of individuals assumed constant in time), that is, S + I + R = N, at every time. The assumptions concerning the initial and final conditions are, S(0) ≠ 0, S(∞) > 0, I(∞) = 0 and R(0) = 0, R(∞) = N − S(∞).

This model does not have an analytical mathematical solution even assuming that the two parameters a and β which are constant during the spread of the epidemic. An approximate solution is derived only assuming, βR/N < 1, that is, when the epidemic essentially concerns a small number of recovered individuals compared to the total number ones. In this case a Taylor's expansion to an exponential function is used. In our study, we working with the general case without the above approximation.

Synergy in the parametric space

In our basic model, C(t), the integral of the function c(t), must be compared to the total number of infected individuals I, while the parameter A undertakes the scaling of the particular data. The parameter τ does not necessarily coincide with the inverse of the “ mean infection rate”, a, but 1/τ expresses an “effective transmission rate” in our model. The parameter n cannot be equalized to any of the parameters of the SIR model. However, this parameter contributes to the key parameter of the epidemic spread, R 0, which is defined at t = 0 and is equal to R 0 = (a/β)(S(0)/N), where β is the “mean recovering rate”. Because S(0) ≈ N, it becomes R 0 ≈ a/β. The fraction S(0)/N) represents a basic threshold, the so called “population density”, which depending on the value of a/β, would specify if R 0 ≤ 1 (the case of endemic, gradually decreasing) or R 0 > 1 (whereas an epidemic is initiated).

Moreover the “effective reproductive number”, R , as a function of time, is also defined in the same way as follows

Because the condition for creating an epidemic is R  > 1, the corresponding condition should be S/N > 1/R 0. Also, at t = 0 should be R (0) ≡ R 0, at the peaking time t = t p should be R (t p) = 1 and at t = ∞ takes the value . By using the expressions of Eq. (9), including only the value of I and its derivative, one can estimate roughly the R at any time t. It can be also shown, that the S(∞)/N can be determined by solving the following transcendental equation numerically.

Therefore, we can conclude that for R , at the outbreak of the epidemic (rising branch of the curve) we have, 1 < R  < R 0, exactly at the peak of the curve, R  = 1 (because S = N/R 0) at the mitigation stage (leading branch of the curve), R  < 1 tends to a minimum value at the asymptotic tail of the curve which is

Once the above relationship is achieved, the R 0 can be determined by solving the derived algebraic equation. Indeed, this was our initial motivation to perform the following analysis. The methodology for accomplishing it, was based on the idea of exploiting the property of our model at its maximum at the peaking time, which is t p = nτ, as it can be easily proven by differentiation. On the other hand, in the SIR model a peak is expected some time for the function I, as the typical effect of the epidemic's spread. Considering that both models can be fitted to the data and thus within the vicinity of the peak must agree, we can claim that I p = C(t p). Let us first find an expression of the S, R and I at the peaking time, symbolizing them by, S p, R p and I p respectively.

In order to relate S with R, we replace I from Eq. (6) into Eq. (8), we obtainfrom which, taking into account that S(0) ≈ N, we derive the solution

The above equation applied for the peaking time, t = t p, gives us an expression of R p

Also, according to Eq. (7), at the peaking time we might have

From this equation, and using the definition of R 0, we obtain

Based on Eq. (14) we calculate R p as follows

Adding the three functions at the peaking time, S p, I p and R p, we derive the algebraic equation

In order to form an equation independent of size N, we express I p/N as a function of the model parameters, that is, by using the maximum value of the model curve (Maltezos, 2021) and the Eq. (8) of the SIR model by integration setting as upper limit the infinity, obtainingwhere the symbol Γ represents the gamma function, n 0 and τ 0 are the particular values obtained by a fitting) and τ 0′ = 1/β. Replacing the above expression in Eq. (18) and setting s  = S(∞)/N we obtain

This transcendental equation can be solved only numerically for R 0 in which the combined unknown s N is also found numerically by using again another transcendental Eq. (10), by using multiple iterations leading to a converging accurate solution within 12 loops. The parameter n of the model is essentially the expresser of R 0, while the obtained value of R 0 concerns a hypothetical SIR model fitted to the data of the daily reported cases (DRC).

From the obtained solution for R 0 we can also calculate the parameter a of SIR model, a = βR 0, where β can be calculated from the peak value of the daily reported recovered individuals by the formula, β = R p/I tot, p, where I tot, p represents the integral of the DRC curve with upper limit the peaking time t p. In particular,

Implementing the above methodology, by using appropriate software codes written in Matlab platform (Mathworks, 2015), we obtained the fitting of the DRC curve for Greece at the mitigation stage, shown in Fig. 1 and Fig. 2 . The fitted parameters, n 0 = 4.57 and τ 0 = 5.96 and the solutions R 0 = 2.90 and s N = 0.06. For the parameter β we used a typical-average value found in the literature where, β = 0.10 and we used the same value for the analysis of the other countries. Two characteristic parametrizations of very large normalized size and very small one (64 times smaller), that is, of Belgium and Malta, respectively, are given in Fig. 3 and Fig. 4 . Obviously, without seeing the vertical scale, one cannot distinguish which corresponds to a large or to a small normalized size. The only visible difference at a glance, is the peaking time (29 and 16 days respectively).

Fig. 1

Indicative plot of the daily reported cases in Greece until May 9, 2020 (black circles) and the optimal curve of the fitted model (red solid line) obtained by the total reported cases.

Fig. 2

The corresponding total reported cases in Greece (black circles) with the fitting (blue solid line). The uncertainty zone concerns 1σ and is shown by the red dashed lines.

Fig. 3

Indicative parametrization curve (red solid line) of the daily reported cases (black circles) in Belgium, as a typical example of a very large normalized size (see Section 4).

Fig. 4

Indicative parametrization curve (red solid line) of the daily reported cases (black circles) in Malta, as a typical example of a small normalized size (see Section 4).

Study of the end-to-end epidemic growth

Correlation searches

For the data analysis we selected the 29 countries of EU, including Switzerland and UK obtained from (https://www.worldometers.info/coronavirus/, 2021). The characteristics of the countries relevant to our study are summarized in Table 1 . In particular, we used the population density, the estimated normalized total number of infected individuals (determined by the number of deaths by using a typical constant factor and per million of population) and the Gross Domestic Product (GDP), per capita. The degree of correlation among the above characteristics and the modelling parameters, was studied by the “theoretical pearson linear correlation coefficient” given by the following formulawhere X and Y are considered normal random variables, σ and σ are the corresponding standard deviations and Cov(X,  Y) is their covariance. However, as it is done in practice, we calculated the “sampling pearson coefficient” (SPC), r(X,  Y), for n observed random pairs (X ,  Y , …,  X ,  Y ), where the X can represent the first selected variable and Y the second one.

Table 1

Summary of the obtained parameters for 29 countries in Europe, after analysing all the available reported data. By we symbolize the epidemic size normalized to one million people. The epidemic curve in Bulgaria and UK had not yet reached the turning point at the time of the present study.

Country	D [p/km2]	N˜size/1Mp	GDP [EUR]	n0	τ0 [days]	R0	a = βR0 [days−1]
Austria	109	12,000	44,900	5.81	3.61	4.55	0.455
Belgium	383	115,200	41,200	4.77	6.13	2.78	0.277
Bulgaria	64	1600	8680	–	–	–	–
Croatia	73	2400	13,330	2.99	6.14	3.47	0.347
Cyprus	131	2400	24,920	2.86	5.99	3.65	0.365
Czechia	139	4000	20,640	4.21	5.57	3.23	0.323
Denmark	137	14,000	53,430	2.49	10.3	2.34	0.234
Estonia	31	7000	21,160	1.12	9.79	3.39	0.339
Finland	18	6400	43,480	2.18	12.5	2.14	0.214
France	119	60,200	36,060	6.49	5.22	2.80	0.280
Germany	240	13,800	41,350	6.95	4.83	2.92	0.291
Greece	81	2400	17,500	4.52	6.01	2.89	0.289
Hungary	107	5400	14,720	4.38	7.41	2.45	0.245
Ireland	72	41,000	59,910	4.38	5.09	3.49	0.349
Italy	206	86,000	29,610	4.31	8.37	2.26	0.226
Lithuania	43	3000	17,340	4.38	8.38	2.25	0.225
Luxembourg	242	27,200	102,200	2.46	5.41	4.59	0.459
Malta	1380	1800	26,350	2.04	7.90	3.20	0.32
Netherlands	508	50,000	46,820	4.54	7.28	2.46	0.246
Poland	124	2600	13,780	3.65	8.83	2.30	0.23
Portugal	111	16,800	20,660	2.55	7.65	2.98	0.298
Romania	84	6000	11,500	2.85	10.0	2.28	0.228
Serbia	100	3200	6590	6.90	4.67	3.03	0.303
Slovakia	114	600	17,270	3.21	7.84	2.64	0.264
Slovenia	103	7600	22,980	1.40	9.02	3.34	0.334
Spain	94	96,400	26,440	2.88	9.50	2.36	0.236
Sweden	25	42,600	46,180	2.93	12.9	1.95	0.195
Switzerland	219	36,800	69,760	3.79	5.34	3.57	0.357
United Kingdom	281	57,400	34,190	–	–	–	–

The correlation study took into consideration eight pairs, as is illustrated in Table 2 . The conclusions of the linear correlation study are the following:

Table 2

Summary of the correlation results.

Correlation pair	SPC	P-value	Statistically significant, C.I.
D−N˜	0.0645	0.749	No
GDP−N˜	0.311	0.113	No
n0 − τ0	−0.645	<0.001	Yes, 99%
D − R0	0.101	0.617	No
N˜−R_0	−0.193	0.335	No
tp − R0	−0.734	<0.001	Yes, 99%
t_p−N˜	0.309	0.117	No
tp − D	−0.169	0.399	No

For the population density D: no correlation was found with other parameters.

For the model parameters n 0 and τ 0: strong anti-correlation was found.

For the peaking time t p: very strong anti-correlation was found with the basic reproductive number, R 0.

The scatter plot of the basic reproductive number R 0 and the peaking time t p are shown in Fig. 6. This correlation gives us the following message: the higher R 0 results in less of a delay in the upcoming peak in the DRC curve. The obtained slope of the linear fit was −8.4 ± 0.2 days. On the other hand, the R 0 among the analyzed countries, present statistical fluctuations from about 2 to 4.6, obeying roughly a gaussian distribution with mean value 2.96 and standard deviation 0.68 (or relative to mean 23%). The parameters n and τ also fluctuate, as we can see in Fig. 7 while the peaking time t p shows stochastic characteristics obeying similarly a gaussian distribution with a mean 25.7 days and standard deviation 7.8 days (or relative to a mean 30%).

Fig. 6

A scatter plot of the determined R0 with the peaking time tp. The red solid line represents the linear fit.

Fig. 7

A scatter plot of of the determined n and τ. The red solid line represents the linear fit.

Since R 0 fluctuates (and in turn the time t p due to their linear correlation) among the different countries randomly without presenting any correlation with their associated parameters, we can conclude that the normalized size of the epidemic can be explained only by taking into account other reasons and aspects related to the way citizens interact and behave as well as the degree of social distances and mobility or transport within a country's major cities. Also, a crucial role was played by the degree of quarantine imposed and likely some individual biological differentiations (genetic and other related characteristics).

A plot of the obtained values of R 0 for the different European countries in descending order is presented in Fig. 5 . In this plot we see that R 0 in the different countries shows a variety of values between 1.95 and 4.59 without identifying any significant “step” in the histogram. The only find on which we could focused is the fact that the lowest value refers to Sweden which was the only country in which no strict lockdown was applied similar to that in the other European countries. However, the normalized size of the epidemic in this country was slightly smaller that the average value.

Fig. 5

Plot of the obtained values of basic reproductive number for the different European countries in descending order.

Quantitative surveillance of the epidemic

The capability to surveying the epidemic spread during the three main stages is very important and could be based on the daily data and the mathematical modelling we presented. In the mitigation stage the surveying is even more useful and crucial when the restrictive measures are starting to be relaxed. The crucial condition for a new epidemic reappearance is based on the effective reproductive number, R , as well as, on the corresponding population threshold. However, because of the large statistical fluctuations caused by the poor statistics of data as well as because of the low slope of the epidemic curve in this stage, it is very hard to achieve accurate numbers, but only a qualitative estimate as follows. The R can be estimated from the expressions in Eq. (9) and using average numerical approximations of the slope, dI/dt. An alternative and practical formula based on the parametrization model SG-LPE can be easily proven and is, R  = 1 + (nτ/t − 1)/β. This formula is valid only in the vicinity of the peak, namely in the narrow range from 0.5t to 1.5t , because the fitted model and the SIR one differ in the slopes at both side tails. Once R is estimated, the population density threshold, in turn, can be calculated and should be S(t)/N = 1/R , assuming that the normalized size N can also be estimated by a similar level of accuracy. Therefore the crucial condition in the mitigation stage is written as

The derivative has to be calculated as an average slope, ΔI/Δt, preferably at least within one week. Assuming that this slope is I w′ and the corresponding average cases in a week is I av the crucial condition becomeswhere we used the typical values, β = 0.1 days−1 and for having a practical result as a case study. This simplified formula combined with one-week measurements should be very useful because the relative fluctuations of the DRC are expected to be very large.

Conclusions

A systematic analysis of the epidemic characteristics of the new COVID-19 disease spread is presented in this work. For the mathematical analysis, we used a model that we called LPE-SG which facilitates the parametrization by an analytical mathematical description. We also presented a methodology of its coupling with the SIR-based models aiming to determine the basic and effective reproductive numbers based on the fitted parameters. Analysing the daily reported cases of European countries, we found no correlation between the population density, normalized size or GDP per capita of the countries with respect to the spreading characteristics. Another important finding of our study was a strong anti-correlation, statistically significant, of the basic reproductive number and the peaking time. Moreover, we found that the basic reproductive number in the epidemics studied showed a uniform distribution with a wide range of values. This means that it is mainly influenced by many other factors and genetic characteristics of the society in a country.

Declaration of competing interest

None.