The daily computed weighted averaging basic reproduction number R 0 , k , ω n for MERS-CoV in South Korea.

In this paper, we propose the daily computed weighted averaging basic reproduction number R 0 , k , ω n for Middle East respiratory syndrome coronavirus (MERS-CoV) outbreak in South Korea, May to July 2015. We use an SIR model with piecewise constant parameters β (contact rate) and γ (removed rate). We use the explicit Euler's method for the solution of the SIR model and a nonlinear least-square fitting procedure for finding the best parameters. In R 0 , k , ω n , the parameters n , k , and w denote days from a reference date, the number of days in averaging, and a weighting factor, respectively. We perform a series of numerical experiments and compare the results with the real-world data. In particular, using the predicted reproduction number based on the previous two consecutive reproduction numbers, we can predict the future behavior of the reproduction number.

Introduction

Middle East respiratory syndrome coronavirus (MERS-CoV), first identified in Saudi Arabia in June 2012, is a viral respiratory disease caused by a novel coronavirus  [1]. South Korea experienced the largest outbreak of MERS-CoV infections outside the Arabian peninsula  [2]. There have been 186 confirmed infective cases in South Korea within two months after the first infective person who returned from a trip to the Arabian peninsula was diagnosed on 11 May, 2015, and 38 people died and more than 16,000 people were quarantined due to the spread of the disease. Fig. 1 represents the epidemic curve of MERS-CoV, South Korea, 20 May–17 July 2015  [3], [4].

Fig. 1

Epidemic curve of MERS-CoV, South Korea, 20 May–17 July 2015.

In order to implement appropriate surveillance and control measures, it is very important to predict the future trend of the epidemic. Therefore, it is the purpose of the present paper to show daily computed reproduction numbers for epidemics, in particular, MERS-CoV, so that we can predict the future behavior of the reproduction number day-by-day.

In the mathematical modeling of the spread of an epidemic disease, it is crucial to estimate the parameters (e.g., contact rate and recovery rate in an SIR model), but it is difficult to do so due to the lack of data available for the estimation. In case of the MERS-CoV spread in South Korea, there are almost complete daily data of the spread of the disease, and we choose the MERS-CoV spread case in South Korea as an epidemic model to verify the validity and accuracy of the method proposed in this paper. To the authors’ knowledge, there have been no previous works that present computational methods for estimation of the values of the epidemic parameters based on actual daily epidemic data.

The rest of the paper is organized as follows. In Section  2, we describe the mathematical model. In Section  3, we provide a numerical algorithm for the estimation of the parameters. We perform several numerical experiments in Section  4. In Section  5, we provide a summary and present our conclusions.

The mathematical model

In this paper, we consider the SIR model that was introduced in 1927 by A.G. McKendrick and W.O. Kermack  [5]. The model is simple and has been widely used so far, for instance, in multigroup epidemic modeling  [6], online social network dynamics  [7], the model adopting the delay term  [8], stochastic model  [6], [9], [10], vaccination strategy  [11], [12]. In this model, the population is divided into susceptible , infected , and removed individuals at time . The governing ordinary differential equations for the SIR model are as follows: with initial condition , , and . Here, the parameters and denote the contact (susceptibility to disease) and removed (either dead or recovered) rates from disease, respectively  [13]. We assume that a removed individual can never be infected again. Let be the constant total population size. Therefore, it is sufficient to solve only two Eqs. (1), (2), i.e., . By assuming that and are time-dependent parameters, we generalize Eqs. (1), (2), (3) as follows: For example, the transmission of the vector-borne diseases such as Dengue fever and Malaria has strong seasonal patterns and thus the parameters are estimated as the time-dependent functions  [14], [15]. Another example is a seasonal SIR model for the spread of the whooping cough in which the contact rate is given as a time-dependent periodic function which accounts seasonal changes  [16].

The numerical method

Let us assume that we have daily statistics about , , and . Let be time, where is a reference time. Let , , and be observed susceptible, infected, and removed data at time , respectively. Let , , and be numerical approximations of , , and , respectively. In this paper, we propose a computation of and on a daily basis, which fit best to the real data. To be more specific, for given data and , we want to find piecewise constant parameters and which minimize the following where , , and are the numerical solutions of Eqs. (4), (5), (6) with initial condition . We divide one day between and into subintervals, then is the time step size. Let and , then by applying the explicit Euler’s method to Eqs. (4), (5), we have where we have used the condition . Unless otherwise stated, we use in all numerical tests. If the solutions of Eqs. (8), (9), (10) are obtained for all , then we define , , and . Note that while we can use higher-order numerical methods such as Runge–Kutta schemes to solve Eqs. (4), (5), we use the explicit Euler’s method for the sake of simplicity. Next, to find the best and , we use a MATLAB routine called lsqcurvefit, based on the least-squares sense, which finds coefficients and that solve the problem (7). We use and as an initial guess in all numerical simulations. The MATLAB code is given inAppendix for the interested readers.

Numerical simulations

In all computations, the initial conditions are taken to be , , , and with May 2015. Fig. 2 shows the time series of observed MERS-CoV epidemiological data and numerical solutions , , and , respectively.

