Modeling the COVID-19 Epidemics in Macedonia

We aim here at modeling the spread of the COVID-19 epidemics in Macedonia over time using the data on the infected individuals. The focus in on the second wave, so we use the data from 27th of May on. For modeling the temporal dynamics of the population of infected individuals, we use a Kermack-McKendrick variant of the SIR epidemiological model. We fit the model parameters to the data using the modeling tool ProBMoT. In turn, we use the model to simulate and predict the dynamics of the infected population in the future. From the predicted dynamics, we can infer the expected time point of the infection peak, the number of infected individuals at the peak and the expected time point of the epidemics demise.

The model

The Kermack-McKendrick variant of the SIR model of the temporal spread of the infectuous disease in a closed, homogeneous population of individuals. It comprises three differential equations modeling the change of the three variables S, I and R, corresponding to the size of the populations of susceptible (S), infected (I) and recovered (R) individuals: \[ \frac{d S}{d t} = - \beta S I \\ \frac{d I}{d t} = \beta S I - \gamma I \\ \frac{d R}{d t} = \gamma I . \] The two constant parameters of \(\beta\) and \(\gamma\) represent the rates of infection and recovery, respectively.

Given the population sizes at the initial time point \(t_0\), \(S(t_0)\), \(I(t_0)\) and \(R(t_0)\), along with the values of \(\beta\) and \(\gamma\), we can simulate the dynamical change of the three populations through time. At any simulation time point \(t\), the sum \(S(t) + I(t) + R(t)\) is constant and equals the size of the whole population.

Fitting the model to the data

The data used for fitting the model comprises the time-series corresponding to the dynamic change of the number of infected individuals in Macedonia from 27th of May to the current date. We regularly update the data using the Covid-19 treker Website.

To obtain the values of the model parameters \(\beta\) and \(\gamma\), we fit the model against the observations of the number of infected individuals. To this end, we use the ProBMoT tool for data- and knowledge-driven modeling of dynamical systems. ProBMoT formulates the fitting of the model parameters to observed data as a task of numerical optimization, where the objective is to minimize the discrepancy between the observations \(I\) and simulated values \(\hat{I}\) of the size of the infected population. The discrepancy is measured using the standard root-mean-square-error (RMSE):

\[ {RMSE} = \sqrt{ \frac{1}{n} \sum _{t = t_0} ^{t_{n-1}} (I(t) - \hat{I}(t))^2}, \] where \(n\) is the number of observations at the time points \(t_0, t_1, \ldots t_{n-1}\). Note that \(RMSE\) is measured using the simulated behavior of the model \(\hat{y}\): ProBMoT simulates the SIR model using the CVODE solver of ordinary differential equations as implemented in the Sundials suite. The parameter fitting of the SIR model leads to a simple two-dimensional numerical optimization problem. To solve it, ProBMoT employs the Differential Evolution algorithm with standard settings and a population size of 100. The obtained solution of the optimization problem assigns values to the fitted parameters that lead to the minimal RMSE discrepancy between the observations and the simulations.

Simulating the model

Once we fit the values of the initial population sizes and the two rate parameters, we can simulate the model beyond the last observation point \(t_{n-1}\) and therefore obtain prediction of the infection dynamics.

Note however, that the model parameters are being inferred from a very limited number of observations: we observe only few data points at the very beginning of the second wave of the epidemics with a (rapidly) growing size of the population of infected individuals. Second, the observations are based on a limited sample of individuals that have been tested. On such a small set of observations, we face a high risk of overfitting. While the bias component of the fitting error can be kept at bay by assuming the SIR model, the variance of the parameter estimates is out of control.

To account for the high variance of the parameter estimates, we perform a bootstrap simulation of the model as opposed to a single model simulation with the optimally fitted parameter values. The bootstrap simulation procedure samples the parameter values in two ways. First, we sample all the parameter values considered by the Differential Evolution algorithm during the optimization: given the algorithm settings, 2,000 evaluations are performed. We constrain the sampling to those parameter values that lead to RMSE, which is within a 5% margin of the RMSE of the optimal parameter values, i.e., lies within the interval \([RMSE_{opt}, RMSE_{opt} \cdot 1.05]\), where \({RMSE}_{opt}\) is the error obtained with the optimal parameter values.

Second, given this initial set of parameter values, we sample their neighborhoods from the Gaussian distribution \(N(p_0, \sigma_p)\), where \(p_0\) and \(\sigma_p\) are estimated from the ProBMoT evaluations selected in the first phase. Again, we constrain the sampling to those parameter values that lead to RMSE, which is less than a half of the RMSE of the simplest model that predicts the mean of the observed values.

