Dynamic data-driven algorithm to predict cumulative COVID-19 infected cases using susceptible-infected-susceptible model

Introduction

The use of epidemiological models to control the spread of disease and predict the course of an outbreak has a long history. In 1760, Daniel Bernoulli proposed a mathematical model for smallpox (Hethcote 2000). At the beginning of the 20th century, William Hamer and Ronald Ross studied the epidemic behavior using the law of mass action (Hamer et al. 1928; Ross 1911). In recent times, the use of epidemiological models is inevitable for better management of infectious diseases (Agarwal and Jhajharia 2021; Gumel et al. 2021; Ramezani, Amirlatifi, and Rahimi 2021).

We have seen the use of various epidemiological models to combat the recent outbreak of Coronavirus disease 2019 (COVID-19). COVID-19 was first reported in Wuhan city of China, but soon spread to other parts of the world (Al-Raeei, El-Daher, and Solieva 2021; Gagliardi et al. 2020). Many authors have used some version of Susceptible-Infected-Recovered (SIR) models to predict the COVID-19 outbreak in different countries or regions (Gounane et al. 2021; Leonardi et al. 2020; Magnoni 2021; Ray et al. 2020; Wangping et al. 2020). The basic SIR model assumes that infected individuals are either recovered (and immune) from the disease or died (Keeling and Rohani 2008). It also assumes the number of deaths from the disease is negligible compared to the total population. However, the World Health Organization (WHO) mentioned that “there is currently no evidence that people who have recovered from COVID-19 and have antibodies are protected from a second infection” (WHO 2020). For example, health authorities in South Korea noticed that 163 patients became COVID-19 positive again after a full recovery (NPR 2020). Brouqui et al. 2021 stated that “in COVID-19, it quickly became apparent that naturally acquired immunity would not, in all cases, provide protection for the months following the first infection”. Iwasaki 2021 claimed that “reinfection cases tell us that we cannot rely on immunity acquired by natural infection to confer herd immunity”. Several other studies have found that individuals who are infected by the COVID-19 may build short-term immunity against the disease, and there is no long-lasting guaranteed protection (Edridge et al. 2020; Liu et al. 2020; Tillett et al. 2021). In this context, when there is no long-term protection from the disease after infection, the Susceptible-Infected-Susceptible (SIS) model is appropriate. In an SIS model, people who recover from the disease are added to the susceptible compartment as they can be infected again. In this work, we consider the SIS model to predict the COVID-19 outbreak.

In an SIS model, the main focus is to determine the number of infected people at a given time point. However, it is also essential for planning purposes to know the cumulative number of infected people at a given time point. One cannot directly find the cumulative number of infected people from the SIS model’s present structure. In this work, our main contribution is to provide an SIS model-structure which can give the cumulative number of infected people easily. We incorporate a death due to disease compartment in the SIS model to estimate the model parameters accurately. We develop a dynamic data-driven algorithm to estimate the model parameters efficiently to predict the cumulative infected cases. In this process, we show how to select the optimal training phase to build the model. Finally, the developed algorithm is implemented using COVID-19 data from Delhi, India’s capital city. We also provide an R-package so that users can easily implement the developed model with their data.

Susceptible-Infected-Susceptible (SIS) model

In an SIS model (Hethcote 1989), there are only two compartments, Susceptible and Infected. An SIS model assumes that an individual has not developed any long-term immunity against the disease after infection and thus is at risk of re-infection; hence, it gets added back to the susceptible population. In other words, as shown in Figure 1, after recovering from an infection, an individual again becomes susceptible. Examples of such infections are the common cold and influenza.

Figure 1:

Pictorial representation of an SIS model, where S stands for susceptible and I for infected.

These equations can well describe an SIS model (Keeling and Rohani 2008),

