1.Introduction

The ongoing COVID-19 epidemic represents a global challenge for health-care systems. The SARS-COV-2 coronavirus causing the COVID-19 is highly transmissible ^1,2. Nevertheless, the epidemic propagation in a region/country is influenced by several factors, including governmental measures ³. The Cuban government has implemented several strategies of social isolation to control this epidemic. As a result, a small number of cases is reported in the province of Santiago de Cuba. Mathematical models are useful for understanding the COVID-19 transmission and the influence of government interventions ^4-7. Currie et al. ⁴ suggest that system dynamics (differential-equation-based modeling) is suitable to model government decisions that cause social distancing and isolation.

The most common system dynamics in epidemic modeling are SIR (susceptible - infectious - removed) and SEIR (susceptible-exposed-infectious-removed) type. The classic SIR and SEIR models assume that individuals change classes (e.g., from susceptible to infected) at constant rates. This assumption is not realistic in most situations ⁸. Additionally, it does not explain why the effective reproduction number of COVID-19 changes over time ⁹. Some authors introduce time-dependent parameters to formulate realistic models of the current pandemic ^{5 10-12}.

A challenge for COVID-19 modelers is to describe the difference between total cases and diagnosed cases. One approach is to exclude this difference and fit the model to reported data ^16-8. Another option is to consider undetected cases as a compartment of the differential equations system ^6,7,12,19,20. Finally, some authors avoid the complication of adding a new compartment, and they model the error with different ad hoc methods ^5,10,21,22.

Mathematical models have been used to forecast the dynamics of COVID-19 in delimited geographical areas. These models usually fit big-number data after several days of epidemic evolution ^7,10,23-26. The official data reported from Santiago de Cuba show a maximum value of 36 active cases reached 31 days after the first confirmed case. This peak is followed by a decreasing trend. These small and noisy data make difficult the model fitting. A common strategy to describe COVID-19 dynamics is to update one or more model parameters based on key days of government actions ^7,10,12.

Ordinary differential equation (ODE) based models are widely used to capture the kinetic of biological systems and to predict the behaviour of time-dependent variables. The minimization of a cost function (as a strategy for fitting a given dynamic model) does not take into account the uncertainty of parameters. Additionally, this approach can lead to unrealistic fittings ²⁷.

Bayesian inference offers a platform that deals naturally with both problems, but it requires the calculation of likelihood functions ^13-15,28. In many practical problems, likelihood functions become analytically intractable; nevertheless, data can be simulated from a parameter vector by applying some algorithm that depends on the model. Thus, several techniques have been developed to infer parameters without using likelihood functions. These techniques are usually called approximate Bayesian computation (ABC) or free likelihood inference ^28,29.

ABC can be used to estimate parameters when the available data of an infectious disease is coarse ³⁰. Besides, ABC accomplishes parameters estimation for models that include dedicated laws to describe the error between total cases and reported cases.

The aim of this paper is to adapt two models reported in the literature to capture the dynamics of the noisy and small-number data of COVID-19 in Santiago de Cuba. For this, we modify an extended SEIR model (10) by introducing a dynamic-mode gamma distribution to describe the error between total and reported cases. Additionally, we adapt the classic SIR model reported in ³¹ by introducing Poisson noise in the observation of removed cases. For both models, we update a parameter on critical days of government actions and estimate the parameter values with ABC. We are not aware that these modifications have been reported in the literature.

2.Methods

3.Results

Model-I simulation (Fig. 1) was performed with MAPs calculated with ABC SMC. All parameters and results of model-I fitting are summarized in Table I. The Mann-Whitney U test performed on the residues of model-I fit did not reject the null hypothesis based on a 5% significance level.

Figure 1 Estimated evolution of COVID-19 epidemic in Santiago de Cuba with model-I. The selected parameter values for simulation correspond to the MAP of the approximate posterior distribution. The model-predicted and reported infected individuals a), and the cumulative b), both predicted in simulation (C _T and C _R predicted) and observed (C _R ) cases are plotted. The shaded area is limited by the extremes of 95% credible intervals of posterior distributions obtained with the ABC SMC sampler. 

Figure 2 was obtained using the MAPs of posterior distributions of model-II parameters that were computed with ABC SMC. These MAPs and other values from model-II are summarized in Table II.

