1. INTRODUCTION
"When will the peak of the pandemic hit? When will it be over?"
Those are unarguably among the most asked questions during the ongoing coronavirus disease 2019 (COVID-19) crisis, i.e., a disease outbreak of atypical pneumonia that originated from Wuhan, China [1]. The disease spread to over 100 countries in a matter of weeks [2], with most internationally imported cases reported to date having history of travel to Wuhan [3], [4]. The pandemic has made governments all over the world take serious responses [5]. The governmental measures result in a significant disruption in the lives of their people that raised such questions above.
Diseases grow rather exponentially at the initial transmissibility of outbreak [6], [7]. When there is no intervention and the proportion of infections starts to become comparable to the entire population, the growth will slow down as susceptible is fuel to diseases. This type of logistic growth will yield the peak of a pandemic
*Corresponding author
Received Mei 21st , 2020, Revised Mei 21st , 2020, Accepted for publication June 12th, 2020. Copyright ©2020 Published by Indonesian Biomathematical Society, e-ISSN: 2549-2896, DOI:10.5614/cbms.2020.3.1.4 that everybody is interested in and its arrival can be forecast using, e.g., the susceptible-infected-removed (SIR) compartment model [8].
However, in the presence of pandemic, human beings adapt. Governments intervene. As such, using the SIR model to predict the peak, while new cases are mainly outcomes of national policies and/or community behaviour, would be similar to forecasting what policymakers would do or the effectiveness of their response [9], [10], [11], which is dynamic and can be unprecedented. On top of that, there is a lack of knowledge of epidemiology characteristics and a high rate of undocumented cases [12]. A brute force analysis by fitting reported data to the SIR model is therefore prone to a false prediction if not done carefully (see, e.g., Fig. 2 of [13] that incorrectly predicted the peak time as well as the total infection of COVID-19 in Italy when compared to the latest data). We will show below how data-driven forecasts are sensitive to the time-series information that we input in the model.
Considering the limitations and obstacles, it is therefore important to determine instead the so-called disease reproduction number or reproductive factor [14], which is the number of people that are infected by one infectious person during an epidemic outbreak [15], [16], [17], [18], [19]. It depends on the duration of the infectious period, the probability of infecting a susceptible individual during one contact, and the number of susceptible people contacted per unit of time.
There are generally two types of such a number, i.e., basic [15] and effective (or instantaneous) [16]. While basic reproduction number is the average expected number of cases directly generated by one case in a population where all individuals are susceptible, effective reproduction number is the number of cases generated in the current state of a population. This paper is intended to give a brief review of these numbers to undergraduate students and a broad science-educated audience in general. We also hope that the paper can be an expository article of epidemiology to policyholders in making public health measures.
To quantify directly the actual reproduction number is difficult, if not impossible, and as such, we can only estimate it indirectly. One common approach is to fit a model to epidemiological data that will provide values of some parameters [6]. Here, we use the SIR compartment model as our model reference, where the reproduction number is associated to the threshold point for stability of the disease free equilibrium.
There are three estimation methods that we will discuss. As a case study, we apply the methods to discuss and forecast viral transmission of COVID-19 in Italy. The first one is by parameter fit to the SIR model [20], which is probably the most popular analysis to the study of COVID-19 [13]. The computed parameters will then be used to obtain the reproduction number. The second method is to use the reported data of infected and removed people [21]. Comparing the number with that obtained using the parameter fit shows a similar trend in the decrease of the infection rate in Italy. The third method is using the ratio of increment of infections from two subsequent days [22], [23]. However, such a quantity is usually highly fluctuating as we demonstrate for the case of Italy. The trend is obtained using, e.g., a parameter fit of the Richards curve [6], [24] to the cumulative cases.
As the methods presented here are all based on the SIR model, they are limited by assumptions commonly made within the SIR model. An important assumption is that the presented data are an accurate representation of what happens in the population, although this can be relaxed for some methods in this paper. Another assumption or limitation is that it does not include people that are infected but not infectious, which can be overcome by incorporating another compartment, such as the commonly used Exposed group.
We conclude the paper with a brief review of improvements to the methods by including randomness (stochastic processes).
2. SIR MODEL AND THE REPRODUCTION NUMBER
The SIR model equations are given by
\[\frac{dS}{dt} = -\beta \frac{SI}{N},\tag{1}\]
\[\frac{dS}{dt} = -\beta \frac{SI}{N}, \qquad (1)\] \[\frac{dI}{dt} = \beta \frac{SI}{N} - \gamma I, \qquad (2)\] \[\frac{dR}{dt} = \gamma I. \qquad (3)\]
\[\frac{dR}{dt} = \gamma I. \tag{3}\]
Here, S and I denote the total number of susceptible and infected individuals. Variable R represents the removed compartment that can consist of recovered (and become-resistant) and deceased individuals. The total population is N=S+I+R. Note that dN/dt=d(S+I+R)/dt=0, which implies that N is constant. The parameters \(\beta\) and \(\gamma\) are the transmission and removal rate constants, respectively. The average length of time an infected individual remains infective, i.e., the infectious time, is \(1/\gamma\). Note that this still applies even when the parameters \(\beta\) and \(\gamma\) are functions of time. Additionally, we denote the cumulative (total) case as T=I+R, which satisfies the equation
\[\frac{dT}{dt} = \beta \frac{SI}{N}.\tag{4}\]
Equation (2) can also be written as
\[\frac{dI}{dt} = \gamma \left( \mathcal{R}_t - 1 \right) I,\tag{5}\]
where
\[\mathcal{R}_t = \frac{S}{N} \mathcal{R}_0, \quad \mathcal{R}_0 = \beta/\gamma.\] (6)
\(\mathcal{R}_t\) is the effective reproduction number and \(\mathcal{R}_0\) is the basic one. Note from (5) that depending on the value of \(\mathcal{R}_t\), i.e., whether \(\mathcal{R}_t > 1\) or \(\mathcal{R}_t < 1\), the infections I will increase or decrease in time, respectively. It is therefore important to track this number to forecast the spread of an infection in an area.
As data are collected and reported regularly in a certain time interval, it is instructive to consider instead the discrete model
\[\Delta S_n = -\tau \beta \frac{S_n I_n}{N},\tag{7}\]
\[\Delta I_n = \tau \beta \frac{S_n I_n}{N} - \tau \gamma I_n, \tag{8}\]
\[\Delta R_n = \tau \gamma I_n,\tag{9}\]
where \(\Delta K_n = K_{n+1} - K_n\), K = S, I, R. \(\tau\) is the time interval, which in the limit \(\tau \to 0\), make the model (7)-(9) become (1)-(3). We take \(\tau = 1\) day, which is the standard time interval to report updates on cases for COVID-19. The reproduction numbers (6) can be checked to be still applicable here. In the following, we will mainly use the discrete SIR model (7)-(9).
As the main part of this report, we will consider three different methods to approximate the effective reproduction number \(\mathcal{R}_t\).
2.1. Method 1: Parameter fit
To calculate \(\mathcal{R}_t\) (6), one needs to determine first the parameters \(\beta\) and \(\gamma\), as well as the number of susceptible \(S_n\) and hence the population size N. Because infection data are given in terms of the number of infected and removed (i.e., recovered or deceased) people, we can find the parameter set for which the model has the best agreement with the data. In that case, we fit the deterministic epidemiological model (7)-(9) by employing a generalized least squares scheme, i.e., we search for the minimum of an unconstrained problem specified by
\[\min_{\{S_1,\beta,\gamma\}} \sum_n \left( I_n - Idata_n \right)^2 + \left( R_n - Rdata_n \right)^2, \tag{10}\] where \(Idata_n\) and \(Rdata_n\) are reported infected and removed cases at day n. Here we only limit ourselves to minimization using three parameters \((S_1, \beta, \text{ and } \gamma)\) only. Note that N is implicitly part of the estimated parameters because \(N = S_1 + I_1 + R_1\). We take \(I_1 = Idata_1\) and \(R_1 = Rdata_1\) at the initial step.
Table 1: Parameters obtained from the minimization procedure in Method 1 for the case of COVID-19 in Italy.
| - | Fitted Day 1-53 | Fitted Day 53-87 | Fitted Day 1-87 |
|---|---|---|---|
| \(\overline{S_1}\) | 77301.031 | 154875.683 | 161479.047 |
| \(\beta\) | 0.282 | 0.119 | 0.246 |
| \(\gamma^{(1)}\) | 0.021 | 0.0169 | 0.017 |
| \(\gamma^{(2)}\) | 0.017 | 0.008 | 0.009 |
Figure 1: Predicted evolution of the COVID-19 outbreak in Italy based on the fits of the discrete SIR model (7)-(9) for (a) active cases; (b) recovered; (c) deceased. Shaded regions represent the official data retrieved from the JHU CSSE repository [25]. There are three different predicted trends, based on the length of the fitted data, see the legends. Day 1 = 31 Jan 2020, Day 39 (national lockdown) = 9 Mar 2020, Day 87 = 26 Apr 2020.
Figure 2: Estimated effective reproduction number using Method 1 (dashed lines) and Method 2 (stars). The lines are necessarily split into two parts following the national lockdown, see the text.
The search is done using fminsearch function of MATLAB that implements the Simplex search method. To illustrate our computation, we consider COVID-19 cases in Italy, which was one of the world's worst-hit countries. Data were retrieved from [25] on 26 April 2020, which in the analysis will be denoted as Day 87 (i.e., Day 1 is 31 Jan 2020). We present in Fig. 1 the reported data and the fitting SIR dynamics. We slightly modify the model by splitting the removed compartment (R) into recovered and deceased ones with rates \(\gamma^{(1)}\) and \(\gamma^{(2)}\) respectively and as such, \(\gamma = \gamma^{(1)} + \gamma^{(2)}\).
From our computations, using all the available data, i.e., Day 1-87, to be fitted into the SIR equations yields rather bad agreement. It turned out to be necessary to split the data into two periods, separated around the national lockdown that was implemented on 9 March 2020 (Day 39). To be precise, it is the threshold date of the lockdown effect that manifested two weeks later, i.e., from Day 53.
The split is to incorporate the intervention and behavioural change of the population in the model that requires the parameter values to vary over time. Note that the SIR model assumes the parameter values to be constant during the fitted period, which is not necessarily correct. Assuming constant value of those parameters implicitly assumes that the decline in active cases is because herd immunity (i.e., substantial decline in susceptible population) has been formed, which has not been detected anywhere, even at places with high death counts. Splitting the graph and fitting the parameters separately are therefore to solve the assumption violation, where an extra care must be taken in the procedure.
Using the splitting, we obtain good agreement as can be seen in Fig. 1. It is important to note that using data from Day 1-53, we obtained a predicted peak at the end of March 2020, which clearly is not correct, i.e., parameter fits depend sensitively on the fitted data. This explains the incorrect prediction of [13].
In Table 1, we list the fitting parameters. Using the values, we plot in Fig. 2 the resulting estimated effective reproduction number \(\mathcal{R}_t\) (6). It is clear that the national lockdown effectively decreased the number. The curve crosses the axis \(\mathcal{R}_t = 1\) at the peak on 19 Apr 2020.
2.2. Method 2: Using infected and recovered data
In the second method, instead of getting the parameters \(\beta\) and \(\gamma\) from fitting, we will derive them from the governing equations (8)-(9) directly. Writing \(\beta \equiv \beta_n\) and \(\gamma \equiv \gamma_n\) on the right hand side of the equations, it is straightforward to obtain
\[\beta_n = \frac{\Delta(I_n + R_n)}{\tau S_n I_n} N, \quad \gamma_n = \frac{\Delta R_n}{\tau I_n}.\] (11)
From the definition (6), we have R0 = βn/γn and as such,
\[\mathcal{R}_t = \frac{S_n}{N} \mathcal{R}_0 = 1 + \frac{\Delta I_n}{\Delta R_n}.\] (12)
We therefore obtain that the effective reproduction number is related to the ratio between the change of the infected and the removed compartments. Because ∆Rn > 0, then Rt < 1 if and only if ∆In < 0.
We show in Fig. 2 the estimated reproduction number using the second method, depicted in stars. Because it uses the increase in infected and removed counts which tend to be highly fluctuating, the resulting curve is also wavering. This could be simply solved by smoothing the data using a moving average filter. Nevertheless, in our case here, we still can observe that it is following the same trend as that obtained from Method 1, i.e., the dashed curve.
2.3. Method 3: Using new cases
The third method is to exploit the daily reported new cases, which in terms of the SIR model will be given by the daily difference of the cumulative cases ∆Tn. Integrating (2) in time between t and (t + τ ) gives us
\[I_{n+1} = I_n e^{\gamma \int_t^{t+\tau} (\mathcal{R}_t - 1)dt} \simeq I_n b(\mathcal{R}_t), \tag{13}\] where
\[b(\mathcal{R}_t) = e^{\gamma \tau(\mathcal{R}_t - 1)}. (14)\]
Here, we denote I(t + τ ) = In+1 and I(t) = In. In the last equation, we have assumed that Rt is constant within the time interval.
On the other hand, we have from (4) a discrete approximation
\[\Delta T_n = \tau \beta S_{n+1} I_{n+1} / N = \tau \beta S_{n+1} I_n b(\mathcal{R}_t) / N \simeq \tau \beta S_n I_n b(\mathcal{R}_t) / N. \tag{15}\]
The last step is expected to be valid for emerging diseases, i.e., I varies slowly.
At the same time, we also have from (4)
\[\Delta T_{n-1} = \tau \beta S_n I_n / N. \tag{16}\]
Combining (15) and (16) gives us the effective reproduction number
\[b(\mathcal{R}_t) = \Delta T_n / \Delta T_{n-1}, \quad \mathcal{R}_t = 1 + \frac{1}{\tau \gamma} \ln \left( b(\mathcal{R}_t) \right). \tag{17}\]
Because b(Rt) is a monotonically increasing function in Rt, it can be enough to plot b itself to determine whether the disease decreases or not.
In Fig. 3(a) we plot b(Rt) from the data of COVID-19 cumulative cases in Italy shown in stars. However, the curve is highly fluctuating that may hide the trend. As mentioned in Method 2 above, this could also be solved by smoothing the data using a moving average filter. Another possible way is by approximating the reported cumulative cases with a continuous function. A natural candidate is certainly the generalised logistic function, also known as Richards' curve,
\[T_n = \frac{A}{\left(B + e^{-C(n - n_0)}\right)^{1/\nu}}. (18)\]
Again using a least square method to fit the reported cumulative cases to the function, we obtain that the best parameter values are A = 119520, B = 0.9858, C = 0.0665, n0 = 0.0086, ν = 0.0232. We plot the fitted data and the approximation in Fig. 3(b) in stars and dashed line, respectively.
Using the approximation (18), we are now able to plot a smooth curve (shown in dashed in Fig. 3(a)) that is expected to indicate the trend of b(Rt) calculated from the reported data.
Figure 3: (a) Plot of b(Rt) in time from the reported data (stars) and the approximation obtained using Richards' curve (dashed line). (b) Cumulative cases from data (stars) and Richards' approximation (dashed line).
3. CONCLUSION
We have presented three simple (or actually simplistic) methods to estimate the reproduction number of the COVID-19 pandemic based on the SIR equations as the underlying model. We applied the methods to the data of COVID-19 cases in Italy, where we saw that the implemented national lockdown had positive impacts that appeared about two weeks later.
To extend the deterministic methods reviewed herein, one may consider complex models that include more compartments [6], [26]. However, to be more realistic, one should include statistical randomness and probability in the calculations.
In the spirit of Method 1, Cintron-Arias et al. [ ´ 20] combined parameter fits with statistical asymptotic theory and sensitivity analysis to give approximate sampling distributions for the estimated parameters. Method 3 has been improved in [22], [23] to include a probabilistic description such that the probabilistic formulation for future cases is equivalent, via Bayes' theorem, to the estimation of the probability distribution for the reproduction number.
In addition to estimating the reproduction number based on a model, it is also possible to approximate the reproduction number from the serial interval (the time between the onset of symptoms in a primary case and the onset of those in secondary cases) without assuming a model [27], [28], [29].
ACKNOWLEDGEMENT
H.S. is extremely grateful to his wife, dr. Nurismawati Machfira, who has happily taken a new additional job as 'head teacher' of their children at home during school closure, while maintaining her job as their primary carer, so that he could still #workfromhome and wrote this paper. Part of this research is funded by Program Pengabdian Masyarakat ITB 2020.