Here t denotes time. In this work, a day is the smallest unit of time t. However, one can choose other suitable units as necessary. S and I are the susceptible and the infected number of people in the population, respectively. The total population size is N, which is the sum of susceptible (S) and infected (I) populations (see Table 1). From Eqs. (1) and (2), it is evident that d S d t + d I d t = 0 . The parameter β, the transmission rate, is the product of the contact rates among infected, and transmission probability (Bjørnstad 2019; Keeling and Rohani 2008). In other words, the parameter β is the average number of individuals infected per unit time (a day) from an infected person. It can also be interpreted as contact rates between infected and susceptible. Here, by the assumption, I infected individuals can contact some individuals randomly; a fraction of S N of them will be susceptible. The parameter γ is the recovery rate. It is assumed to be 1 T where T is the average duration for which infection lasts in an infected person (Ferguson et al. 2020; Ghosh, Ghosh, and Chakraborty 2020). Equation (1) shows the rate of change in the susceptible (S) population; whereas Eq. (2) depicts the rate of change in the infected (I) population. The term β S I N denotes the number of susceptible people infected daily and is removed from the susceptible compartment and added to the infected compartment. The γI denotes the number of people recovered daily and is added back to the susceptible compartment and removed from the infected compartment. Figure 2 shows a simulated SIS model with S + I = N at all time points. Notice that one cannot get the cumulative infected cases directly from the above SIS model. The death due to the disease is not adjusted into this model. It may affect the efficiency of estimating model parameters when the number of deaths due to disease is not negligible (as observed in COVID-19). In the next section, we consider a modified SIS model to address these two issues.

Figure 2:

A simulated SIS model with the initial values of S, I as S(0) = 1,000 and I(0) = 1, respectively. Here, β = 0.2 and γ = 0.05.

Model equations for modified SIS model

The proposed model (Figure 3) can be well described in these equations,

(3) d S d t = − β S I N + γ I ,

(4) d I d t = + β S I N − γ I − μ I ,

(5) d C d t = β S I N ,

(6) d D d t = μ I ,

(7) with,  d S d t + d I d t + d D d t = 0 ,

where S, I, N, β and γ are the same as defined earlier. The C is the cumulative infected cases from the beginning (see Table 1). It includes every person who is infected or was infected. Note that if a person is infected twice (re-infection), it will be counted as 2 cases rather than 1 in C. The D is the deceased population due to the disease. Note that D does not include death counts from other causes. We assume that the death rate from other causes not involving the concerned disease is the same as the birth rate.

Figure 3:

Pictorial representation of a modified SIS model, where S stands for susceptible, I for infected and D for deceased population due to the disease.

Equation (3) is the same as (1). Equation (4) represents the effective change in the infected compartment. As explained earlier, β S I N number of susceptible individuals get infected daily and are added to this compartment. The μ is the mortality rate of infection. Thus, γI infected individuals are recovered from infection daily, and μI infected cases resulting in death, are removed from this compartment. Equation (5) represents the rate of change in the cumulative infected cases, which is equal to daily infected cases ( β S I N ). Equation (6) represents the rate of change in the deceased compartment which is equal to μI. Note that, d C d t is not included in Eq. (7) as the addition of β S I N is already done using d I d t . In other words, d C d t is the first term of d I d t .

Dynamic data-driven algorithm to estimate model parameters

In general, the SIS model parameters are constant for the entire duration of the study period. When the disease under consideration is present in the community for a longer time, the estimated parameter based on the entire study period may not give the right picture. For example, the COVID-19 disease outspread is highly unpredictable in the long term because contact rates and transmission probabilities are changing over time. They vary due to various reasons like control measures implemented by respective governments. Therefore, it may be appropriate to train an SIS model with a shorter training phase and make short-term predictions. Here, the “dynamic data-driven algorithm” means the training phase, used to estimate the model parameters, is dynamic (not fixed) and optimally chosen based on the appropriate historical data.

The two phases of the study period are the training and the prediction phases. Figure 4 shows how the study period is divided into different parts for estimation purpose. We define the four-time variables as follows:

1. T _Current: Denotes the date when the training phase ends. After this date, the prediction phase starts.

2. T _Start: Denotes the length of the default training phase (in days). The default training phase is the interval, [T _Current − T _Start + 1, T _Current].

3. T _Pred: Denotes the length of the prediction phase (in days). The prediction phase is the interval, [T _Current + 1, T _Current + T _Pred].

4. T _Limit: Denotes the upper limit of the number of additional days that can be added to the default training phase to optimally choose the training phase. The length of the training phase keeps increasing with a step of 1 day. Therefore, the maximum training phase interval can be, [T _Current − T _Start − T _Limit + 1, T _Current].

where T _Current − T _Start + 1 ≤ i ≤ T _Current; C _o [i] denotes the observed cumulative infected cases on ith day and C _p [i; t, β, μ] denotes the predicted cumulative infected cases on ith day considering [T _Current − T _Start + 1 − t, T _Current] as the training phase. The optimal value of t is obtained (using (8)) as

Finally, we obtain the optimal training phase that can be used for future prediction as [T _Current − T _Start + 1 − t _opt, T _Current].