Figure 2 Estimated evolution of COVID-19 epidemic in Santiago de Cuba with model-II. Predicted infected individuals follow the dynamics of the reported ones a). A long-term scenario of both reported and predicted cumulative cases, and the area limited by the extremes of the 95% credible intervals of posterior distributions obtained with the ABC SMC sampler b) are plotted. 

In all trials made, the ABC H sampler consumed less time and reached similar parameter regions to the ones obtained by ABC SMC (Figs. 3-6) for the same prior distributions.

Model-I predicted a total of 73 infected at the end of the COVID-19 epidemic in Santiago de Cuba, with 21% of undocumented cases. The final value of M(t) was 57, being the model-I prediction of confirmed cases. Model-II forecasted 51 confirmed cases at the end of the epidemic. Additionally, the ML ratio calculated was 0.941; and the values of determination coefficients were R² = 0.9450 for model-I and R² = 0.9781 for model-II.

4.Discussion

COVID-19 dynamics in Santiago de Cuba is difficult to capture empirically without the introduction of stepwise functions for parameter 𝛼 in model-I and parameter β in model-II. Results show that model-I (Fig. 1) and model-II (Fig. 2), fitted with ABC algorithms, can predict an early small-number peak of infected cases and the approximate ending of the COVID-19 in Santiago de Cuba. The non-rejection of the null hypothesis on the residuals allows us to continue considering the validity of the error law for model-I.

All MAPs calculated with ABC H are within credible intervals obtained with ABC SMC, except parameter λ of model-I. Even for this parameter, the histogram generated by ABC H partially overlapped the one generated by ABC SMC. We deduce that both algorithms find the same parameter region for the chosen prior. A visual comparison of Fig. 3 with Fig. 4, and Fig. 5 with Fig. 6, reveals that ABC H does not sample from the same distribution as ABC SMC. Although this does not affect the use of ABC H in this work, a modification of the ABC H sampler is required to make it a fully functional ABC algorithm. The scheme w={1,0,0,0} for sampling choice corresponds to a ABC RS process ^33,40. The scheme w={0,0,0,1} resembles an ABC MCMC method ⁴¹. ABC RS and ABC MCMC are algorithms that sample from adequate posterior distributions but in an inefficient time, so other intermediate schemes may be evaluated in further works.

Figure 3 Posterior density estimation of model-I parameters using ABC SMC sampler and COVID-19 data from Santiago de Cuba. The estimated parameters are β₀, k, λ, α₁, α₂, α₃, γ, E(0), and the scale parameter of the gamma distribution (α _GAM ). 

Figure 4 Posterior density estimation of model-I parameters using ABC H sampler and COVID-19 data from Santiago de Cuba. The estimated parameters are β₀, k, λ, α₁, α₂, α₃, γ, E(0), and the scale parameter of the gamma distribution α _GAM . 

Figure 5 Posterior density estimation of model-II parameters using ABC SMC sampler and COVID-19 data from Santiago de Cuba. The estimated parameters are β₁, β₂, β₃, γ, and the interval expectation of the Poisson distribution (λ_Poisson). 

Figure 6 Posterior density estimation of model-II parameters using ABC H sampler and COVID-19 data from Santiago de Cuba. The estimated parameters are β₁, β₂, β₃, γ, and the interval expectation of the Poisson distribution (λ_Poisson). 

The computed percentage of undocumented cases is in contrast with the 86% of unreported infected cases estimated for China ⁴³. This result suggests that the limited progression of the COVID-19 epidemic in Santiago de Cuba may be due to opportune epidemiological investigations, effective control measures in each source of infection, and a low number of initially exposed individuals (E(0) = 6). Additionally, case clusters prevail over transmission clusters for COVID-19 in Santiago de Cuba province. The policymakers used these estimates to complement an investigation into the initial infectious load of SARS-COV-2 in Santiago de Cuba.

The MAPs of a posterior distributions were estimated at α = 0.5245, α₂ = 0.2855, and α₃ = 0.9665 (Table I). Value of 𝛼 is expected to grow as government measures increase but, surprisingly, the value after the closure of transportation between provinces is lower than before (α₂ < α₁). We consider two possible explanations for this result:

First, the interpretation of α as government action is somewhat forced. In a previous study ⁴², α was interpreted as seasonality of transmission associated with the school calendar for a similar formulation of β(t). Consequently, the interpretation of α remains open.

