The COVID-19 pandemic already has reached practically the whole planet. According to the World Health Organization (WHO, Situation Report-41 [1]), although around 80% of the infected people present mild symptoms (equivalent to the common flu), older people and those with a history of other diseases like diabetes, cardiovascular, and chronic respiratory syndrome can develop serious problems after being infected by the SARS-CoV-2 virus. It was first identified in December 2019 in the City of Wuhan, Province of Hubei, China. From there it spread across Asia, Europe, and the other continents [2]. On March 11th, the WHO declared it to be a pandemic [3].

In this work, we investigate the spread of the pandemic in 4 different scenarios: Germany, Brazil, the Brazilian State of Paraíba, and the City of Campina Grande. We use the data from Germany mostly as a benchmark test for our model, since very likely their COVID-19 data is one of the most accurate in the world, which is due to the widespread testing of its population. Furthermore, the social distancing, isolation, and use of personal protection equipment (PPE) there has been far more efficient than the measures taken in Brazil in containing the spread of COVID-19 infections.

The first confirmed case in Brazil occurred in February 25th of 2020, when a 61-year-old man who had returned from Italy tested positive and the first death due to the pandemic occurred on March 12. It is very unlikely all cases of COVID-19 contaminations evolved from this patient. As reported in other countries, there were multiple imported infections, which usually drive the contagion rate to very high values in the initial stages of the epidemic. In Brazil and in the majority of countries around the world social distancing policies have been adopted in order to decrease the rate of contagion and thus allow that health systems do not collapse and have conditions to treat the gravest cases of the disease.

In the State of Paraíba, according to the Paraíba Department of Health, the first confirmed case of COVID-19 was registered on March 18th of 2020. It was a man who lived in the City of João Pessoa that had returned from a trip to Europe on February 29th. On March 31st the first death due to the pandemic was recorded in the State.

Campina Grande is the second largest city in the State of Paraíba. The social distancing policy was implanted in this city on March 20th of 2020, with the closure of universities, schools, and non-essential stores. The social distancing was implemented in a preventive form since the first case of infection by the SARS-CoV-2 virus only came to be registered 1 week afterwards, on March 27th. The first death due to COVID-19 in this city was only registered on April 19th.

In this work, we adapted the simple and well-known SIR epidemiological model, developed by Kermack and Mckendrick [4], to study the evolution of the dissemination of COVID-19. The SIR model is a well-known and tested epidemiological compartmental model that has been applied to very diverse epidemics (see for example [5–10]). The model receives its name for dividing the population in three groups: susceptible, infected, and removed. Although, SIR and SIR-like models, such as the SIRD and SEIR(D), are simple models compared with far more detailed alternatives, such as proposed by Ndairu et al. [11], Tang et al. [12], or Gatto et al. [13], the simplicity of SIR-like models is a strength when it comes to parameter estimation. More complex models have more parameters to be determined, and, hence, more uncertainties that may render them non-viable for large scale applications to many populations. Usually, the available epidemiological data is not complete enough to provide estimates for all the parameters of more complex models. The model we propose here has three independent parameters that need to be estimated from the epidemiological data: contagion rate, average time duration of infection, and lethality rate. We validated the proposed theoretical model with comparisons of official data of the numbers of confirmed, recovered, death, and active cases due to COVID-19 infections in Germany, Brazil, the Brazilian State of Paraíba, and the City of Campina Grande, located in Paraíba. We used a compartmental model called SIRD, which replaces the removed by the recovered and the deceased cases. Furthermore, in the present work we use a time-dependent contagion rate and time-dependent lethality and recovery rates.

SIR(D) or SEIR(D) models with time-dependent contagion rates are not new [14], even for the COVID-19 pandemic there are already several SIR and SIR-like models with time-dependent parameters. The most difficult part of this approach is the design of a parameter estimation method that is automated and renders the model accurate. The early work by Fang et al. [15] used a SEIR model in which they estimated the parameters based on the epidemiological data using statistical methods, but at the time of publication there were only about 40 data points and scant comparison of data with theoretical model time series. Zhong et al. [16] used a time-dependent SIR model in which they estimate the parameters from the data of active and recovered cases. The authors obtained big fluctuations in the contagion rate and in the removal rate (1/τ). This likely occurred because they approximated the time derivative of the infected (dI/dt) on a daily basis and also, perhaps, because recovered cases data is often less reliable than the confirmed cases data. Chen et al. [17] also used a parameter estimation technique of the contagion and the removal rates, but they did not provide information on the parameters obtained in their results explicitly. Dehning et al. [18] estimated parameters based on a Bayesian inference with the Markov chain Monte Carlo technique, but not many data points were available at the time of publication and no estimates for active cases were provided. Linka et al. [19] also used a Bayesian parameter estimate in a SEIR model. None of the above articles provide information on estimates for death cases. A comparison of models (SIR, SEIR, and a branching point process) highlighting the strengths and difficulties of each model, and stressing the importance of nonpharmaceutical public health interventions was made by Bertozzi et al. [20].