Figure 4:

Different time points of training and prediction phases.

There are three parameters in the modified SIS model, namely, β, γ, and μ. As argued earlier, the recovery rate, γ = 1 T , where T is the average duration for which infection lasts in an infected person. In other words, T refers to the average duration of recovery. In the case of COVID-19, T is taken as 14 (Ferguson et al. 2020; Ghosh, Ghosh, and Chakraborty 2020). However, in some specific cases, T can be more than 14 also (Arifin et al. 2020). Given a training phase, Algorithm 1 is used to find the estimates of β and μ by minimizing L(t, β, μ).

Algorithm 1:

Dynamic data-driven algorithm to estimate β and μ.

Algorithm 2:

Prediction algorithm for the cumulative infected cases (C).

Prediction of cumulative infected cases

The β ̂ opt and μ ̂ opt are the optimum values of β and μ, respectively, using Algorithm 1. Using Algorithm 2, predicted values of cumulative infected cases (C) are obtained for every day starting from T _Current+1 to T Current + T Pred . Here, the prediction period’s length depends on the user-supplied value of T _Pred. For 1 ≤ i ≤ T _Pred, the root mean squared error for prediction is

Equation (8) in Section 3 refers to the root mean squared error for the optimal criterion to choose the appropriate training phase. Whereas Eq. (10) refers to the root mean squared error for prediction. Denominators and parameter values in these two equations are different. Note that the equation of cumulative cases in the Algorithm 2 is approximated.

Application on real data

Predicting cumulative infected COVID-19 cases for Delhi

We consider the COVID-19 data from Delhi, India’s capital city with a population size of around 20 million, to demonstrate the proposed algorithm’s prediction performance. Delhi observed more than 600 thousands of cumulative COVID-19 infected cases at the end of 2020. The data is publicly available from https://www.covid19india.org/.

In Figure 5, we have considered four different T _Current as 29 May 2020 (in (A)), 24 July 2020 (in (B)), 29 December 2020 (in (C)) and 15 January 2021 (in (D)). This set-up can check the proposed algorithm’s prediction performance using the modified SIS model concerning different time periods. Table 2 shows the β ̂ opt , μ ̂ opt , root mean squared error (RMSE) for prediction and other related information. From all the four graphs in Figure 5, it is evident that the proposed algorithm is working well to predict the cumulative infected cases with two different prediction periods 30 (for (A) and (B)) and 40 (for (C) and (D)). From Table 2, we see that the chosen optimal training periods’ lengths can be different with different values of β ̂ opt , μ ̂ opt . Notice that β ̂ opt is decreasing over time from 0.12 to 0.09 for Delhi, whereas μ ̂ opt increases for the first three scenarios (from 0.007 to 0.091) then dropped a little to 0.070. The reported RMSEs, specially in the last two scenarios, are quite good considering the scale of cumulative cases above 600,000. These observations support the idea of estimating the modified SIS model’s parameters based on the optimal training phase instead of the entire history as the training phase.

Figure 5:

Modified SIS model predictions based on optimal training phases for Delhi at different periods.

Table 2:

Summary of the prediction process for Delhi.

Scenarios of T Current	β ̂ opt	μ ̂ opt	T Current	Assessment	Optimal training	Prediction	Root mean
			(YYYY-MM-DD)	period (days)	length (days)	length (days)	squared error
A	0.12	0.007	2020-05-29	15	36	30	1,639.27
B	0.12	0.086	2020-07-24	15	31	30	2,678.58
C	0.11	0.091	2020-12-29	30	31	40	245.43
D	0.09	0.070	2021-01-15	30	34	40	976.25

Figure 6 shows what could happen if we include the entire history as a training phase to estimate the model parameters. The blue line is the fitted trained line based on the estimated parameters from Algorithms 1 and 2. The red line indicates the fitted trained line based on the estimated parameters from the entire history data, not using Algorithms 1 and 2. The 30-day prediction curve based on the entire history (125 days) is exponentially higher than the observed curve of the cumulative infected cases (root mean squared error=175,884.1). The difference between the two curves is getting much bigger for the latter part of the prediction period. However, the prediction curve based on the optimal training phase (a total 23 days with 15 days of the assessment phase) is closer to the curve of observed cumulative infected cases (root mean squared error=3,090.25). The estimated model parameter β ̂ = 0.24 , μ ̂ = 0.09 based on the entire history, whereas the same parameter estimates are 0.10 and 0.06 using the optimal training phase, respectively. The β ̂ is quite higher in the case of entire history compared to the same ( β ̂ opt ) based on the optimal training period. It suggests that the estimation of β should be based on an optimal training period to capture the most recent trend rather than the overall trend using entire history.

