Monitoring and Forecasting COVID-19: Heuristic Regression, Susceptible-Infected-Removed Model and, Spatial Stochastic

The COVID-19 pandemic has had worldwide devastating effects on human lives, highlighting the need for tools to predict its development. The dynamics of such public-health threats can often be efficiently analyzed through simple models that help to make quantitative timely policy decisions. We benchmark a minimal version of a Susceptible-Infected-Removed model for infectious diseases (SIR) coupled with a simple least-squares Statistical Heuristic Regression (SHR) based on a lognormal distribution. We derive the three free parameters for both models in several cases and test them against the amount of data needed to bring accuracy in predictions. The SHR model is ≈ ±2% accurate about 20 days past the second inflexion point in the daily curve of cases, while the SIR model reaches a similar accuracy a fortnight before. All the analyzed cases assert the utility of SHR and SIR approximants as a valuable tool to forecast the disease’s evolution. Finally, we have studied simulated stochastic individual-based SIR dynamics, which yields a detailed spatial and temporal view of the disease that cannot be given by SIR or SHR methods.

The consequences of a pandemic like COVID-19 caused by the virus SARS-CoV-2 cannot be overstated (Nature, 2021).Accurate mathematical tools allowing to monitor and forecast the evolution of the contagious disease are useful to guide social, economic and public health decisions made by governments.Nevertheless, despite the availability of powerful mathematical models (Anderson et al., 2020), initial forecasting by some organizations underestimated the evolution of the epidemics, hampering the immediate taking of necessary actions (Economist, 2020;Herzberge and Hecketsweller, 2020).
This study aims to take advantage of available worldwide data on COVID-19 (Roser et al., 2021) (Dong et al., 2020) to benchmark and assign error bars to minimal models, like the susceptible-infected-recovered (SIR) with different sophistication levels (Annas et al., 2020;He et al., 2020a;Kermack and McKendrick, 1927;Khan et al., 2020;Weiss, 2013;Yang and Wang, 2020), a straightforward least-squares best-fit (LS) Statistical Heuristic Regression (SHR) based on a lognormal distribution (Lam, 1988), or basic Monte-Carlo simulation (Gang, 2020;Girona, 2020).It is well-known that finding a global minimum of non-linear least-squares problems for p free parameters requires, at worst, a brute force search in p-dimensional parameter space.If each parameter can take m values inside a given interval, it is a non-polynomial task that scales like m p and becomes non-practical for large or moderate values of p. Correspondingly, there are no general polynomial bounds on the time complexity given in the number of samples and the search space dimension.These models are gauged against two variables measured daily: (i) the number of deaths, and (ii) the number of new infections.Such indicators both possess advantages and disadvantages.
Deaths are usually counted using a consistent methodology 1 and, undeniably, it is an observable proportional to the spread of the disease, but the tally of deaths carry a delay of about one month on the actual dynamics of the disease.On the other hand, the number of infections is timely, but incorporates more uncertainties since it depends on details not related to the disease, e.g. on the number of tests performed.We show the simultaneous monitoring of both observables supplemented with relatively simple mathematical approaches can be used to follow and forecast the evolution of the disease with enough accuracy to help decision-making processes and we discuss the associated error bars.
1 However, some countries have varied these criteria during the course of the epidemic.Moreover, different countries can use different rules to define these variables.Those inconsistencies can be compensated and do not significantly affect the conclusions of our analysis.However, a comparison between countries is only an approximated one.
Other efforts to modeling the pandemics include sensitivity and meta-analysis to estimate averaged values for the reproduction number, incubation time, infection rate and fatality rate (He et al., 2020b), wavelet-coupled random vector functional link networks (Hazarika and Gupta, 2020), machine learning (Dhaka and Singh, 2020), and Advanced Autoregressive Integrated Moving Average Model (Singh et al., 2020).Approaches via learning algorithms are usually compared via corresponding tests (Demsar, 2006), where we recall the significant differences to statistics (Breiman, 2001).
The paper is organised as follows.In Section 2 (results and discussion), we introduce the Statistical Heuristic Regression model (SHR, section 2.1), the Susceptible-Infected-Removed model (SIR, section 2.2), and the Spatial Stochastic Individual-Based model (MC, section 2.3).After each topic, we analyse the corresponding application to different countries or regions, most notably Spain and Germany.Finally, in the section of Conclusions, we summarise and review our approaches and further discuss the application of the cases analysed and future venues.