The present work introduces a simple method for estimating the parameters of the epidemiological model dynamical system. The proposed method combines dynamical systems, probability theory, and statistical analysis of data. In addition to the adapted SIRD model used, we develop a probabilistic model to obtain the recovery and the lethality probabilities. Based on this model, we make estimates for the active, death, and recovered cases based only on data of the confirmed cases. The best fit for the active cases provides us with the average time duration of infection, which is held fixed in our model. With the help of this probabilistic model and the epidemiological data, we also obtain the time-dependent parameters of the model: contagion, recovery, and lethality rates. The estimate of the time-dependent contagion rate is based on a fairly simple statistical analysis of the statistical estimate of the active cases. We estimate the time-dependent lethality and recovery rates based on a statistical analysis of the real death cases data. In this way, one avoids the time-consuming task, for the programmer, of obtaining the contagion rate function that leads to the best fit of the data. In this work, the best fit of the data is done automatically in one pass of integration. Thus, this greatly reduces the time taken to fit the available data with the theoretical curve. Furthermore, we present and validate a forecasting method based on Markov chains and on the parameter estimation method we use to make short-term predictions of the evolution of the epidemiological variables. This allows for short-time forecasting from 1 to 2 weeks in advance. In order to measure the error made when fitting the parameters to the real data, we used the Root Mean Square Error (RMSE). More precisely, in this work, we use the normalized RMSE, which allows comparisons between the fittings of different datasets and also comparisons between the fittings of different time intervals of the same dataset. Moreover, we also take into account the fact that the counting of new epidemiological cases is usually corrupted by delays in the reporting and by missed detection of infections, deaths, and recoveries. To account for that, we add to the daily epidemiological data variation zero-mean Gaussian deviates with standard deviation of 5% of the daily number of cases variation. We run this many times to verify if the curve-fitting ability of our theoretical model is robust to such random perturbations.

Additionally, we prove the well posedness of this model (existence, uniqueness, and continuous dependence on initial data). Also, an important theoretical issue we investigate is the dependence of the solutions on the parameters present in the model. Although the continuous dependence of solutions on parameters is a classic result in abstract dynamical systems (see, for example, [21], Theorem 3.), as far as we know, this has not been proved yet for SIR-type models. Hence, we prove the continuous dependence of the solutions on the system parameters. Here, we use techniques similar to those used in da Silva [22] and Diekmann and Getto [23] to accomplish this result.

We hope that this study of epidemiological dynamics be useful in stressing the importance of public health policies about the application, maintenance, and reinforcement of social distancing measures with the objective of avoiding the collapse of the health system.

This work is organized in the following way: In Section 2, we propose our epidemiological model and we prove results on existence and uniqueness of solution and on the continuous dependence of the solution with respect to the initial data and the parameters present in the model. We develop, in Section 3, the probabilistic model with estimates of the probabilities of recovery and death. We also describe how to make statistical estimates for the various parameters used in our epidemiological model. In Section 4, we investigate the evolution of COVID-19 in Germany as a benchmark test for our model. Afterwards, we present our results and discuss about the evolution of the disease in Brazil, in the State of Paraíba, and in the City of Campina Grande. We also present a method of forecasting for short terms the evolution of the pandemic based on Markov chains. We validate the predictions of our theoretical model by fitting the epidemiological data. In Section 5, we use the root mean square error (RMSE) to measure how well our time-dependent compartmental model fits the epidemiological data. In Section 6, we add noise to the epidemiological data to see if our model is robust enough to this sort of perturbation. Finally, in Section 7 we draw our conclusions.

The evolution of the epidemiological model we propose is determined by the following ordinary differential equations (ODE) system

where S(t), I(t), R(t) are, respectively, the normalized susceptible, infected, and recovered populations at time t. We define M(t) as the normalized number of accumulated deaths due to the epidemic. The introduction of M(t) is mostly for convenience, so that we do not have to perform a separate integration from the numerical routine used to integrate the differential equations of our model. The normalization of the variables was achieved by dividing them by P₀, the total initial population of the region considered. For the sake of simplification, we assume that the population is homogeneous such that all the susceptible individuals have the same probability of being contaminated and the infected individuals have the same probability of recovery or death due to the infection. We also suppose that the population evolves in such a way that the newborn babies are all susceptible and the recovered are all immune. The parameter ν is the population birth rate, μ is the death rate before the onset of the pandemic, κ(t) is the contagion rate function, ρ(t) is the recovery rate, and λ(t) is the lethality rate due to the epidemic. Although ρ(t) and λ(t) are changing over time, τ⁻¹ = ρ(t) + λ(t) is held constant.