For each sample of parameter values, we simulate the model for a period of time from the 27th of May 2020 to the 15th of October 2020. The simulations are aggregated using minimal and maximal simulated values at each time point leading to the graph in Figure 1.

Figure 1: The bootstrap simulation (blue polygon) of the SIR model fitted to the observations of the number of infected individuals in Macedonia (vertical bars) from 27th of May to 15th of October 2020. The solid black line within the polygon depicts the simulation of the model with optimal parameter values (lowest RMSE).

Analysis of the simulations

Using the simulations, we predict the time position of the extreme number of active cases, which we refer to as peak. When considering the peak, both, the magnitude and the time point of the peak are of interest. After observing the robustness of the model predictions during the first ten updates of the model, we put the predictions into perspective in Figures 2 and 3.

Figure 2: Temporal change of the predicted peak time point with the model updates: the black line corresponds to the median predictions, while the red polygon is based on the first and the third quartile of the predictions.

The majority of the bootstraped simulations predict the median time point of the peak to be between end of June and beginning of July. After the initial model updates with unstable and decreasing predictions, the more recent updates lead to a quite stable prediction (with a slight monotonic increase towards the end) of the median peak value at the end of June or 1st of July. The current prediction, based on the model updated with the data up to 10th of July, is 2nd of July: the model predicts that Macedonia is past the peak point of the second wave. The actual peak in the measured data is on 4th of July.

The number of active cases in the next days is expected to decline, which actually happened on 30th of June, the peak being actually observed on 4th of July. For days now, the model prediction of the peak data is overly optimistic. This might be due to the changing dynamics of the number of active cases. While around 20th of June the daily increase of the number of active cases stared to slow down, from 30th of June the daily change varied from -137 (decrease) to 93 (increase). The simple SIR model fails to follow the dynamics accurately.

Figure 3: Temporal change of the predicted peak magnitude with the model updates: the black line corresponds to the median predictions, while the red polygon is based on the first and the third quartile of the predictions.

With the model updates, the median predictions of the magnitude of the peaks vary between 6,000 and 2,800 active cases of COVID-19 cases in Macedonia. The trend of the prediction is almost monotonically decreasing. With the recent updates, the median prediction is quite stable with a moderate increasing trend from 2,800 and 3,600. The current median prediction of the peak magnitude, based on the model updated with the data up to 10th of July, is 3,529 (note that this is a median of the bootstrapped model predictions and is slightly lower than the observed number of active cases).

Important notes

Why this simple variant of SIR (and not, e.g., SEIR or more complex ones)? Because it is simple and has two parameters only, which we have to fit against a very limited amount of data. The variance of the parameter fit is high enough even with these two parameters, fitting complex models with higher number of parameters is out of scope in this situation.

How reliable are these predictions? Please note that this might be considered an optimistic predictions since the observed size of the infected population might be unrealistically small due to the policy of limited, non-systematic testing for Covid-19 in Macedonia. Using the same model in Slovenia (for the first and, so far, the only wave) lead to precise predictions of both, the peak data and the peak magnitude. The good results in Slovenia, however, do not guarantee that the model will show the same utility in the case of modeling the second wave of the epidemics in Macedonia.

What happened on 18th of June? The model could not be properly updated with the data up to 18th of June. The fitted value of the \(\gamma\) parameter was 0, leading to a unrealistic model simulation. That is why we extracted the predictions of the model for 18th of June from the graphs in Figures 2 and 3. This raises additional doubts in the reliability of the model predictions mentioned above.

A model overhaul on 30th of June We realized that the model assumptions related to the range of possible values of the model parameters were not realistic, since the fit of the model to the data started to deteriorate after 20th of June. To obtain plausible fit to the new observations (observations after 20th of June), we extended the parameter ranges from the usual \([0,1]\) to the ones used for fitting the data in Slovenia \([0,5]\). The new ranges lead to models with better fit of the observed data, also for the dates prior to 20th of June. We refitted all the models and revised all the graphs in the report on 30th of June.

Parameter values are out of the ranges prescribed with the SIR model, especially the value of \(\gamma\), which should be at most 1. We assume that the improper parameter values compensate for the measurement bias, due to the limited number of tested individuals and the varying government measures for coping with the epidemics.

---

This is an R Markdown Notebook. The first version of this report was prepared on 31st of March 2020 for the data on COVID-19 epidemics in Slovenia. The last update took place on 11th of July 2020. Modeling and simulation experiments performed by Nikola Simidjievski, Matej Petković and Ljupčo Todorovski.