A. Statistical Heuristic Regression (SHR)
Epidemics can efficiently be modelled as a geometric process related to independent random events (Lam, 1988).This method yields a regression curve that describes the temporal variation of a contagious disease for the number of deaths, infections or some other relevant observable variable.Such a statistical heuristic approach results in a lognormal function, which is the probability distribution function of a random variable whose logarithm, u = ln (t), is normally distributed around its mean value µ with a dispersion σ (Johnson et al., 1994).The beginning of the propagation is determined by t 0 and the value of the single maximum c M = c(t M ) happens at t M = t 0 + µ − σ 2 .
Starting from the model dc(t) dt = α(t)c(t) and imposing general requirements on α(t) (which follow from the observed behaviour of the number of daily cases) also leads to the same expression (Wenbin et al., 2013).Using entropy-related arguments, these authors have estimated that σ ≈ 0.4, which compares well with the averaged values for ten different western countries for deaths and infected, 0.6 ± 0.2 and 0.5 ± 0.2 respectively, cf.Tables I   and II.Finally, we notice that such lognormal distribution derived on (Lam, 1988) was proved useful to model the SARS outbreak in 2003 (Chan et al., 2006).
The corresponding accumulated cases are, Given arbitrary precision, C(t) and c(t) carry the same information about the set of three parameters, F = {c M , µ, σ}, since c(t) is simply the temporal derivative of C(t).However, in practical terms, c(t) is heavily affected by noise in the collection of the data series {c i }, and least-squares fits to the functions C(t) and c(t) are expected to determine slightly different values for F .Therefore, we chose to report values related to C(t), which are less affected by noise.Still, we notice that the information contained in c(t) is equally valuable and sometimes simpler to obtain, in particular, the position and value of its inflexion points and single maximum.
Next, we aim to prove that the ansatz in Equations 1 and 2 reproduces the behaviour of COVID-19 in ten different western countries using actual data up to the time of submission (revised version, March 2021).We observe a first wave that is relatively homogeneous among all countries (if properly normalized) and that could be considered strongly mitigated around May-June everywhere (see Section I.A.5).Other waves have later developed, which are more heterogeneous because they reflect each country's different responses to the epidemics.A superposition of individual elementary peaks has been used to model these ulterior waves.Even if SHR merely amounts to a precise fit of the data, we observe that it carries significant advantages over the mere manipulation of the data series, {c i }, as: (i) it can be extrapolated to the near future (extrapolations should be treated with great care, but an informed extrapolation about the behavior in the future is always better than a wild guess) and, (ii) it reduces long lists of numbers to an analytical expression which only depends on three parameters.Such an analytical function can be then easily manipulated to get integrals, derivatives of any order, or to search for extrema/inflexion points, etc.