It is important to point out further differences between the original model (1) and the proposed SIR model [4]. This could become relevant if the pandemic lasts for over a year. This is relevant also as a source of comparison with the death rates due to the COVID-19. One important difference we introduced is that we now allow for time variation in the contagion rate κ(t), which reflects changes in confinement, social distancing, and mask use adopted by the population. In addition, we allow for time variation in the recovered and lethality rates, which might reflect possible improvements in the treatment efficacy of COVID-19 and/or changes in the demographics of the infected. Here, we also allow variations in the total population, taking into account the contributions of birth and death rates to the population evolution.

In this section, we explain the derivation of the reproduction number and develop the estimation procedures for the parameters present in the model. In addition, we describe our forecasting method.

We used the Odeint function of the Python's scientific library package SciPy [30] to integrate the ODE system of Equation (1) with the integration time-step dt = 1.0/24, which corresponds to an hour when the time unit is a day. In the cases investigated, we took τ 14 days. We chose the value that provided the best fit of the active cases with the delay estimate given in Equation (16) when the active cases data was available, such as in the cases of Germany and Brazil. When the active cases data was not available, we chose τ that would give the best fit between the theoretical prediction for the active cases and the corresponding estimate of Equation (16). We suppose this parameter does not change appreciably during the time scale of the outbreak of the pandemic until now, or at least until more efficient treatments are discovered. The initial values used are: S(0) = 1−C₀/P₀, I(0) = A₀/P₀, R(0)=R0/P0, and M(0) = D₀/P₀, with P₀ being the population just before the outbreak of the pandemic, C₀ is the initial number of confirmed cases, A₀ is the initial number of active cases, R0 is the initial number of recovered cases, and D₀ is the initial number of death cases, usually either 0 or 1. We now apply our model to the four cases of COVID-19 dissemination: in Germany, in Brazil, in the State of Paraíba, and in the City of Campina Grande.

In Figure 2, we show results of numerical simulation for a range of values of κ. Unlike the other results, κ is held constant during each time integration of the equations of motion given in Equation (1). This result is important in conveying the message of the paramount importance of the contagion rate on the possible outcomes of the pandemic. Not only we observe an increase in the number of deaths when the contagion rates increase, but we also see a sharp transition. When there is a growth in the contagion rate from 0.1 to 0.15, this gives rise to an extremely sharp increase in the number of deaths. This implies that there is a critical value of κ, around which there is a rapid increase in the total number of deaths due to the pandemic. Beyond the critical value, we see a saturation in the total number of deaths. The value of κ that corresponds to R₀ = 1 in our model is κ = μ+1/τ = 0.0714. We believe, this reinforces the great relevance of social distancing, since increasing the average distance between people, we will be decreasing the rate of contagion and, consequently, decreasing also the number of deaths due to COVID-19.

Although there are several papers that use root-mean-square error (RMSE) and RMSE-like measures of goodness-of-fit such as Fernandes [32], Kırbaş et al. [33], Uba et al. [34], Al-qaness et al. [35], and Gupta et al. [36], the populations are different to the ones we investigated or, as in the case of Germany, the datasets are very short compared with our datasets. The closest reference we found for comparison with our results is the work of Uba et al. [34] that only uses statistical fitting of confirmed cases data for the first 120 days of the Covid-19 epidemic in Brazil. And even then, our model is on par or better than their best fits for this time period. We also found that better measures of goodness-of-fit than RMSE are the normalized RMSE. The normalized RMSE allows comparisons of fits of different datasets and also different periods along the same dataset.

Below we provide results for the normalized RMSE, defined as

where

where N_p is the number of data points, f(t_n) is the estimate from the model, y(t_n) represents the data, and ȳ is the average of y on the same time interval of the measurement of RMSE. We show results for both the confirmed and the death cases. We used two ways to calculate the NRMSE, one is the running window NRMSE with a window length of N_w days. In the cases shown, we use N_w = 90. As a comparison, we also calculated the cumulative NRMSE starting from the first 90 days of pandemic data until a recent date.

