Systematic description of COVID-19 pandemic using exact SIR solutions and Gumbel distributions

Nonlinear Dyn https://doi.org/10.1007/s11071-022-07907-4 ORIGINAL PAPER Systematic description of COVID-19 pandemic using exact SIR solutions and Gumbel distributions J. E. Amaro Received: 31 May 2022 / Accepted: 13 September 2022 © The Author(s) 2022 Abstract An epidemiological study is carried out in several countries analyzing the first wave of the COVID-19 pandemic using the SIR model and Gumbel distribution. The equations of the SIR model are solved exactly using the proper time as a parameter. The physical time is obtained by integration of the inverse of the infected function over proper time. Some properties of the solutions of the SIR model are studied such as time scaling and the asymmetry, which allows to obtain the basic reproduction number from the data. Approximations to the solutions of the SIR model are studied using Gumbel distributions by least squares fit or by adjusting the maximum of the infected function. Finally, the parameters of the SIR model and the Gumbel function are extracted from the death data and compared for the different countries. It is found that ten of the selected countries are very well described by the solutions of the SIR model, with a basic reproduction number between 3 and 8. Keywords COVID-19 coronavirus ·SIR model · Differential equations ·Gumbel distribution J. E. Amaro (B) Departamento de Física Atómica, Molecular y Nuclear and Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada, 18071 Granada, Spain e-mail: [email protected] 1 Introduction Since the declaration of the COVID-19 pandemic by the World Health Organization in March 2020, studies by the mathematical epidemiologist community have intensified and models of various kinds have been developed in order to provide insights and make predictions about the spread of the disease [1–4]. Epidemiological models have been used as basic tools in epidemiology for over a century [5–13] and have been extensively used prior to the COVID-19 pandemic [14– 17]. A wide range of models have been proposed and tested in an attempt to describe the COVID-19 data and forecast the future evolution of the pandemic in different regions of the planet. From very simple models [18,19] to numerous variants of compartmental models based on the SIR model (susceptible, infected, removed) have been proposed [20–23]. Additional models have been employed such as the SEIR model [24], which adds the exposed compartment of individuals, the uncertain SIR model [25], and others that include different parameters that statistically describe the many factors that may influence the pandemic dynamics. It is worth mentioning the SITR model [26], which also includes the treatment process T, and the SITRS model [27], which includes the possibility that recovered people can lose their temporary immunity against the virus and subsequently join the susceptible compartment again. These models may also include two or 123 J. E. Amaro more susceptible groups of individuals S1,S2,... to take into account the different susceptibilities due to age or other factors. Disease simulations are obtained by solving a system of nonlinear differential equations, or by means of discretization methods of different types, such as the discrete fractional model [27]. More sophisticated approaches also take into account the spatial spread of the disease, not just the temporal one. Models may involve partial differential equations in space and time coordinates, or, alternatively, stochastic methods [28–31]. Nonetheless, it has been argued that complex models with numerous parameters may not necessarily be advantageous without having enough data for a meaningful validation [32]. More recent studies on the mathematical modeling of the COVID-19 pandemic can be found in the Refs. [33– 39]. After going through the sixth wave in many countries, the world data [40] show that a comprehensive description of the entire time series appears to be an impossible task, since each country presents its own characteristics (e.g., diverse lock-down and socialeconomic measures). Therefore, in this work we proceed by studying the first wave for those countries that present data with a similar structure and that can unambiguously be described mathematically using an epidemiological model of the SIR type. The reason for using the SIR model specifically in this work is because it is the simplest possible model to describe the evolution of an infectious disease. In our case we intend to describe mortality regardless of the detailed description of the intermediate compartments through which individuals may pass. One of the objectives is to determine which countries, if any, can be described with the SIR model, to obtain its parameters and, finally, to compare between the countries to see if some kind of universality of the model can be observed. By inspection, from the recorded worldometer data [40] we found that there are only nine or ten such countries (we leave aside the case of China that has been exhaustively studied and where the pandemic apparently died out without the need for vaccines). Among those countries there are eight Europeans— Spain, France, Italy, UK, Germany, Belgium, Switzerland and Sweden—together with Canada and USA. Moreover, we have also added the cases of India and Brazil. In this work we will carry out a systematic study of the pandemic in each of them by studying the series of cumulative deaths and daily deaths. Our hypothesis is that deaths, D(t), can be considered a fraction of removals, R(t), both cumulative and daily, and therefore, both functions follow an epidemiological curve that will differ essentially in a normalization constant and a shift in time. Our purpose is to investigate if data can be described with simple epidemiological curves using the SIR model and the even simpler Gumbel function [41]. A fundamental question about the first wave of COVID19 is whether the lockdown limitations had an effect in reducing the number of deaths. Non-pharmaceutical interventions (NPI) are still under debate. A recent meta-analysis review [42] fails to confirm that lockdowns have had a large, significant effect on mortality rates. If the daily mortality curves fit well with a basic SIR model, it would be interesting to conclude affirmatively or negatively regarding the effect of NPI on them. The structure of the paper is as follows. In Sect. 2 we review the solutions of the SIR model that will be considered here. We describe in detail how to obtain numerical solutions as a function of the proper time, depending on the parameters βand the basic reproduction number ρ=R0=λ/β. In Sect. 3We examine how well the Gumbel function fits exact SIR solutions with only one parameter, barring normalization and a temporal shift. In Sect. 4we will discuss the time scaling of SIR solutions and define an asymmetry parameter that depends linearly on ρand therefore can be used to characterize the value of the basic reproduction number from a set of data. In Sect. 5we present our results of fits of death data with the exact SIR model and Gumbel functions. In Sect. 6we draw our conclusions. 2 Solutions of the SIR model and the proper time In this section we briefly describe the SIR model and discuss its analytical solution in terms of, what we will call here, proper time,τ, which is a natural variable to measure time through the proportion of removals, where the SIR equations have a trivial and easily interpreted solution. The real time is then obtained by integrating the exact solution. In the SIR model the individuals of a closed population Naffected by a contagious disease are divided into three types: susceptible, S, infected, I, and removed (recovered or dead), R. As functions of time, the num123 Systematic description of COVID-19 pandemic ber of individuals in each compartment is assumed to verify the following equations dS dt=−λ NIS,(1) dR dt=βI,(2) I=N−R−S.(3) In this transmission-dynamics system the first equation means that the variation of susceptible individuals decreases by infection and is proportional to the number of susceptible and the number of infected individuals. The constant λmeasures the rate of infection. The second equation describes the removal variation as proportional to the number of infected individuals, the removal rate being β. By the third equation, the difference between the total number minus the susceptible individuals minus the removed ones must be the number of infected at each instant t. We will consider the initial values S(0)=S0<N and R(0)=0. Therefore, I(0)=N−S0>0. There must be a number, albeit small, of infected in the system initially for the epidemic to begin. For convenience below we will work with the percentages of susceptible, infected and recovered individuals, over the total number of the population, which are obtained by dividing by N: s(t)=S N,i(t)=I N,r(t)=R N,(4) with s(0)=s0,i(0)=1−s0, and r(0)=0. The proper time, τ, is defined as the temporal variable that describes naturally the evolution of the epidemic, by counting the evolution of the recovered individuals, r(t), which is always an increasing function with time. From the second SIR equation in differential form, it is defined by dR=βIdt≡Ndτ, (5) So, the definition of proper time is dτ=drand the interpretation of this variable is that we measure the change in time using the change in recoveries as the biological clock instead of using the physical clock. If we demand that τ=0fort=0, we trivially have τ=r(t). (6) The idea of proper time is based on other equivalent approaches described, e.g., in [43,44], where the susceptible function sis used as variable instead of r. The proper time is nothing more than a change of the time variable into a more convenient one. In our case, time is measured by counting the number of recovered (in percent), since r(t)is an increasing function, although it does not depend linearly on time. Note that the recovered function verify 0 ≤R≤Nand therefore 0≤r(t)≤1. Thus, by definition the proper time has a range limited by 0≤τ≤1.(7) From Eq. (5)wehave I=N β dτ dt,i=1 β dτ dt.(8) Thus, the change of the physical time is given by dt=1 β dτ i(τ),(9) where i(τ) are the infected percent expressed as a function of proper time. To obtain the susceptible function note that we can write, inserting Eq. (5) into Eq. (1) dS=−λ NISdt=−λ βSdτ=−ρSdτ, (10) where ρis the so-called basic reproduction number ρ≡R0≡λ β.(11) The parameter ρhas here the meaning of being the decay constant of the susceptible population in units of proper time. Equation (10) is readily integrated giving S=S0e−ρτ,(12) Thus, Sfollows an exponential decay law as a function of proper time. The constant τ1/2=ln 2 ρ=βln 2 λ(13) represents the half-life in proper time units, i.e., the length of proper time after which the susceptible population is reduced to half. Finally, the third SIR equation (3) gives directly the infected population as a function of proper time I=N−Nτ−S0e−ρτ (14) i(τ) =1−τ−s0e−ρτ,(15) with the condition 0 ≤i≤1. For τ=0 we obtain the initial number of infected population i0=1−s0. The end of the epidemic is reached when i=0. This happens for a value of the proper time τ=τf>0 that is the solution of the transcendental equation i(τ f)=0⇒ τf+s0e−ρτf=1 (16) 123 J. E. Amaro Fig. 1 Solution of the SIR equations as a function of proper time τfor initial susceptible s0=0.99 and for several values of the basic reproduction number ρ=R0 Fig. 2 The peak values of i,τas a function of the basic reproduction number ρ.Fors0=0.99 the peak values are almost independent of s0 In Fig. 1we show some numerical examples. The analytical solutions of the SIR equations are plotted as a function of proper time for various values of the basic reproduction number ρ=1.5,2, 3 and 5. In all cases we assume that s0=0.99; i.e., that one percent of the population is initially infected. In Fig. 1we see that the number of infected individuals as a function of τfirst grows to a maximum and then decreases to zero. The maximum of i(τ) is the peak of the epidemic. It is reached for di/dτ=0. Then −1+ρs0e−ρτ =0,(17) and therefore s0e−ρτ =1 ρ.(18) From here the peak of the epidemic verifies τpeak =ln s0ρ ρ(19) speak =s(τpeak)=1 ρ(20) ipeak =i(τpeak)=1−1+ln s0ρ ρ.(21) A condition for the existence of this maximum for τpeak >0, from Eq. (19), is that s0ρ>1. If we assume that, at the very beginning of the epidemic, s0is very close to 1, then it is enough that ρ>1 for the epidemic to begin to grow [45]. In that case, i(τ) starts growing up to a maximum reached at τpeak,whereitstartsto decrease to zero. In Fig. 2we show the peak values of the infected rate ipeak(ρ) in the epidemic as a function of the basic reproduction number ρ,fors0=0.99. Since s0is very close to one and the dependence on s0is logarithmic, these peak values are almost independent on the precise value of s0≃1. The peak value of infected individuals grows with the basic reproduction number ρ. When ρ is very large, above 10, the logarithmic dependence on the numerator makes the peak to grow more slowly. For ρ=10 we have ipeak =0.67, that is, at the peak of the epidemic two thirds of the population will be infected simultaneously. For ρ=50, more than 90% of the population will be infected simultaneously at the peak. In Fig. 2we also show the value of the proper time at the epidemic peak, τpeak(ρ). It presents a maximum for dτpeak dρ=1−ln s0ρ ρ2=0(22) ln s0ρ=1⇒ ρ=e s0 ≃2.72,(23) for s0≃1. The maximum value of τpeak is then (τpeak)max =s0 e≃0.367 (24) for ρ=e=2.72. Finally, starting from the analytical solution of the SIR equations as a function of the proper time, i(τ), 123 Systematic description of COVID-19 pandemic Fig. 3 Solution of the SIR equations as a function of physical time, t, for initial susceptible s0=0.99, for β=0.3d −1,and for several values of the reproduction number ρ=R0 we will proceed to obtain the solution as a function of the physical time, t. It is obtained from Eq. (9)by integrating between 0 and τ. Assuming that t=0for τ=0, we obtain t(τ) =1 βτ 0 dτ i(τ) =1 βτ 0 dτ 1−τ−s0e−ρτ (25) This integral is not analytical, but it can be calculated numerically with precision using any numerical algorithm, such as Simpson’s rule or Gaussian integration, since the function inside the integral is quite smooth. The final SIR solution is then obtained in parametric form by tabulating t(τ ),i(τ ),τand s(τ ). By plotting i(τ),τand s(τ) as a function of t(τ) we obtain the results of Fig. 3, corresponding to the exact (numerical) solution of the SIR equations as a function of physical time, for the same parameters of Fig. 1and for β=0.3 d−1. Note in Fig. 3, firstly, that the height of the maximum of i(t)coincides with that of the maximum of the analytical solution i(τ) of Fig. 1, as it should be, since we have only made a change of variable—from the proper time to the physical time. Of course, the dependence on physical time has drastically changed. For a constant recovery rate β=0.3d −1, by increasing the value of the basic reproduction number, ρ, the epidemic passes more quickly and lasts less time, almost explosively for ρ=5—for which it lasts only 20 days—compared to ρ=1.5—where it lasts almost two months. Note that the time in Eq. (25) is inversely proportional to β. Making βsmaller, time becomes larger and with it the duration of the epidemic. This will be seen in more detail in Sect. 3.1 when we discuss the time scaling property. Therefore, assuming that the recovery rate, β,inan epidemic is somewhat constant, the basic reproduction number ρlargely determines the evolution of the pandemic. The analytical solution allows estimating the maximum number of simultaneous infections at the peak, ipeak, from this number. This does not depend much on s0, as long as this number is close to one, due to its dependence on the logarithm of soρ. Finally, note that the exact solution i(τ) of the SIR model is asymmetric. It rises very fast at first— exponentially actually—and then falls more slowly. The most explosive, severe outbreaks occur for ρ≥5, where the end of the epidemic is reached for τf≃1 (see Eq. (16)); thus, most of the initial susceptible individuals have been infected. For lower reproduction number the epidemic curve becomes more symmetrical and it is less explosive. Note that in the first wave of the COVID-19 pandemic daily deaths rose very rapidly and fell more slowly, indicating a high reproduction number, so that the data could not be fitted with logistic-type functions, which are symmetric, and a linear combination of two logistic functions was required to fit the data [18,19,28]. In the next section we will see that some more appropriate functions to describe this situation are the Gumbel functions, because they have the adequate asymmetry to fit almost correctly the solutions of SIR equations. 3 Approximation to the SIR solutions with Gumbel functions In Ref. [41] a Gumbel distribution was used to forecast the series of daily deaths with COVID-19 positives. The Gumbel distribution describes the probability of maximum (or minimum) values from the data of many observations [46–48]. It is not theoretically clear why the distribution of deaths or infections roughly recreates the distribution of maxima. In this section we will compare the Gumbel distribution to the exact solutions of the epidemiological SIR model. Both functions of 123 J. E. Amaro Fig. 4 Gumbel distribution compared to the exact Solution of the SIR equations for s0=0.99, β=0.3d −1, and for several values of the reproduction number ρ time presents similar asymmetry and their parameters can be related numerically. Although a Gumbel distribution does not exactly verify the SIR equation, the exact solution can be approximated quite well by the Gumbel distribution by choosing the parameters appropriately. In this section we will provide formulas that relate the parameters of SIR model with the parameters of an appropriate Gumbel function. For this end we will use the proper time method introduced in the previous section. The Gumbel function is defined here as a timedependent function with three parameters G(t)=ae−e−(t−t0)/b(26) and the derivative gives the Gumbel distribution g(t)=dG dt=1 be−(t−t0)/bG(t). (27) First, we will see how well the Gumbel distribution approximates the infected function of the SIR model, g(t)≃i(t). In Fig. 4we compare the Gumbel distribution g(t) with the exact SIR solution i(t). The parameters of the Gumbel distribution have been fitted with leastsquares method to obtain the optimal distribution that best describes the exact solution. The parameters of the SIR model are s0=0.99, β=0.3, and we use four values of ρ=1.5,2,3 and 5—the same as in Fig. 3. The fitted parameters of the Gumbel distribution are given in Table 1. The fits of Fig. 4provide fairly good approximations to the exact SIR solutions with the Gumbel distribution. The results in Fig. 4provide a validation of the Gumbel distribution to describe the evolution of an epidemic. Later, in Sect. 5of results, we will see that this assessment is corroborated in the case of the description of the COVID-19 data. Although the fit is not perfect, the Gumbel distribution clearly shows similar asymmetry as the SIR model. The similarity is greater when the epidemic is explosive, with a high reproduction num123 Systematic description of COVID-19 pandemic Table 1 Parameters of the Gumbel distribution obtained by least squares fit to the exact solution of the SIR model, for four different values of the basic reproduction number ρ(see Fig. 4) s0ρβat 0b 0.99 1.5 0.3 2.04 20.510.9 22.6714.96.29 33.068.81 3.75 53.125.30 2.44 ber ρ∼3–5. Where Gumbel fails most in these fits is in the initial and final stages of the Pandemic. To investigate the connection between the Gumbel distribution and the solution of the SIR equations, it is convenient to express them as functions of proper time. First, we substitute i(t)by g(t)in the definition of proper time dτ=βidt=βgdt=βdG dtdt(28) Thus integrating between 0 and twe have τ=t 0 βdG dtdt=β[G(t)−G0](29) where G0=G(0). Hence G(t)=G0+τ β=ae−e−(t−t0)/b(30) Taking the logarithm on both sides, we obtain e−(t−t0)/b=ln a G0+τ/β (31) from where we obtain the following approximate relation between tand τ t=t0−bln ln a G0+τ/β,(32) and we can write the Gumbel distribution as a function of τ g=1 be−(t−t0)/bG=1 bG0+τ βln a G0+τ/β (33) Starting with this expression as a function of proper time, we can consider several alternative ways of adjusting the Gumbel parameters from the parameters of the SIR model, assuming that the proper time is the same in both models. 3.1 Proper time fit 1 The idea of the proper time fit consists in imposing that the maximum of g(τ) coincides with the maximum of i(τ).Infit1,wewillalsoassumethatG0is very small can be neglected, G0≃0, i.e., we will not try to impose any additional condition on the initial value. This allow us to estimate the parameters aand beasily, but we will not be able to obtain the value of t0, which will later be adjusted to fit the temporal peak of the SIR solution. We start by writing Eq. (33)forG0=0 g(τ) =τ bβln aβ τ.(34) To find the maximum of this function, we compute the derivative g(τ) =1 bβln aβ τ−1 bβ(35) from where we find the maximum condition ln aβ τ=1.(36) Thus, the value of the peak position, τpeak,is τpeak =aβ e(37) and the height of the peak (maximum of g)is g(τpeak)=1 e a b(38) Comparing to the peak values of the SIR solution, Eqs. (19,21), and equating the values we obtain aβ e=ln s0ρ ρ(39) 1 e a b=1−1+ln s0ρ ρ.(40) From here we obtain the values of the parameters aand bof the Gumbel distribution in terms of the parameters of the SIR model a=elns0ρ λ(41) a b=e1−1+ln s0ρ ρ(42) We see that the values of aand bcan be calculated directly from the parameters and initial conditions of 123 J. E. Amaro the SIR model. The value of t0would be fitted later to the position of the peak as a function of time. Note that if the initial number of infected individual is very small, then s0can be approximated by one in the logarithm ln s0ρ≃ln ρ. Then the parameters aand bdo not depend appreciably on the precise value of s0, but only on ρand λby the simple relations a=elnρ λ(43) a b=e1−1+ln ρ ρ(44) Note also that Eq. (44) can be solved numerically to obtain the reproduction number ρas a function of a/b, and then Eq. (43)givesλ. So, the constants λand ρof the SIR system can also be computed from the values of the constants aand bof the Gumbel distribution. 3.2 Proper time fit 2 A second fit can be done by fixing also the number of initial infected i(0), in which case a value for t0can be theoretically obtained. We proceed as in previous section, by computing the maximum of the Gumbel distribution, as a function of proper time, but in this case we use the exact expression, Eq. (33). The maximum of g(τ) is now obtained for ln a G0+τ/β =1 (45) from where the peak (maximum of g) is reached at the proper time τpeak =aβ eG0β(46) and the value of the maximum is g(τpeak)=1 e a b(47) this is the same value for the peak obtained in Eq. (38). Comparing to the peak values of the SIR solution, Eqs. (19,21), and equating the values we obtain aβ e−G0β=ln s0ρ ρ(48) 1 e a b=1−1+ln s0ρ ρ(49) Equation (49) gives the same value for a/bobtained in fit 1. For s0≃1, this value only depends on the basic reproduction number ρ. A third condition is obtained if we demand g(0)= i(0)=1−s0.UsingEq.(33) for the Gumbel distribution as a function of proper time, we have g(0)=1 bln a G0 G0=1−s0(50) where G0=ae−et0/b=ae−x(51) where we have defined the parameter x=et0/b.(52) In terms of the x-parameter, the initial value condition can be written as g(0)=a bxe−x.(53) The procedure of fit 2 follows the following steps: 1. Compute a/bfrom Eq. (49) a b=e1−1+ln s0ρ ρ.(54) 2. Write the initial condition g(0)=i(0)as xe−x=b a(1−s0)≡0(55) and solve this equation numerically for x>0 3. Once we know the values of xand a/b. compute the value of ausing Eq. (48) a=ln s0ρ λ(e−1−e−x)(56) 4. Compute b=(b/a)a, and finally 5. Compute t0as t0=bln x.(57) In Figs. 5and 6the solution of the SIR model is compared with the Gumbel distributions corresponding to the three fits considered in this work: the least square fit, and the proper time fits 1 and 2. The least square fit, with the parameters of Table 1, is represented with dotted lines as a function of proper time in Fig. 5and is very close to i(τ ). We see in Fig. 5 that the proper time fit 1 and fit 2 are essentially the same and that their maximum coincides (by construction) with the maximum of the SIR solution i(τ).The parameter t0in the case of fit 1 is taken from Table 1, because it cannot be obtained theoretically. The parameter t0in fit 2 is computed from Eq. (57) in terms of the computed values of band x. We see that fit 1 and 2 in general do not give the end time τfof the epidemic correctly. This is because 123 Systematic description of COVID-19 pandemic Fig. 5 Three Gumbel distributions as a function of the proper time, compared to the exact Solution of the SIR equations for s0=0.99, β=0.3d −1, and for several values of the reproduction number ρ. The parameters of the Gumbel distributions correspond to the mean-square fit, proper-time fit 1 and proper-time fit 2. In proper-time fit 1, the value of t0is taken from Table 1 Fig. 6 Three Gumbel distributions compared to the exact Solution of the SIR equations for s0=0.99, β=0.3d −1,and for several values of the basic reproduction number ρ. The parameters of the Gumbel distributions correspond to the mean-square fit, proper-time fit 1 and proper-time fit 2. In proper-time fit 1, the value of t0is taken from Table 1 123 J. E. Amaro Fig. 13 Total deaths D(t) during the first wave of the COVID-19 pandemic for several countries, compared to the Gumbel function G(t). The parameter bof the Gumbel function fitted to thedataisgiveninTable2. Data are from Ref. [49] related to the duration of the epidemic. The USA data have been fitted until day 150, when the second wave starts to appear. In the case of Brazil and India the fit is performed until day 250. Note that Spain, France and Belgium and Switzerland have similar b-values in the range b∼12−13 days. Italy, UK and Canada are well fitted with b∼18 days. Germany is in between with b∼15, and Sweden is the European country with the highest value of b∼27 days. Cumulative deaths, D(t), are fairly smooth distributions and very similar among different countries because we are dealing with sums —or integrals on the continuous limit. More detailed information is obtained by describing the daily death data D(t),shownin Figs. 14 and 15. Although these data show more fluctuations, they can be fitted well with the Gumbel g(t) distribution, although the fit parameters differ somewhat from those obtained by fitting the Gumbel function G(t), since different functions and data are being involved. the fitted parameter bis in the third column of Table 2. Again we only tabulate the parameter b, because aand t0give simply the relative height and the position of the peaks. in the case of the USA, Brazil and India, we fit two Gumbel distributions, since it is apparent that there are at least two overlapping waves. 123 Systematic description of COVID-19 pandemic Fig. 14 Daily deaths D(t)during the first wave of the COVID-19 pandemic for several countries compared to the Gumbel function G(t).The parameter bof the Gumbel function fitted to the data is given in Table 2.Dataare from Ref. [49] In these cases, in Table 2we tabulate the two values b1 and b2of the two fitted waves. We see that all countries are well fitted with one or two Gumbel distributions, so this function is an optimal candidate to quantitatively describe an epidemic of these characteristics with only one parameter, b,plus the normalization and position of the peak. The fact that the nine countries considered with an isolated first wave (Spain, France, Italy, UK, Germany, Canada, Belgium, Switzerland and Sweden) only require a timeindependent parameter is remarkable. This does not happen in the following waves or in other countries, where the data show different behavior with large overlaps and stochastic fluctuations. Since Gumbel provides a good analytical approximation to the SIR model solution, it is natural to wonder if the exact SIR solution would give an even better description of the data. So in Figs. 16 and 17 we compare data with exact SIR solutions given by the equations of the previous sections. We know from the last section that there is a linear relationship between the asymmetry of data and the basic reproduction number, ρ. This has allowed us to obtain an approximate value of ρand then we have fitted the value of βand a normalization factor to the width and height of the data, respectively, and we have added a shift in time to get the position of the peak. The parameters ρand βare given in Table 2.InFig.17 we have fitted the US data only to where the first peak is clearly seen. For this reason the other two countries, Brazil and India have not been fitted. In Figs. 16 and 17, we also plot the results of the Gumbel distribution, but this time the parameters have not been fitted to the data, but to the respective SIR 123 J. E. Amaro Fig. 15 Daily deaths D(t)during the first wave of the COVID-19 pandemic for several countries compared to the Gumbel function G(t).The parameter bof the Gumbel function fitted to the data is given in Table 2.Dataare from Ref. [49] solution. Note that, from Eqs. (43) and (44) that βb=ln ρ ρ−1−ln ρ(70) From inspections of the numbers given Table 2, columns 4, 5 and 7, we see that this equation is approximately verified. For instance, for ρ=7 (Spain and UK), the right-hand side of Eq. (70) gives 0.48, and 1/β ≃2bis roughly verified from Table 2.Forρ=3 the equation gives 1/β ≃b/1.2 and this is also verified from Table 2for the case of Belgium. We see that in general the fit of the Gumbel distribution to the SIR solutions in these countries is good, although it begins to fail, for high values of ρ, in the tail part, as we have already seen in the previous section. Note in Fig. 17 that the USA data are described with the SIR model up to day 150 and cannot be described further because the second wave begins almost immediately. The Gumbel distribution, as we have already mentioned, fails to describe the tail of the SIR solution as a result of the least squares fit, which tends to fit best in the region around the maximum. From Figs. 16 and 17 we conclude that the data of the countries considered are globally well described with an exact SIR solution, without the need for any time dependence of the parameters. Again we underline that this only happens in the first wave of the countries that we have considered here and not in the other countries or in the remaining waves. The cause of this requires a detailed study of the epidemiological causes that is beyond the scope of this paper. Finally, in Fig. 18 we plot the values of the parameters of the SIR solution in the (ρ, 1/β) plane. If we fit 123 Systematic description of COVID-19 pandemic Fig. 16 Daily deaths D(t)during the first wave of the COVID-19 pandemic for several countries compared to the Gumbel function G(t).The parameter bof the Gumbel function fitted to the data is given in Table 2.Dataare from Ref. [49] a right line to these data, we see that the general trend seems to be for 1/β to increase with ρ. If we exclude the ’outsider’ countries, Belgium and Switzerland (which have small values of ρ=3 or 4) and Sweden and USA (with large values of ρ=7, 8), we are left with six countries with similar parameters, in the central zone, around ρ=6 and 1/β =20 d —Spain, France, Italy, UK, Canada and Germany. In these countries a relationship between βand ρis not found. On the other hand, the fact that βis similar in these six countries indicates that the probability of recovery is similar in all of them. The recovery probability βfor an individual should in principle be independent of location. But if we consider that there could be important effects due to medical treatments and hospital capacity to treat severe cases in different countries, this could explain the differences between the βvalues. Note that in most countries the value of ρexceeds five, indicating an explosive increase in the first stage of the epidemic in each country. An important consequence of the SIR equations for these high values of ρ is that the total number of people infected during the epidemic reaches almost 100%. Indeed, this is mathematically given by the value of the proper time τat the end of the epidemic, which is the solution of Eq. (16). Numerically it is easy to verify that this end value is practically one, for ρ>5. This can also be seen in the lower panel of Fig. 7, where the value of τ(t)=r(t)for large tis practically one. As a consequence, our results indicate that the data from these countries where ρ>5 are compatible with an epidemic where practically all initially susceptible individuals were infected, according to the SIR model. 123 J. E. Amaro Fig. 17 Daily deaths D(t)in the first wave of coronavirus pandemic for several countries, compared to the Gumbel function G(t). The parameter bof the Gumbel function fitted to thedataisgiveninTable2. Data are from Ref. [49] Fig. 18 Locations of the different countries studied in the plane of SIR parameters (ρ, 1/β) At this point we can already link with the question raised in the introduction of this work, which is whether the lock-downs and other restrictions over the population had any effect in reducing mortality. According to the meta-data study of Ref. [42], the NPI had practically no effect on mortality. This seems to be corroborated in our study for three reasons: (i) that the data are compatible with SIR solutions with high value (R0=ρ≃6) of the basic reproduction number, which implies that all the susceptible individuals were infected; (ii) that the SIR parameters do not depend on time, but if there were some NPI effects the parameters should be time dependent and the epidemiological curves should differ from the SIR solutions; (iii) there does not seem to be a relationship between the intensity of the lockdown measures and the basic reproduction number. For example, Spain, which had very harsh restrictions, is fitted with the same reproduction number as Sweden, which practically did not have, and is greater than Italy, where the measures were introduced a week earlier. 6 Conclusions In this work we have systematically studied the data from the COVID-19 pandemic using the simplest epidemiological model, the SIR model. With the data from the first wave already consolidated and with a perspective of more than two years, we are in a position to analyze the time series of daily mortality data, verify the validity of the SIR model in this pandemic and extract the parameters of the different countries. The SIR model describes the evolution of an epidemic based 123 Systematic description of COVID-19 pandemic solely on statistical laws of proportionality in a sample, and with sufficient data the curve of infected follows a precise and characteristic mathematical law—in fact, in the countries studied we have seen that it describes the data quite accurately. In the first part of this work we have reviewed the mathematical formalism of the SIR model. First we have solved the differential equations of the SIR model in a parametric way using the proper time as a parameter, defined as the relative number of recovered individuals τ=r(t). As a function of τ, the SIR solution is analytical, which allows us to study some of its properties, such as, for example, calculating the maximum number of infected ipeak and the asymptotic number of recovered at the end of the epidemic r(∞). Secondly, we have studied the possibility of approximating the SIR solutions using Gumbel distributions g(t), because this family of functions presents a similar asymmetry as the SIR solutions, and only depends on one parameter, plus the normalization and the position. We have proposed various methods of fitting Gumbel distributions to exact SIR solutions. In particular, using the proper time, we have found simple relationships between the Gumbel parameters and the SIR model parameters. Third, we have discussed the scaling properties of the exact SIR solutions when plotted as a function of βt, where βis the probability of removal per unit of time. Next we have defined an asymmetry parameter, as the ratio between the right and left half-widths at halfheight, of the SIR solutions. We have shown numerically that the asymmetry A(ρ) grows almost linearly with the reproduction number ρand that it is independent of β. The asymmetry uniquely characterizes the value of ρand vice versa. Therefore, a measure of the asymmetry of some data of i(t)at the middle of the height allows to extract the value of ρ. In the results section we have applied the SIR model and the Gumbel distribution to study the daily deathdata in the first wave of the COVID-19 pandemic in a dozen countries. The countries have been selected because they are the only ones that present a peak that closely resembles a SIR solution. The data from Spain, France, Italy, UK, Canada, Germany, Belgium, Switzerland and Sweden can be fitted quite well with a SIR solution and also with a Gumbel function. Except for Belgium and Switzerland, data in the rest of the countries are compatible a reproduction number ρ>5. This seem to indicate that in practically all susceptible individuals were infected and eventually recovered. This raises questions about the effectiveness of nonpharmaceutical interventions, such as lock-downs, in many countries. In short, the success of the SIR model to describe the first wave of the COVID-19 pandemic in the countries analyzed has not only allowed us to extract the two parameters that govern the temporal evolution, the basic reproduction number and the constant removal rate, β. It has also made it possible to carry out a comparative study between the different affected countries. Acknowledgements The author thank Dr. Nico Orce for critical reading of the manuscript. This work has been supported by the Spanish Agencia Estatal de investigacion (D.O.I. 10.13039/501100011033, Grant No. PID2020-114767GB-I00) and the Junta de Andalucia (Grant No. FQM-225). Funding Funding for open access charge: Universidad de Granada / CBUA. Data availability The datasets analyzed during the current study are available in the “ourworldindata” repository, https:// ourworldindata.org/coronavirus. Declarations Conflict of interest The author declare that he has no known competing financial interest or financial relationships that could have appeared to influence the work reported in this paper. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/. References 1. Hui, D.S., Azhar, E., Madani, T.A., et al.: The continuing 2019-nCoV epidemic threat of novel coronaviruses to global health—the latest 2019 novel coronavirus outbreak in Wuhan, China. Int. J. Infect. Dis. 91, 264 (2020) 2. Tang, B., Wang, X., Li, Q., Bragazzi, N.L., Tang, S., Xiao, Y., Wu, J.: Estimation of the transmission risk of the 2019nCoV and its implication for public health interventions. J. Clin. Med. 9(2), 462 (2020) 123 J. E. Amaro 3. Wu, J.T., Leung, K., Leung, G.M.: Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study. Lancet 395(10225), 689 (2020) 4. Kraemer, M.U.G., Yang, C.H., Gutierrez, B., Wu, C.H., Klein, B., Pigott, D.M., Brownstein, J.S.: The effect of human mobility and control measures on the COVID-19 epidemic in China. Science 368(6490), 493 (2020) 5. Hamer, W.H.: The Milroy Lectures on Epidemic disease in England—the evidence of variability and of persistency of type. Lancet 167(4305), 569 (1906) 6. Ross, R.: Report on the Prevention of Malaria in Mauritius. Waterlow and Sons, London (1908) 7. Ross, R.: An application of the theory of probabilities to the study of a priori pathometry. Part I. Proc. R. Soc. Lond. A 92, 204 (1916) 8. Ross, R., Hudson, H.P.: An application of the theory of probabilities to the study of a priori pathometry. Part III. Proc. R. Soc. A 93, 225 (1917) 9. Kermack, W.O., McKendrick, A.G.: A contribution to the mathematical theory of epidemics. Proc. R. Soc. A 115, 700 (1927) 10. Kendall, D.G.: Discussion of ‘Measles periodicity and community size’ by M. S. Bartlett. J. R. Stat. Soc. A 120,64 (1957) 11. Bartlett, M.S.: Deterministic and stochastic models for recurrent epidemics. In: Berkeley Symp. on Math. Statist. and Prob., Proc. Third Berkeley Symp. on Math. Statist. and Prob., vol. 4, pp. 81–109. Univ. of Calif. Press (1956) 12. Bartlett, M.S.: Measles periodicity and community size. J. R. Stat. Soc. A 120(1), 48 (1957) 13. Flanders, W.D., Kleinbaum, D.G.: Basic models for disease occurrence in epidemiology. Int. J. Epidemiol. 24(1), 1 (1995) 14. Weiss, H.: The SIR Model and the Foundations of Public Health, MATerials MATematics no. 3 (2013) 15. Chauhanl, S., Misra, O.P., Dhar, J.: Stability analysis of SIR model with vaccination. J. Comput. Appl. Math. 4(1), 17 (2014) 16. Chao, D.L., Dimitrov, D.T.: Seasonality and the effectiveness of mass vaccination. Math. Biosci. Eng. 13(2), 249 (2016) 17. Rodrigues, H.S.: Application of SIR epidemiological model: new trends (2016). arXiv:1611.02565 18. Amaro, J.E.: The D model for deaths by COVID-19 (2020). arXiv:2003.13747v1 19. Amaro, J.E., Dudouet, J., Orce, J.N.: Global analysis of the COVID-19 pandemic using simple epidemiological models. Appl. Math. Model. 90, 995 (2021) 20. Cooper, I., Mondal, A., Antonopoulos, C.G.: A SIR model assumption for the spread of COVID-19 in different communities. Chaos Solitons Fractals 139, 110057 (2020) 21. Kudryashov, N.A., Chmykhov, M.A., Vigdorowitsch, M.: Analytical features of the SIR model and their applications to COVID-19. Appl. Math. Model. 90, 466 (2021) 22. Postnikov, E.B.: Estimation of COVID-19 dynamics “on a back-of-envelope”: does the simplest SIR model provide quantitative parameters and predictions? Chaos Solitons Fractals 135, 109841 (2020) 23. Fanelli, D., Piazza, F.: Analysis and forecast of COVID-19 spreading in China, Italy and France. Chaos Solitons Fractals 134, 109761 (2020) 24. Radulescu, A., Williams, C., Cavanagh, K.: Management strategies in a SEIR-type model of COVID-19 community spread. Sci. Rep. 10, 21256 (2020) 25. Chen, X., Li, J., Xiao, C., Yang, P.: Numerical solution and parameter estimation for uncertain SIR model with application to COVID-19. Fuzzy Optim. Decis. Mak. 20, 189 (2021) 26. Guerrero Sanchez, Y., Sabir, Z., Guirado, J.L.G.: Design of a nonlinear SITR fractal model based on the dynamics of a novel coronavirus (COVID-19). Fractals 28(08), 2040026 (2020) 27. Elsombaty, A., Sabir, Z., Ramaswamy, R., Adel, W.: Dynamical analysis of a novel discrete fractional SITRS model for COVID-19. Fractals 29(8), 2140035 (2021) 28. Amaro, J.E., Orce, J.N.: Monte Carlo simulation of COVID19 pandemic using Planck’s probability distribution. Biosystems 218, 104708 (2022) 29. Xie, G.: A novel Monte Carlo simulation procedure for modelling COVID-19 spread over time. Sci. Rep. 10, 13120 (2020) 30. Allen, J.L.S.: An introduction to stochastic epidemic models. In: Brauer, F., van den Driessche, P., Wu, J. (eds.) Mathematical Epidemiology, pp. 81–128. Springer (2008) 31. Andersson, H., Britton, T.: Stochastic Epidemic Models and Their Statistical Analysis. Springer, Cham (2000) 32. Roda, W.C., Varughese, M.B., Han, D., Li, M.Y.: Why is it difficult to accurately predict the COVID-19 epidemic? Infect. Dis. Model. 5, 271 (2020) 33. Nauman, A., Elsonbaty, A., Raza, A., Rafiq, M., Adel, W.: Numerical simulation and stability analysis of a novel reaction–diffusion COVID-19 model. Nonlinear Dyn. 106(2), 1293 (2021) 34. Khan, M.A., Atangana, A.: Mathematical modeling and analysis of COVID-19: a study of new variant Omicron. Phys. A Stat. Mech. Appl. 599, 127452 (2022) 35. Masum, M., Masud, M.A., Adnan, M.I., Shahriar, H., Kim, S.: Comparative study of a mathematical epidemic model, statistical modeling, and deep learning for COVID-19 forecasting and management. Socio-Econ. Plan. Sci. 80, 101249 (2022) 36. Srivastav, A.K., Tiwari, P.K., Srivastava, P.K., Ghosh, M., Kang, Y.: A mathematical model for the impacts of face mask, hospitalization and quarantine on the dynamics of COVID-19 in India: deterministic vs. stochastic. Math. Biosci. Eng. 18, 182 (2021) 37. Tiwari, P.K., Rai, R.K., Khajanchi, S., Gupta, R.K., Misra, A.K.: Dynamics of coronavirus pandemic: effects of community awareness and global information campaigns. Eur. Phys. J. Plus 136, 994 (2021) 38. Rai, R.K., Khajanchi, S., Tiwari, P.K., Venturino, E., Misra, A.K.: Impact of social media advertisements on the transmission dynamics of COVID-19 pandemic in India. J. Appl. Math. Comput. 68, 19 (2022) 39. Majumder, M., Tiwari, P.K., Pal, S.: Impact of saturated treatments on HIV-TB dual epidemic as a consequence of COVID-19: optimal control with awareness and treatment. Nonlinear Dyn. 109, 143 (2022) 123 Systematic description of COVID-19 pandemic 40. https://www.worldometers.info/coronavirus/ 41. Furutani, H., Hiroyasu, T., Okuhara, Y.: Simple method for estimating daily and total COVID-19 deaths using a Gumbel model. Researchsquare (2020). https://doi.org/10.21203/rs. 3.rs-120984/v1 42. Herby, J., Jonung, L., Hanke, S.H.: A literature review and meta-analysis of the effects of lockdowns on COVID-19 mortality. Stud. Appl. Econ. 200, 1 (2022) 43. Harko, T., Lobo, F.S.N., Mak, M.K.: Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates. Appl. Math. Comput. 236, 184 (2014) 44. Miller, J.C.: A note on the derivation of epidemic final sizes. Bull. Math. Biol. 74(9), 2125 (2012) 45. Diekmann, O., Heesterbeek, J.A.P., Metz, J.A.J.: On the definition and the computation of the basic reproduction ratio R0in models for infectious diseases in heterogeneous populations. J. Math. Biol. 28, 365 (1990) 46. Gumbel, E.J.: Les valeurs extrêmes des distributions statistiques. Annales de l’Institut Henri Poincaré 5(2), 115 (1935) 47. Gumbel, J.E.: The return period of flood flows. Ann. Math. Stat. 12, 163 (1941) 48. Gumbel, E.J.: Statistical theory of extreme values and some practical applications. U.S. Department of Commerce, National Bureau of Standards. Applied Mathematics Series, vol. 33, 1st edn (1954) 49. https://ourworldindata.org/coronavirus 50. https://covid19.who.int/ Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123