Fig. 2

Time series of observed MERS-CoV epidemiological data and numerical solutions of (a) , (b) , and (c) , respectively. Here, is used.

The basic reproduction number is the average number of secondary infections caused by one infectious individual in a completely susceptible population. It is an important measurement in that one can determine the stability of the disease-free equilibrium by computing the value of ; implies that the disease-free equilibrium is locally asymptotically stable  [17], and if , it is unstable  [18], [19]. Finding basic reduction number can be done through the next generation matrix  [20], [21] whose spectral radius is defined as . In the SIR model, the basic reproduction number is computed as   [22], [23], which will be used for our method in the next section.

Weighted averaging basic reproduction number

In this section, we propose the weighted average values for infected and removed rates, and then define weighted averaging basic reproduction number. Let for . We define -days () weighted average as where and the superscript for represents an exponent. Also, is defined similarly as Eq. (11). Finally, we define the daily computed weighted averaging basic reproduction number as . Fig. 3 (a)–(c) represent , , and for , respectively, with , 0.3, and 1. Fig. 4 shows the results of the parameters for . Note that when , we have for all . If is large, then we have more smooth profile and the transition point () moves to the right. If is close to zero, then we have a similar profile to . Note that we have extraordinarily large values of from to approximately since the recovered individuals were almost zero () in early times and the infected ones increased rapidly.

Fig. 3

Weighted averaging , , and for . In each case, we compare the values with , 0.3, and 1.

Fig. 4

Weighted averaging , , and for . In each case, we compare the values with , 0.3, and 1.

Predictability of a linear extrapolation,

In this section, we investigate the predictability of a linear extrapolation of . Let us define the linear extrapolation where means the predicted reproduction number based on the previous two consecutive reproduction numbers and the error is defined as

Fig. 5(a) shows and . Up to about three weeks, there are large differences between two reproduction numbers. However, we can observe a good agreement interval from to , see Fig. 5(b). Fig. 5(c) and (d) are corresponding errors of (a) and (b). The errors from to are listed in Table 1 . Except some days, the absolute errors are less than one. We note that we have a spike in Fig. 5(b) and (d). That is due to no population data changes for two days, i.e., 11 July and 12 July as shown in Fig. 1.

Fig. 5

(a) Comparison of the with which are denoted by the symbols () and (), respectively. (b) A magnified view of a rectangular box in (a). (c) Error for from 20 May 2015 to 17 July 2015. (d) A magnified view of a rectangular box in (c) from 18 June 2015 to 17 July 2015.

Table 1

Error which is expressed by from 18 June, 2015 to 17 July, 2015, i.e., .

n	E^n	n	E^n	n	E^n
30	−2.4532	40	−0.1110	50	−0.0002
31	+0.3215	41	+0.0000	51	+0.0001
32	−0.6608	42	−0.0001	52	−0.0001
33	+0.3751	43	−0.1999	53	−0.0002
34	−0.2251	44	+0.2571	54	−1.6766
35	+0.5175	45	−0.4144	55	+3.3534
36	−0.0958	46	+0.6574	56	−1.6766
37	−0.1967	47	−0.1004	57	−0.0000
38	−0.0000	48	−0.1997	58	+0.0009
39	+0.1110	49	+0.0002	59	−0.0000

The SIR model with disease-related death

In this section, we consider the modified SIR model with disease-related death. Let the total population be divided into susceptible , infected , recovered , and dead individuals at time , then the SIR model with fatality rate is given by with initial condition , , , and . Here, the time-dependent parameters , , and denote the contact, recovered, and dead rates from disease, respectively  [22]. Now, by using the proposed method, we evaluate the numerical solution of the modified SIR model (14), (15), (16), (17). We take the same initial conditions: , , , , and . The initial guesses for finding coefficients , , and are 1.0e−6, 1.0e−1, and 1.0e−1, respectively. Using this model, we can predict the ratio of individuals who died by disease to the infected population at a given time. To investigate the ratio, we introduce the fatality percentage rate which is defined as

In Fig. 6 , we show the fatality percentage rate from the real observation data and the numerical approximations. Here, the fatality percentage rate from the numerical approximation is obtained by the numerical solutions , , , and . The numerical solutions , , , and are evaluated by using the best-fitting parameters , , and at time . As shown in Fig. 6, the both results are in good agreement in later times.

Fig. 6

Fatality percentage rate by the observation data (circle marker) and the numerical approximations (solid line).

Conclusion

We have proposed the daily computed weighted averaging basic reproduction number for MERS-CoV outbreak in South Korea, May to July 2015. The linearly extrapolated reproduction number, , has important implications for the future prediction of the trend of the reproduction numbers. We showed that the results computed by our method match very well with the actual data of the MERS-CoV spread in South Korea. Using the method proposed in this paper, one can predict the progress of the spread of infections on the use of real-time epidemic data, so that proper control policies can be performed in an appropriate time for reducing the spread of the infections. As future research, the methodology introduced in this paper can be applied to various diseases with other epidemic spreading models.