Spain.
Spain is a country where the disease was particularly virulent in its first wave, spreading with remarkable strength.The SHR model agrees well with the data for both deaths and infections (Fig. 1 and Tables I and II ; a common colour code is applied in Figs. 1,2,3,4,8,9,10 and 11 to facilitate comparisons: red is used for daily cases (empty dots for actual data and dashed/dotted lines for models), black for 7-days moving averages of daily cases, and blue is used for accumulated cases (again, empty dots for actual data and dashed/dotted lines for models).In Fig. 11 we have used red vs magenta for daily cases and blue vs cyan for accumulated cases to allow easy comparison between models.Other details are given in the figure captions.Together, these two variables provide a better idea of the epidemic's course by identifying two critical items: the impact on the population via infections and, the impact on the health system via deaths.Three simple features defining the epidemics that will be rationalized later in the context of the SIR model are: (i) the exponential behavior near the origin, (ii) the position and value of the single maximum in the daily number of cases and, (iii) an asymmetric decay towards the future w.r.t the past.The ratio between total infections and deaths has evolved from about 1% in March to a maximum of 12% in August, but it has significantly decayed for the second wave to about 4% at the end of September (inset in left-hand side in Fig. 1).
Three regions are identified in the plots, both for deaths and infections.The first region (I, wave 1) finishes approximately on the 1st of May, 2020 (d = 152) and it clearly shows the three aforementioned features marking its association with an infectious disease.The second region (II, inset in right-hand side in Fig. 1), goes approximately between the 150th and 200th day, and its hallmark is to sustain a fairly constant level of daily infections, c(t) =< c i >, which reflects in a linear increase of the accumulated number of cases, C(t).
Region II approximately terminates near the end of the general lockdown in Spain, on the 21st of June (d = 173).
Neither the SHR nor the SIR models can account for a sustained period of constant infections, although they can accommodate this regime via the slowly decaying queue of the distribution where derivative of the function is very low.In contrast, such behavior can be naturally described via Monte-Carlo simulation.Likewise, while MC can describe several waves by producing more than one local maxima due to spatial inhomogeneous dynamics, SHR and SIR can only describe such a scenario via a linear combination of individual waves, each one governed with its own parameters.
Finally, in a third region (III, wave 2) the collective transmission displays again a similar behavior to the region I, marking the evolution of an out-of-control disease.The superposition of these multiple regimes, plus other waves if needed, describes the overall function well.We notice that the accuracy in the fit for any wave is not expected to be reasonably stable until at least the corresponding maximum is well developed (Section I.A.6).However such incertitude, the model predicts that in Spain the number of infections due to the second wave should be reaching its maximum in December 2020, at most in January 2021.2In addition, the model predicts that the strength of the second wave is approximately weaker than the first one by a factor two, as measured by the number of accumulated certified deaths from SARS-Cov-2.Although these predictions may be affected by large error bars since the maximum in the second wave is not yet well developped, those values offer sound guidance about the course of the disease.We have used this model to extrapolate the shape of the curve by a fortnight after the last day of the corresponding available data;3 the resemblance to the ulterior course of the disease will be seen in the next weeks.
The accompanying number of registered infections yield a picture of the likely evolution of deaths in the following days, even if the variation in the absolute numbers from the first to the second wave is dominated by the change in the number of tests performed.Given the large dispersion of raw data due to difficulties to collect them it is clear the necessity to perform moving averages and the advantages of working with least-square approximants that can be extrapolated a few days ahead, a statement that is true for the behaviour of other countries.While deaths only show two waves so far, infections identify at least four local maxima that can be correlated with different events, like the end of the summer vacations or the occurrence of several bank holidays in Spain where the population has been moving and mixing in great numbers.

Germany.
Compared with other countries with large populations, Germany has managed the pandemics quite well, as it is observed by comparing the number of cases per inhabitant.Moreover, its evolution has been recorded with consistency both for deaths and infections.Therefore, it is an appropriate benchmark for any model.
Similarly to Spain, the SHR model can be used to accurately represent the disease evolution using only three parameters per wave (Fig. 2).Curiously enough, best-fit values for µ and σ are quite similar to Spain (Tables I and II), indicating that, independently of the absolute strength, there are common underlying features in both cases.Therefore, it is interesting to explore the ability of a single normalized averaged curve to represent such contrasting cases as Spain and Germany, using C(∞) ∝ c M as the single only free parameter.
Such a curve is represented in Figs. 1 and 2 by the green dashed line having µ = 3.53 and σ = 0.56 (Section I.A.5), and it is clear that despite having such a limited freedom for fitting (since it only depends on one parameter), it provides a very reasonable approximation to the data.In contrast to Spain, the ratio between total infections and deaths in Germany evolved from about 1% in April to 4% in September, which is about three times lower than for Spain (inset in Fig. 2).

Other countries.
We also prove the capabilities and versatility of the SHR ansatz to reproduce the observed data by applying the same methodology to a pool of western countries: Great Britain (GBR), Italy (ITA), United States (USA), France (FRA), Switzerland (CHE), Denmark (DNK), Austria (AUT) and Finland (FIN), cf.Figs. 3 and 4. In general, the agreement is quite good, both for deaths and infections.Among other advantages, this procedure allows a quick and simple monitoring of the evolution of the disease in the different countries.In particular, it is a useful tool to identify and forecast the appearance of a second wave.At the moment of writing, only the USA has fully developed the maximum associated with the second wave and, from the combined behaviour of deaths and infections, it could be argued that the country is clearly heading towards a third wave.Since this is the only case so far, it is not possible to characterize well such a second wave by a proper average of different countries, although it seems fair to say that it is represented by a wider distribution of daily deaths (e.g. the second component represented by the dotted red curve corresponds to having µ = 5.5 and σ = 1.2) and a lower value at the peak by about a factor 2.
Prominent places where the infection spreads quickly are densely populated regions, which constitute the core of the propagation of the disease.Therefore, it is interesting to compare the distribution of cases in those regions.We have juxtaposed the performance of New York City (9.1 M-people, NYC) and the Community of Madrid (6.7 M-people, CAM) in the first wave (Fig. 5).To highlight the similarities rather than the differences they are superimposed in such a way that the position (day) of the maximum coincides.Furthermore, CAM has been scaled by the ratio of respective populations, which makes the value at the maximum very similar for both regions.Despite all the differences between these regions, it is clear that a typical pattern emerges, which leads us to investigate the advantages of working with averages.an answer with a fractional error of ≈ ±5%, which is an excellent initial guess taking into account that it only depends on a single parameter, c M .Such parameter c M can be easily obtained from a single point: the maximum value in the daily distribution of cases for each wave, which we derive from a moving average of a few days (seven days makes appropriate averages that account for regular weekly routines and removes most of the noise for all the cases we have analyzed).

Accuracy of SHR.
To be able to confidently use a least-squares statistical regression to a given data set the main question is how many data points, n, are needed to yield a reasonable estimation of the evolution of the epidemics based solely on the extrapolation of the fitted functions.Such question is relevant considering how unreliable extrapolations usually are (Press et al., 2007).Indeed, any simple algorithm to forecast the evolution of an epidemics can only be valuable if reasonable error bars can be assigned to predictions.
A simple target to quantify the error is to study the behavior of the expected total number of cases, C(∞), as a function of n.Fig. 7 shows the variation of the predicted asymptotic value as a function of the available amount of data after the second inflexion point.In most cases, a fractional accuracy of ±15% is achieved a fortnight after the second inflexion point, which is further decreased to ±5% in another fortnight.
B. Susceptible-Infected-Removed (SIR) So far, we have shown that SHR qualifies as a quick and straightforward way to describe the evolution of an infectious disease.If adequately used, i.e. attached with appropriate error bars, it can be extrapolated to make predictions in the near future, since the functional forms associated with Equations 1 and 2 adapt so well to the observed data.
However, a better understanding of the dynamics of the epidemics can be obtained from a set of differential equations which describe its time evolution.The simplest model for the evolution of a contagious disease is to postulate that the rate of new infections is proportional to the number of infected people itself, dI(t) dt = I(t) τ 0 , which results in an unbound exponential growth, I(t) ∝ e t τ 0 , and makes a characteristic mark for the onset of a pandemic.Such a simple model does not take into account how the rate of infections decreases as the number of infections approaches the total population.Therefore, a refined version is to divide a given population of size N into three classes (S, I, R): (i) susceptible entities who can catch the disease, S(t), (ii) infected ones who have the disease and transmit it, I(t), and (iii) removed ones who have been isolated, died, or recovered and become immune, and are therefore not able to propagate the disease, R(t).In this model, individuals pass from the susceptible class S to the infective class I and finally to the removed class R with rates determined by a set of ordinary differential equations(ODEs) (Anderson, 1991;Hethcote, 2000;Kermack and McKendrick, 1927;Weiss, 2013).The ODEs derive from the interactions of the entities in the different classes, which can be represented as where we assume generalized mass-action kinetics (Müller and Regensburger, 2012) (with slightly different scaling with respect to N).First, it is assumed that the number of susceptible individuals decreases at a rate proportional to the density of infected, i(t) = I(t) N times the number of susceptible individuals, S(t), where τ 0 is an adjustable parameter that represents a typical time to transmit the disease, and n is a parameter that influences the ability of the disease to infect susceptible individuals in a non-linear way (e.g. it might represent the effect of the viral load).Its main effect is to alter the temporal scale of the epidemics, which in some circumstances facilitates the fitting of the model to real data.The standard SIR model is recovered with n = 1.
Removed entities originate from infected; therefore, its variation is assumed to be proportional to the number of infected, where τ 1 is an adjustable parameter that represents a typical time to recover from the disease.
This equation merely helps to count the total number of removed from the beginning of the infection up to a given day t, Lastly, the infected vary according to the inflow of susceptible individuals who become infected minus the outflow of infected that have been removed, The derivative dI dt moves from positive to negative depending on the balance between both terms in the equation and it determines a single peak in I(t) (for n = 1, I ′ (t) = 0 for S(t) N = τ 0 τ 1 ).The task at hand for a given population of N elements is to determine the parameters, τ 0 , τ 1 and n, that best reproduce the behavior of the epidemics by solving the coupled system of differential equations 3 and 6, subject to some initial conditions, e.g.S(0) = N − 1, I(0) = 1.Good agreement with data can be used to lend an interpretative value to τ 0 and τ 1 (unlike parameters µ and σ which only have a statistical meaning).The ratio ℜ 0 = τ 1 τ 0 is called the effective reproductive number; values ℜ 0 >> 1 characterize a virulent disease where R(∞) = S(0).
First, we focus on the task of simulating a population where s(0 1.For this particular case, ℜ 0 >> 1 and the entire susceptible population is removed at the end.The proposed algorithm goes as follows, 1. We use the daily number of deaths to identify the position and maximum value in the infections/deaths data: t * M and i * M . 2. τ 0 is the main parameter that determines the position of the peak in i(t).We estimate a value for τ 0 that brings the maximum in i(t) near t * M .
3. We get an approximate value for τ 1 from the expression i *

The value i *
M yields N in the particular case of R(∞) = N = c M .We adjust the value of N to agree with i * M .
5. We minimize the root-mean square deviation, χ N = , between the number of accumulated cases predicted by the model, R(t), and the recorded data, C i , to find optimal values for τ 0 , and τ 1 .Similarly to the SHR analysis we have presented above, we illustrate the performance of the SIR model by first looking at the distribution of deaths and infections in Spain (Fig. 8).The lower left panel shows how numerical solutions to SIR equations match very well the temporal behavior of the epidemics under the condition s(0) = r(∞) = 1 for optimized values of τ 0 and τ 1 (Table III).Dispersion of data in the daily reported cases is usually smaller before the peak is reached (the quasi-exponential region) and fluctuations grow in importance after the maximum is reached-which is a general observation holding for most of the countries we have studied.We assign it to the balance between different currents transferring individuals between the three classes, the phenomenon responsible for the appearance of a single maximum in daily cases for a given wave in the pandemics.
The proposed procedure works for the deaths subset as follows.First, the curve of daily cases is followed up to the appearance of its maximum, which to circumvent the noise is identified from a smoothed curve obtained by a five-day moving average, I * M (t * M = 96) = 17.1 per million people (the single daily maximum value is I M (t M = 94) = 20.2).The SHR model for accumulated deaths using only data up to six days past the maximum yields a prediction of total deaths of N = 383, which is off the final mark by about 40%.
Once the maximum is identified, the quasi-exponential behavior near the origin is used to estimate an initial value for τ 0 (supplementary material, Eq.S6).For Spain, the first case happens at t = 65, and the first inflexion point is at t 1 = 82.Therefore, the first 10 points (about halfway to t 1 ) are used to get an exponential fit to the accumulated number of cases that yields τ 0 ≈ 2.8 ± 0.3.Such a value, combined with an initial guess ℜ 0 = 10 produces a maximum in the curve of daily deaths at t M = 103.Accordingly, τ 0 is decreased until we locate the maximum closer to the right position.For τ 0 = 2.2 we get t M = 95 and I M = 11.7 (per million inhabitant).Therefore, we update the value of N using the ratio 17.1 11.7 and start an efficient local Levenberg-Marquardt minimization of the root-mean-squared deviation between the actual data and the computed values.This is done to simultaneously optimize N, τ 0 and τ 1 (Fig. 8, left-lower panel).Taking into account that only data up to six days past the maximum have been used, it is remarkable that this self-consistent procedure reduces the fractional error between the prediction of the SIR model and the data from 40% to ±3%, being the root-mean-squared deviation (RMSD) between the accumulated data and the predicted function χ N = 0.6%.Such a low RMSD value matches the good visual agreement observed.We believe that the logic behind the steps proposed above amounts to more than a recipe to get a best fit, yielding meaning to the values obtained and their interpretation.
Next, we explore how the SIR model represents the evolution of the number of infections.
The number of infections is a magnitude that carries larger error bars, but it can provide timely information on the evolution of the epidemics (Fig. 8, right-lower panel shows the case for Spain).As expected, infections start earlier than the deaths(t = 32 vs t = 65), but need more time to attain its maximum value (t = 54 vs t = 25 after the first case).
A prominent feature is the existence of the second wave of infections separated from the first one by a region of sustained constant number of cases, as we have discussed in Fig. 1.To fit the data, we superpose the two waves, each with its own defining parameters.However, the constant region between waves cannot be easily accommodated in these models and it is a clear indication of a different stage in the epidemics with low but sustained transmission of the disease at a pace similar to the one at which individuals are removed (while in the SIR model usually it is assumed that τ 1 > τ 0 ).We shall come back to this point in the context of Monte-Carlo simulation.Finally, we notice that this second wave of infections has finally overlapped with a third one, as it is noticeable in Fig. 1.

Germany.
We have applied the same procedure to Germany, a country which had in the first wave about four times less casualties per million inhabitant than Spain.The left panel of Fig. 9 shows the final iteration for the daily and accumulated number of deaths, which again predicts the total number due to the first wave with accuracy ≈ ±3% of the final true value, even if we have only used data up to six days past the maximum.The RMSD between the accumulated data and the predicted function is χ N = 0.5%, which reflects the good visual agreement observed too.
Regarding the infections, it is interesting that the region of sustained infection is also observed, although a second wave is only weakly apparent up to the present day (t=200).
Again, infections start earlier (t = 28) w.r.t deaths (t = 70).Furthermore, maximum values are attained after a longer amount of time (t = 63 for infections (t = 63) than for deaths (t = 30), counted after the first case, following a similar procedure to the one for Spain.

Other countries.
Finally, similar results have been observed in four more countries: France, Italy, Great Britain and Switzerland (Fig. 10 and Table III).

C. Spatial individual-based model.
To gain further insight into the spatio-temporal evolution of COVID-19, we consider next a stochastic discrete-time individual-based model in which the spread propagates on a two-dimensional N × N lattice, where each node represents an individual.The dynamics are Markovian, and we will use Monte Carlo (MC) to sample from its distributions in time, which is a technique known to handle well difficult collective effects in many-body systems, like e.g. the magnetic phase transition in the 2D Ising model (Peliti, 2011).The N 2 individuals can be in any of the three states of the SIR model, making transitions between them with two probabilities: (i) for someone susceptible to be infected S → I, p i and, (ii) for someone infected to recover and be removed I → R , p r .At each time-step, individuals make transitions between classes according to the corresponding probabilities.We consider various scenarios of uniform and spatially dependent Markov dynamics.
First, we start with a single isolated case of infection per 10 4 individuals, and we use τ 0 for S → I in close analogy to the SIR model, while we assign a constant value p r = 1 τ 1 to the second transition probability, I → R. Comparing MC simulations for N = 100 with p i (t) = 1 2 i(t) and p r = 1 10 to the deterministic SIR with τ 0 = 2.1 and τ 1 = 9 yields an excellent agreement between both approaches for same initial values, which confirms the adequacy of Monte-Carlo techniques (left panel in Fig. 11, where both results cannot be distinguished).By way of example, we modify the model to increase the probability of infection of individuals in next-neighbors contact with members already infected to p i = 3 4 i(t).As expected, infections grow faster near the onset, the daily maximum happens earlier and results in a larger and narrower peak (while keeping the final total number the same, Fig. 11 blue dotted line compared to thick dotted line).
On the other hand, a scenario where the infection probability is kept constant (p i = 0.1, p r = 0.05) results in a wider and smaller maximum (the infection and recover constant probabilities have been adjusted to yield the peak near the same MC steps on the previous cases, Fig. 11 red dotted line compared to thick dotted line).For these conditions, a typical temporal evolution of individuals (pixels) is shown in Fig. 12.A weak tendency to clustering is observed, although the system is seen to reach a quasi-homogeneous state fast.
Unlike SIR, this model can sustain in a natural way a constant background of infections if at some point in the epidemics p i becomes very similar to p r , establishing a transient regime which we categorize as qualitatively different from the region where the daily distribution derived from SHR or SIR is simply too low.This is a feature that can be observed in real data (Fig. 1 inset in the right-hand side).
Finally, we checked how statistical properties of the model perform and scale under different lattice sizes and parameters via simulation.The distributions over time for N = 100 and N = 1000 are virtually indistinguishable as long as the initial infectious individuals are kept in the same ratio.In order to further visualize the stochasticity under the chosen scale, we show in the right panel of Fig. 11 ten randomly chosen realizations out of one hundred runs with random initial positions of infectious in the lattice.As the starting day where the infection expands is random, we have rigidly displaced the time of the samples such that they peak on the same day.Then, the ten different realizations and their averaged value lie nicely on the same curve and the standard deviation displayed in the inset is seen to be acceptably small.

CONCLUSIONS AND FUTURE VENUES
We have analysed and compared three mathematical approaches of increasing complexity to investigate the dynamics of COVID-19.A take-home message is that all three approaches have enough flexibility to embody the pandemics' actual behaviour for ten arbitrarily chosen countries.However, they display different error bars and have different abilities to be extrapolated into the future to produce valuable predictions We have proved that a least-squares SHR-model based on the lognormal distribution is well suited to describe the epidemic's evolution using only two free parameters per infection wave.Confronted against real data up to the second inflexion point, the values determined for these parameters are well converged and stable, guaranteeing fractional error bars of ±5%.Therefore, the SHR-model is suitable to extrapolate tendencies to the next one or two weeks, even in the presence of noisy data.A simpler averaged version depending only on a single free parameter per wave has been shown to be adequate to be used as a first approximation, albeit with larger associated incertitudes.We have also considered a generalised deterministic SIR dynamics to analyse the temporal evolution of the disease.
In this case, the corresponding two free parameters are well converged and stable once the maximum in the daily distribution of cases is passed, i.e. about a fortnight before the SHR reaches a similar accuracy.Besides the two deterministic models, we have considered stochastic individual-based dynamics reflecting the daily changes in individuals' classes.We examined both the case of uniform and neighbour-dependent transitions via a Monte-Carlo simulation, which has an excellent correspondence with the analogue SIR model's temporal evolution.
While such simple dynamics ignore individual, spatial or further inhomogeneities (e.g., genetic, socioeconomic, or other differences) we have proven that they can reproduce, predict and forecast relevant features of the actual COVID-19 dynamics.In particular, they provide reasonably robust ways to monitor and forecast the actual temporal evolution of contagious diseases in different environments, while only requiring basic mathematical tools.
The analysis of ten different countries makes us conclude that the SHR model can be extrapolated into the future with at most a 5% fractional error after a fortnight passed the second inflexion point.On the other hand, the SIR model, which includes two free parameters only too, seems more stable and can be used with a similar accuracy about one week passed the maximum.Finally, the MC model is helpful to study the interactions between separated regions developing the epidemics.
By comparing SHR and SIR we find an excellent correlation between functions c σ,µ (t) and i τ 0 ,τ 1 (t), and their respective cumulative distributions C(t) and r(t), which suggest that an analytical parametrized solution for SIR might be possible by trying a variational-like approach: Our results strongly suggest that useful bounds can be found for δ(t).Such promising venue will be explored in the future.
On the other hand, the excellent agreement between SIR and MC (provided the transition probabilities are chosen in accordance with the hypothesis behind the SIR model) opens new prospects to whrite spatially resolved SIR-like models that might be solved applying Markov chains techniques.In this appendix, we summarize a few known public-health actions suggested by the SIR model (Weiss, 2013) which are useful to follow the lines of reasoning developped in the paper.
1. Since the maximum value for R(∞) is the entire population, N < ∞, the disease always dies out, I(t > t 0 ) = 0. Otherwise, if for some initial conditions we could have I(∞) = 0, Eq. 3 would imply that R(t) could grow without limit, which proves the fact by reductio ad absurdum.
2. For n = 1 the ratio ℜ 0 = τ 1 τ 0 determines whether the disease grows or dies.Since S(t) can only decrease, we have from the second Eq. in 6 that at the onset, for ℜ 0 >> 1, an initial estimation for τ 0 can be obtained from τ ≈ τ 0 .
4 Therefore, if ℜ 0 < 1, I(t) is a monotonically decreasing function and the infection dies quickly.On the other hand, if ℜ 0 > 1, I(t) increases in the region near t = 0, it reaches a maximum value, I M (t M ), and then it goes to zero, as proved in the point above.ℜ 0 is called the basic reproductive number and it sets up a threshold for the expansion of the disease which is not obvious without a mathematical analysis of the underlying differential equations.
3. The maximum number of infected people can be obtained by dividing the two equations 3 and 6, dS dI = − S(t) i(t) which can be integrated to yield for S(0) ≈ N and I(0) ≈ 1, 4. Similarly, dividing the equation 6 by 4, we get i.e., S(t) = S(0)e −ℜ 0 R(t) , assuming R(0) = 0. Notice that for ℜ 0 >> 1, S(∞) might take the value zero, which corresponds to a very virulent epidemics where everybody dies.

2.
Decrease τ 1 to reduce the duration of infection.Increasing values of τ 1 displace t M towards the origin and increase the value of I M .3. Reduce N = S(0) by vaccination or any other kind of immunity.Increasing S(0) displaces t M towards larger values, decreases i M and r M .. 4. Decrease n.Decreasing n < 1 moves t M towards larger values and, it decreases i M and r M .
MY JN JL AG SP OC NV DC JN FB MR

FIG. 1
FIG.1SHR/Spain.Left/Right panels: deaths/infections related to are taken from(Dong et al., 2020;Roser et al., 2021).Dashed curves fit the data using Eqs 1 and 2. Blue: total accumulated cases per million inhabitant.Red: daily cases per one hundred million inhabitants (the factor ×100 is introduced for the sake of better visibility on the scale of total cases only).The black thin line is a 7-day moving average of data.The green dashed line is the averaged representative curve discussed in section I.A.5.Red and blue thin dotted lines give the contributions of individual waves.The inset (left) gives the percentage between deaths and infections from March to September.The inset (right) enlarges region II where a nearly constant number of infections takes place (red: least-squares fit to data and constant mean value.Black: 7-day moving average).Changes in data collection methodology took place on April 19th, April
MY JN JL AG SP OC NV DC JN FB MR MY JN JL AG SP OC NV DC JN FB MR MY JN JL AG SP OC NV DC JN FB MR MY JN JL AG SP OC NV DC JN FB MR MY JN JL AG SP OC NV DC JN FB MR

FIG. 3
FIG. 3 SHR/Other countries (I).Symbols as in Fig. 1.Changes in methodology took place in United Kingdom (GRB) on May 20th and July 3rd, and in France (FRA) on May 28th.
MY JN JL AG SP OC NV DC JN FB MR MY JN JL AG SP OC NV DC JN FB MR

FIG. 5 FIG. 6
FIG. 5 Comparison of the evolution of number of deaths in NYC and the region of Madrid (CAM) during the first wave.NYC (9.1 M people, red, down pointing triangles, ▽, dashed line) and Madrid (6.7 M people, black, upwards pointing triangles, △, dotted line) duringthe first wave.The data and SHR fits for both locations were juxtaposed matching the day with the maximum number of deaths, aiming to highlight the similarities.The values of the CAM were also scaled to the ratio of population between the two regions ( 9.1 7.6 ) to enable a better comparison.

FIG. 8
FIG. 8 SIR/Spain.Left upper panel: exponential fit near the onset.Right upper panel: initial iteration for deaths (see text).Left lower panel: final iterations for deaths.Right lower panel: final iterations for infections.Blue: accumulated cases, R(t) (per million people).Red: daily cases, I(t)
Parameters for SHR model (confirmed deaths, first wave).P: country's population (millions).µ and σ: parameters in the lognormal distribution.C(∞): asymptotic value for accumulated cases (per million person).R 2 and r 2 : R-squared correlation factors for C(t) and c(t), respectively.t M and t 2 : maximum and second inflection point (origin is the 1st of January, Parameters for SHR model (confirmed infections, first wave).P: country's population (millions).µ and σ: parameters in the lognormal distribution.C(∞): asymptotic value for accumulated cases (per million person).R 2 and r 2 : R-squared correlation factors for C(t) and c(t), respectively.t M andt 2 : maximum and second inflection point (origin is the 1st of January, TABLESTABLE I TABLE III Parameters for SIR model (first wave).N (number of individuals), τ 0 and τ 1 (given in days).Upper: deaths per million people.Lower: infections per million people.
AP MY JN JL AG SP OC NV DC JN FB MR MY JN JL AG SP OC NV DC JN FB MR AP MY JN JL AG SP OC NV DC JN FB MR