In Figure 11, we show two different ways of measuring the NRMSE of our epidemiological data. Both the running-window and the cumulative NRMSE showed very similar results. Overall, the cumulative NRMSE values are slightly smaller for the cumulative case than the running-window case, assuming the window size of the running-window case is the initial length of the cumulative case dataset. We also notice that the cumulative case results are smoother (with broader peaks) than the running-window case. The best fitting was obtained with the data from Germany. We also notice a strong correlation between the Paraíba and Campina Grande data. This is to be unexpected because nearly the same non-farmaceutical measures were taken in both populations. In addition to that, Campina Grande is the second-largest city in Paraíba.

Here, we add noise to the epidemiological data to see if our model estimates are robust enough to the presence of uncertainties in the recording of epidemiological data. These uncertainties can arise from delays in reporting cases and also missed detection due to insufficient testing of the populations. One way to reproduce these uncertainties is to add a Gaussian noise to the epidemiological data. At each day, we add a random value obtained from a normal distribution with zero mean and standard deviation of 5% of the day's cases, which could be confirmed, death, and recovered cases. We add these statistically independent deviates to each one of these datasets, noticing that for the case of confirmed, death, and recovered cases the Gaussian deviate with zero mean and standard deviation of 5% of the day's confirmed, death, and recovered cases, respectively. From these modified epidemiological data, we obtain new time-dependent contagion, lethality, and recovery rates and integrate the equations of motion of the corresponding model. From Equation (22), we estimate that the noise effect on the contagion rate should be small. This is confirmed by our numerical results as one can see in Figure 12A, specifically in the inset.

What is not usually small is the effect on the populations as can be seen in Figure 12B. In this figure, we generated 100 datasets (not plotted) with added Gaussian noise with zero mean and standard deviation on the n-th day given by σ_n = 0.05ΔC(t_n). To each of these datasets we plot the corresponding theoretical estimate (colored solid lines) For comparison, we also plot the confirmed cases data without noise (black dots).

In Figure 12C, we investigate the effect of added Gaussian noise to the curve fitting of the epidemiological data by our model via the cumulative normalized for confirmed cases. Observe that a large fraction of the time-series of the NRMSE have higher values than the the NRMSE without added noise. This likely reflects the fact that the contagion process of COVID-19 is not really random and requires contact between infected and susceptible people as assumed by compartmental epidemiological models such as the one we use. Nevertheless, even in the worst cases, the curve fitting of the data with noise is still very good, since the cumulative NRMSE is small, demonstrating that our model is robust to this sort of perturbation.

Here we summarize the main contributions of our epidemiological study of the evolution of the COVID-19 pandemic in four populations. The proposed model is an adaptation of the SIR model [4], a SIRD model [37], with some notable differences. In the SIRD model we replace the removed population by the recovered and the deceased. We allow that the contagion rate varies in time so that it reflects the fact that social distancing and isolation changes over time. We developed two ways to obtain the active population from the data on confirmed cases only. In the first approach, the number of active cases at time t is estimated simply by the difference between the confirmed population at time t minus the confirmed population at time t − τ, where τ denotes the average time span of infection. In the second approach, the number of active cases is estimated by the probabilistic model proposed here. Both approaches result in fairly accurate predictions of the active cases when the official data is available. In cases in which the data on recovered cases is not available, this approach could be the only way to estimate the number of active cases. Furthermore, we propose a simple method to track automatically the contagion rate based on the statistical estimate of active cases. We divided this data in weekly intervals and for each interval we made a linear regression to obtain the slope, and based on the equation in our model for the infected, we could obtain a daily contagion rate. This contagion rate was fed back into the equations of motion and we were able to fit the data with a minimum of or no ad hoc interventions. We noticed that the contagion rate could reach very high values at the initial stages of the spread of the epidemic in all populations investigated. This likely happens because of inherent numerical errors, as described in Section 3.4.2, and multiple infections that are imported from contaminated visitors or people returning from trips abroad. This type of dynamics is more relevant at the beginning, when the local contaminated population is small and before barriers on traveling were imposed by the governments. Here, we also take into account the contribution from the pre-pandemic birth and death rates to the evolution of the populations investigated. This could become relevant if the pandemic lasts for over a year and also it is important as a comparison for the lethality of the pandemic. According to the results exposed in the previous section, with our model we could fit the official case data from Germany, Brazil, the State of Paraíba, and city of Campina Grande quite well.