Second, the increase in PCR tests coincides with the period corresponding to α₂. Neither the deterministic model nor the error law explicitly takes into account the variability of PCR tests. This variability influences the number of detected cases, and therefore the estimation of all transmission-rate parameters.

Thus, in our results, we do not interpret α as government action. We consider 𝛼 as a multifactorial parameter that works as a calibrator of β(t) dynamics. The only dynamic variable the explicitly influences β(t) in Eq. (2) is 𝐷 ^{10 42}. This influence is in agreement with the epidemiological and sociological researches carried out in Santiago de Cuba, which confirm the reduction of the transmission rate by the high-risk perception of both decision-makers and population. Nevertheless, D depends explicitly on I, which in turn depends on E; so the variable D is affected by changes in these state variables; and therefore, they indirectly influence the dynamics of β(t).

The estimated value k = 3430 (k = 1117 in ¹⁰) may indicate a good response from individuals in Santiago de Cuba to COVID-19 epidemic.

The mean infection period for model-I may be explained by the inclusion of unreported cases and the delay in the reports of those recovered. The value of this parameter for model-II is plausible because it is the number of days expected for a diagnosed patient to recover, and the diagnosis is usually confirmed several days after the exposition. The interval expectation of removed, λPoisson=18 days with CI [12,23], seems consistent with the patient-discharge policy implemented by the Ministry of Public Health of the Republic of Cuba.

One possible explanation for the ML ratio is that the error law for model-I (Eq. (3)) can only partially describe the complex variability of the reporting rate. As model-II is more parsimonious than model-I, we did not perform some further analysis like Akaike’s information criterion (AIC) or Bayes factor, which should favour model-II over model-I. On the other hand, the model-I allowed us to estimate the number of individuals initially exposed, and to compute the percentage of unreported cases. We argue that simpler models like model-II should be preferred for accurate prediction, whereas more explanatory models like the model-I are to be used for phenomenological analysis.

The chat on the left in Fig. 2 shows a remarkable tracking of the data variation by the fitted curve. A conventional SIR model cannot achieve this. Model-II accomplishes this tracking by dividing the adjustable parameters of the SIR model into three time periods. Furthermore, the cases recovered on the day estimated by the SIR model are not subtracted, but on the day expected by the Poisson error. In other words, model-II is advantageous over a conventional SIR as it is a piecewise SIR with a Poisson noise in the number of removed cases.

We hold that the difference between the predicted and reported data in both models is acceptable because of reporting delays, which are in agreement with previous studies ^10,32,44.

One advantage of using a Bayesian approach in parameter estimation is to use different prior distributions in order to set favourable and unfavourable scenarios for models. When the data do not show a peak yet, most models become highly unidentifiable, in such a way that one can find several parameter combinations leading to reasonable fits. Therefore, we deem it essential to complement the mathematical analysis with an evaluation of epidemiological investigation, looking for the most plausible prediction.

The error law for the model-I (Eq. (3)) and the assumption of the gamma distribution cannot fully capture the variability of the reporting rate. This approach is a baseline for future improvements.

The influences of β₀, 𝛼, 𝐷 and 𝑘 in β(t)should not be discussed separately. The variables β₀, a, D, and k are dynamically interrelated, taking into account that the society (complex system) is constituted by the interaction and coupling of its elements as a global unit, in agreement with ⁴⁵.

Further studies will evaluate the influence of other variables on the favourable evolution of COVID-19 in Santiago de Cuba: government-induced social isolation, the social response of individuals, and the environmental factors (e.g., temperature and relative humidity). Additionally, future works should study the ABC H performance for intermediate values of w, looking for an acceptable compromise of efficiency and accuracy. Another possibility is to study adaptive schemes of w, in order to accept a group of candidates in a short time, and then ensure sampling of the posterior distribution.

5.Conclusions

Our study indicates that modified SIR and SEIR models, combined with ABC for parameter estimation, can follow the small-number dynamics of the COVID-19 epidemic. More precisely, models follow dynamics with no explosions of new cases and with small-number peaks by including some parameters as stepwise functions of critical days. An opportune epidemiological investigation, along with a low number of initially exposed individuals, may partly explain the favourable evolution of the COVID-19 epidemic in Santiago de Cuba.