Figure 6:

Comparison of the modified SIS model predictions for Delhi based on the optimal training phase and the entire history period. Based on the optimal training period (fitted by the blue line), the green line shows a better prediction than the purple line’s prediction using the entire history as the training period (fitted by the red line).

Importance of the deceased compartment

Incorporating the deceased compartment into the modified SIS model is crucial because death due to disease may not be negligible. For example, in COVID-19, the number of deaths to the number of people infected is significant in many countries. Figure 7 shows the importance of the deceased compartment in the modified SIS model in terms of μ for Delhi. The prediction curve (purple line) with β ̂ opt = 0.08 and pre-fixed μ = 0 (no deaths due to disease) shows a large disparity from the observed cumulative infected cases, and the difference between the two curves keeps increasing over time, with root mean squared error 12,990.46. Here, the prediction is based on the optimal training period denoted by blue line. The prediction curve (green line) with β ̂ opt = 0.11 and μ ̂ opt = 0.071 is much closer to the observed cumulative infected cases’ curve with root mean squared error 5,187.47. Here, the prediction is based on the optimal training period denoted by red line. In both scenarios, the prediction phase and the assessment phase are of 40 and 15 days, respectively.

Figure 7:

Shows the importance of the deceased compartment (μ) in the modified SIS model for better prediction performance. The red line represents the optimal training period with optimal estimates of parameters β and μ. The blue line represents the optimal training period with an optimal estimate of parameter β but μ = 0, fixed.

Discussion

This work has provided a modified SIS model that accounts for deaths due to disease and predicts cumulative infected cases based on an optimally chosen training phase. The estimation process described in this work is beneficial when the disease under study changes its spreading pattern over time. We have developed the modified SIS model considering COVID-19 as the disease under focus. However, the model and algorithms can be applied to predict the cumulative cases of other infectious diseases.

Even though one can predict for any period-length in the future using the developed model, we recommend restricting the prediction to the short-term only. Any prediction with more than 30 days may not be reliable due to continuous changes in the COVID-19 virus’ characteristics and human behavior (e.g., how social distancing norms followed from time to time). For example, in Delhi, considering the current day as 24 July 20 with 30 days assessment period, the root mean square errors (RMSE) were 2,166.65, 6,046.05, 16,305.69, and 32,412.88, corresponding to the number of prediction days 30, 40, 50, and 60 from the current day, respectively. Similarly, with the same setup, for the current day, 29 December 2020, the root mean square errors (RMSE) were 159.16, 245.43, 468.84, and 896.93, respectively. Therefore, an increase of 30 days of prediction phase from 30 to 60 days can increase the RMSE substantially (almost 15 times and 6 times in two examples, respectively).

Note that the dynamic data estimation of the parameters, using Algorithms 1 and 2, is approximated using only one variable. Therefore, computational errors may occur when there are more sources of variation in active cases are present in reality. In this work, we have taken a fixed value of γ = 1/14. However, some studies also reported other values of γ (Arifin et al. 2020). For the current day, 24 July 2020, we have applied Algorithm 1 and 2 considering different values of γ as 1 5 , 1 7 , 1 10 , 1 12 , 1 14 , 1 15 , 1 20 , with fixed 30 days assessment period and a prediction length of 30 days. The RMSE values are 14,157.61, 4,134.77, 4,161.06, 8,208.95, 2,678.59, 4,119.57, 4,167.57, respectively for different γ values. In this scenario, the lowest RMSE value is observed when γ = 1 14 . However, in general, this might not be the case.

The developed open-access R-package (https://github.com/abh2k/sisd) can be helpful to implement the modified SIS model without dealing with mathematical details of the model. One only needs to prepare the input data set as described in the R-package documentation.

The objective of infectious disease prediction is to give the respective Governments an idea of what can happen in the near future (say 30 days) so that they can act promptly to avoid more difficult situations. Depending on the Government approach and the participation of the common people in the next 30 days, the accuracy of the predicted numbers may vary. For example, suppose we predict an increase of 100,000 cases in the next 30 days, and the Government imposes a complete lockdown from tomorrow. In that case, no model can be able to predict accurately based on the history data.