The modeling of the spread of the pandemic in Germany is very emblematic, since it clearly shows that the strict social distancing measures imposed by the government on 03/22/2020 were very effective in containing the disease, reducing the R₀(t) from 2.8 down to roughly 1 about 2 weeks later, and further down to approximately 0.5 in more 2 weeks, according to our model. The increase in lethality from December/2020 through March/2021 was probably due to the spread of COVID-19 variants, such as the alpha variant. The decrease in R₀(t) from April/2021 to June/2021 is very likely due to the increased vaccination of the population. The subsequent increase in R₀(t) may be due to relaxing of preventive measures by the population and also due to spread of the delta variant.

In the case of Brazil, we conclude that, based on the fit of the proposed model, the decrease of the reproduction number has been far more difficult and bumpier. This means the implanted social distancing measures are having an effect in thwarting the spread of the disease, but it has not been enough, since to control the pandemic R₀(t) has to remain consistently below 1 for several weeks, whereas in reality it has been oscillating around 1 for most of the time and the number of active cases is still very high. The increase in lethality after March/2021 has likely to do with the spread of gamma variants across the country.

The evolution of COVID-19 in the State of Paraíba has been similar to the national case. The amplitude of the weekly modulations of the contagion rate is considerably higher than in Brazil. This may be due to the smaller population involved. The higher lethality from February/2021 through May/2021 is probably due to the spread of the gamma variants, while the lower lethality after the last week of June/2021 is likely due to the increased anti-COVID-19 vaccination of the population. The decrease in the effective reproduction number after June/2021 may also be due to the vaccination.

As was commented in the Introduction, Campina Grande adopted a social distancing policy 1 week earlier than the report of the first confirmed case. Despite of this, the rate of contagion did not decrease as fast as it did in Germany, Brazil or Paraíba. It presents two major peaks in κ(t) spaced apart by a week since the outbreak of the disease here. Probably these could be linked to agglomeration events such as in-branch governmental relief payments to unemployed people. Also, one sees higher amplitude of the modulation of the contagion rate than in the other populations studied, this might be linked to the smaller size of the population involved. In a similar way to Paraíba, from mid June/2021 onward there is a substantial decrease in the lethality probably caused by the anti-COVID-19 vaccination campaign.

We also proposed and validated a simple forecasting method based on Markov chains and on our parameter estimation method for the evolution of the epidemiological data for up to 2 weeks. For the populations we investigated, the epidemiological data fell well within the 95% confidence interval of our predictions.

Furthermore, we calculated the NRMSE, which measures the fitting error of the differences between the epidemiological data and the estimates obtained from the evolution provided by our compartmental model with the time-dependent parameters. We use the NRMSE since it allows comparisons between curve fittings to data points at different time intervals or even between different datasets. In addition to that, we also analyzed the robusteness of our model when one introduces random errors in the epidemiological data. This mimics the fact that in real life, a fraction of the new cases of confirmed, deceased, and recovered, are reported with delays or are missed entirely.

Based on the results shown here, we conclude that the public health officials should look into the local dynamics of the spread of the disease as they compare with the theoretical predictions of models such as the one developed here. In this way, they will know where the social distancing and isolation is being more efficiently implemented. The models will be more relevant and accurate if there is more testing of the population. Even random testing should be considered, as one gains statistical information on the spread of the disease and discovers where there is more under-notification. Also, cellphone data of the motion of the population, as used by Peixoto et al. [38] and Linka et al. [19], should be considered as a means of predicting the contagion and of identifying hot spots of COVID-19. This could help identify where there are more contagions: bars, churches, supermarkets, offices, pharmacies, hospitals, bakeries, restaurants, public transportation, delivery services, or family visits, etc. Furthermore, by comparing different local strategies one could gain insight on what works better for slowing the spread of the disease. It would be interesting to see the amount of correlation in population mobility and contagion rate, this could be specially relevant at the city level. In a more refined work, one could couple several nearby cities into a network of populations.

The code used in this paper is publicly available at the gitHub repository https://github.com/aabatista/covid19_parameter_estimation. See also the Supplementary Material for a description of the code.

All raw data supporting the conclusions of this article is publicly available. Covid-19 epidemiological data for Germany and Brazil was obtained from the site https://data.humdata.org/dataset/novel-coronavirus-2019-ncov-cases (accessed on 10/09/2021, with data collected until 09/09). The Covid-19 epidemiological data for the State of Paraíba and the City of Campina Grande was obtained from the site: https://brasil.io/dataset/covid19/files/ (accessed on 10/09/2021, with data collected until 09/09).

AB designed research, developed compartmental model, statistical modeling, automated parameter estimation and forecasting, wrote and executed python code, wrote manuscript text, and made figures. SS designed research, developed compartmental model, wrote manuscript text, and wrote mathematical proofs. Both authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.