Parameter estimation and fractional derivatives of dengue transmission model

1.   Introduction

Dengue fever is one of the most contagious diseases in the tropics and subtropics around the world. The disease is caused by the dengue virus that is transmitted to humans through the bite of Aedes Aegypti female mosquito. After incubation of the virus for 4-10 days, the infected mosquito is capable of transmitting the virus for the rest of its life. The incidence of dengue globally has been widespread in recent decades. About half of the world's population is now at risk of becoming infected with this disease ^[1].

Compared to other diseases and their effects, dengue is a huge burden on human populations, health, and economic systems in most tropical countries in the world. Every year, 100 out of 1000 cases increase, 20,000 of them die ^[2]. There is no specific treatment for dengue, but early detection and access to appropriate medical care reduce mortality. Prevention and control of dengue depend on effective vector control ^[1].

In Indonesia, dengue fever is one of the infectious diseases that until now is still a public health problem. This disease often appears as an extraordinary event because of its rapid and potentially deadly spread. Dengue fever spread throughout Indonesia and attacked all age groups, especially children since it was first discovered in 1968 in Surabaya ^[3]. In 2013, the number of dengue fever patients in Indonesia reported as many as 112,511 cases with the death of 871 people ^[4].

The mathematical modeling has become a powerful tool for understanding the dynamics of the spread of an infectious disease including dengue fever. Some researchers have developed a mathematical model of the spread of dengue fever with a host and vector approach ^[5,6,7]. A mathematical model cannot be interpreted in real cases if the parameter values of the model are unknown. Pandey et al. in ^[8] have applied Bayesian Markov chain Monte Carlo (MCMC) method to estimate the parameters of simple host-vector and SIR (Susceptible, Infectious, Recovery) model using monthly incidence of dengue fever in Thailand. The comparison result show that the SIR model better than the host-vector model to explain the incidence data from Thailand. Gotz et al. ^[9] linked daily data of dengue and rainfall patients in Semarang city, Indonesia to estimate IR (Infectious, Recovery) model with optimal control theory approach. Recently, Agusto and Khan ^[10] have parameterized the dengue model using the classical least-squares method based on the 2017 dengue data in Pakistan.

When a mathematical model was presented in a closed form, then the parameters in the model can be estimated by deterministic optimization methods such as the Newton method or Nelder-Mead method ^[11]. Unfortunately, the Nelder-Mead or Newton method fails to converge into global minimum of a function if the function has many local minima ^[12]. In addition, many mathematical models occur in non-linear ODE (ordinary differential equation) systems, so the exact solution (closed form solution) of the model cannot be obtained. In this case, a heuristic method can be used to estimate parameter values from a non-linear model or a non-linear ODE system. Heuristic method such as Genetic Algorithm (GA) and Particle Swarm Optimization (PSO) method have been previously used to estimate parameters of non-linear ordinary differential equations system model. GA method has been often employed to optimize the parameters of ODE model ^[13,14,15]. So far the PSO method is rarely explored in this field ^[16,17]. PSO method was developed by Kennedy and Eberhart in 1995 ^[18]. This method is inspired by social behavior of organisms in a group e.g. bird swarm, fish swarm and ant colony.

In the present era the fractional order models associated with epidemic models gaining much attention day by day. The development in fractional calculus in terms of fractional derivative and the relevant definitions are proposed, Caputo, Caputo-Fabrizio and Atangana-Baleanu derivative. These derivative have successfully applied to the models of real-life problems and obtained good results ^{[19,20,21,22,23,24,25]}. The Atangana-Baleanu is the recent development in the fractional calculus and getting attention day by day due to its memory effect ^[26]. This derivative is based on nonlocal and non-singular kernel and can handle the problems associated to real life situations. The crossover behavior and the non-local dynamics of the realistic models cannot be utilized efficiently with the integer derivative and to such problems the role of fractional calculus is appreciated. More recently, Qureshi and Atangana ^[27] have compared the integer and fractional derivative dengue models related to dengue fever outbreak in Cape Verde islands around 2009.

In this paper, we propose the parameter estimation techniques with PSO method. This method is applied to estimate the parameters of the host-vector and SIR type dengue fever transmission model based on cumulative data of dengue fever patient in East Java province, Indonesia. Furthermore, we formulate the dengue transmission model of the host-vector and SIR type in Atangana-Baleanu derivative. The numerical procedure for Atangana-Baleanu derivative is then implemented for different values of the fractional order of the model.

2.   Host-vector model for dengue transmission

In this section, we describe a simple host-vector model for dengue transmission. Here we follow the host-vector model from ^[8,19] to describe the dynamics of dengue disease transmission. We also apply the model to explain dengue transmission in East Java province, Indonesia. We included the disease-induced mortality in human population (host). The host-vector model is divided into three human host populations, susceptible ($S$), infectious ($I$) and recovered ($R$), and two mosquitoes (vector) populations, susceptible ($V_S$) and infectious ($V_I$). The system of differential equations describing the host-vector model is written as

where $N = S+I+R$ and $V = V_S + V_I$ denote the total population size of host and vector, respectively. Parameter description of model (2.1) could be seen in Table 1.

The change rate of the total human populations is given by

and the change rate of the total vector populations is as follow

It is noted that in the absence of the disease-induced death ($d = 0$), the total human population, $N \rightarrow \frac{\Lambda}{\mu}$ as $t \rightarrow \infty$. The biologically feasible region of the model (2.1) consisting of

with

and

As the standard analysis of the model, we begin by calculating the basic reproduction number of the model. The basic reproduction number is the expected number of secondary case per primary case in a virgin population ^[28,29]. The model (2.1) has two equilibriums, namely the disease-free equilibrium and the endemic equilibrium. The disease-free equilibrium of model (2.1) is given by $E^0 = (S^0, I^0, R^0, V_S^0, V_I^0) = (\frac{\Lambda}{\mu}, 0, 0, \frac{\Lambda_v}{\mu_v}, 0)$. Using the next generation method ^[30], the basic reproduction number of the model (2.1) is

The disease-free equilibrium of model (2.1) is locally asymptotically stable if $\mathcal{R}_0 < 1$, however the endemic equilibrium locally asymptotically stable if $\mathcal{R}_0 > 1$.

In order to parameterize dengue model (2.1), we assumed that the total human population ($N$) is a constant due to we only fit the data to the model about one year. We also assumed that the total vector populations is the constant. Thus, we choose $\Lambda = \mu N+ d I$ and $\Lambda_v = \mu_v V$. Hence, the mathematical model (2.1) can be written follow

There are six unknown parameters in host-vector model (2.3) i.e. $\beta_h, \mu, \gamma, d, \beta_v$ and $\mu_v$. The six parameters of the model (2.3) were estimated using dengue fever data obtained from East Java province using particle swarm optimization method.

3.   SIR model for dengue transmission

Most of the mathematical models for dengue fever transmission have been developed with the host-vector model approach. However, vector population (mosquito) data in the field is usually difficult to obtain because of the infection dynamics in the vectors are faster than in the hosts (humans). In terms of estimating the model parameters, it is more appropriate to model the spread of dengue using the SIR-type. The dengue transmission model using SIR-type is represented by the following equations

where $S, I$ and $R$ denote the susceptible, infected and recovered human populations respectively. The parameters $\Lambda$ denote the birth rate, $\beta$ denote the rate for dengue transmission and $\mu$ denote natural death rate. By $\gamma$ we denote the recovery rate and $d$ denote dengue death rate. All parameters of the model (3.1) are assumed to be constant and nonnegative. If we assume $N$ denotes the total human population, then we have $N = S+I+R$. The change rate of the total populations is given by

The model (3.1) has the biologically feasible region on $\Omega$ with

The model (3.1) has a disease free equilibrium $E_0 = (S_0, I_0, R_0) = (\frac{\Lambda}{\mu}, 0, 0)$. Using the next generation method as presented in ^[30], the basic reproduction number of the model (3.1) is

Thus, the endemic equilibrium of the model (3.1) is given by $E^* = \left(\frac{(\gamma+\mu+ d)\Lambda}{\mu \beta}, \frac{\Lambda}{\beta}(\mathcal{R}_0-1), \frac{\gamma \Lambda}{\mu\beta}(\mathcal{R}_0-1) \right)$. The disease-free equilibrium of model (3.1) is locally asymptotically stable if $\mathcal{R}_0 < 1$. Meanwhile, the endemic equilibrium is exists and locally asymptotically stable if $\mathcal{R}_0 > 1$.

Next, for the simplicity we assumed the total population ($N$) was constant by using $\Lambda = \mu N+ d I$. Thus, the mathematical model (3.1) can be expressed as follows.

In the next section, we estimate all parameters in the model (3.4). The model has four parameters i.e. $\beta, \mu, \gamma$, and $d$. We apply the genetic algorithm to estimate all parameters in the model.

4.   Parameter estimation using particle swarm optimization method

In this method, solutions are represented by particles. A particle is represented by the position of the particle ($x_i$) and the particle velocity ($v_i$). This method begins with random selection of particle position and particle velocity.

Motion of a particle is influenced by stochastic and deterministic components. Every particle is influenced by global-best position ($g_{best}$) and particle-best position in a particle group ($p_{best}$). In general, the parameter estimation step by using the particle swarm optimization method is as follows:

(a) Evaluate cost function value of every particle. Here we choose mean absolute percentage error (MAPE) as a cost function which is given by

Here $n$ is the number of data, $y_{ij}$ is the real data component, and $y_{ij}^*$ is the component of the solution of the differential equation using the Runge-Kutta method.

(b) Update the particle best ($p_{best}$) and the global best ($g_{best}$). Let the particle swarm optimization method contains $m$ particles. At every iteration, the particle-best position ($p_{(best, i)} (t)$) is replaced by position of current particle $x_i(t+1)$ when the MAPE of $x_i (t+1)$ is lower than the MAPE of $p_{(best, i)}(t)$. Moreover, and the global best position ($g_{best}$) is substituted by position of the current particle $x_i (t+1)$ when the MAPE of $x_i(t+1)$ is lower than the MAPE of $g_{best}(t)$.

(c) Calculate new velocity of every particle by using

Here $r_1$ and $r_2$ are random numbers between zero and one with uniform distribution.

(d) Calculate new position of every particle by using

The steps are reiterated some termination condition is met. Here we choose $c_1 = c_2 = 2$ and $w = 1$.

5.   Parameter estimation of host-vector and SIR model

East Java Province is located in the eastern part of the island of Java, Indonesia. The population is about $N = 41.5$ millions ^[31]. East Java has the widest area among 6 provinces in Java Island and has the second largest population in Indonesia after West Java. The data used for parameter estimation on the SIR-type model is cumulative data of dengue fever case per month from March 2012 - February 2013 presented in Figure 1 ^[3]. The data represent data of dengue epidemic condition in East Java. Cumulative data are generally smoother than actual event data and it is also easier to handle data delays due to holidays and weekends ^[8].

In this section, we estimate parameters of host-vector model (2.3) and SIR model (3.4) by using PSO method. We estimated parameters in the host-vector and SIR such that the mean absolute percentage error (MAPE) is minimum. In the PSO implementation, the number of particles is set to 40 particles. We applied the methods for 20 trials where for every trial the methods were terminated after 250 iterations. The best results of the host-vector model and SIR model are presented in the Table 2 and Table 3 respectively.

From the Tables 2 and 3, minimum of MAPE of the host-vector and SIR model are around 7.77 % and 11.61 % respectively. By using parameter values from the tables, we find the basic reproduction number for the host-vector model and the SIR model are

respectively. The basic reproduction number values from both models indicates dengue endemic condition in East Java province, Indonesia. This prediction result is consistent with the real situation, where dengue infection is not removed yet from the East Java province.

Figures 2 and 3 illustrate both models matching images with real data. From Figure 2 it can be seen that the number of infected humans of the SIR model between estimation results with real data is very close at the end of the observation although there is a small deviation at the beginning. However, the estimated result of the host-vector model for the infected human is very closed for the first 10 months and then there is a small deviation at the end of the observation. Similarly, in Figure 3 indicates that the number of recovered human of both models between estimation results and real data is also quite close until the end of the observation. Hence, we will explore the two dengue models with a fractional model approach in order to see the comparison with real data.

6.   A fractional dengue model

The present section presents the dengue SIR model given in (3.1), will be formulated in Atangana-Baleanu derivative. Therefore, the fractional Atangana-Baleanu representation for the SIR dengue model (3.1) take the shape below:

where, $S$, $I$ and $R$ and their parameters defined in above sections. Further, the initial conditions associated with the model (6.1) are non-negative. Before working on the model (6.1), first, we present the basic definitions related to the fractional Atangana-Baleanu derivative and then, we give the numerical procedure to obtain the graphical results of the model with various fractional values of the fractional order parameter $\alpha$. The Vector host model in the form of Atangana-Baleanu derivative can be expressed as

In the numerical section, we give a numerical procedure for the solution of the SIR model in Atangana-Balenau derivative (6.1) and for the vector host model in fractional derivative (6.2) the readers can easily implement a numerical scheme by following the procedure given for (6.1).

6.1.   Basic concepts of A-B derivative

First, we give some definitions on A-B fractional derivative which will be used later in our study ^[26].

Definition 6.1. Suppose $h\in F^1(a, b)$, $b > a$, $\alpha \in [0, 1]$ then in Caputo sense the newly fractional derivative can be written as follows:

where $B(\alpha) = 1-\alpha+ \frac{\alpha}{\Gamma(\alpha)}$ is a normalized function with $B(0) = B(1) = 1$ and $E_\alpha$ is Mittag-Leffler function

Definition 6.2. Consider $h\in F^1(a, b)$, $b > a$, $\alpha \in [0, 1],$ (not necessary differentiable) then, in Riemann-Liouville (ABR) sense, one can express the newly derivative knows as Atangana-Baleanu fractional derivative as is follows:

Definition 6.3. For the Atangana-Baleanu fractional derivative one can express the fractional integral wit non local kernel as follows:

We obtain the initial function and the classical integral respectively, when the fractional order is zero and 1.

Theorem 6.4. On $[a, b]$, the following inequality holds for $f$ when the function $h$ is continuous on $[a, b]$.

Theorem 6.5. Atangana-Baleanu and Atanagan-Baleanu-R derivatives both fulfill the condition Lipschitz given in the following:

also for ABR derivative we have

Theorem 6.6. For the FDE

one can obtain a unique solution by using the inverse Laplace transform and the convolution result ^[26]:

7.   Numerical procedure

The aim of the present section is to obtain the numerical results of the model (6.1) using the method in ^[32] which is applied efficiently to various kind of real-life problems such as ^[21,22,23]. Now using this technique we apply to our model (6.1) as follows:

For the simplification purpose, we express the model (6.1) in the following form:

System (6.1) can be converted in the following form using the fundamental theorem of integration.

For $t = t_{n+1}$, $n = 0, 1, 2, ...,$ we obtain from Eq (7.2)

Approximating the integral in Eq (7.3), using two-point interpolation polynomial and then simplifying, we finally get the following iterative solution for the model given (6.1). In a similar way for the rest of equations of the system (6.1), we obtained the recursive formula as below

After the successful implementation of the numerical scheme on the fractional model (6.1) as explained above, we find the graphical results of the model (6.1), by considering and assigning values to the fractional parameter $\alpha \in[0, 1],$ and to the model relevant parameters.

The parameters used in numerical simulations are estimated from the real data given in Table 2 for the host-vector model and in Table 3 for the SIR model. The value of recruitment rate of host and vector are $\Lambda = 41,539.567095$ and $\Lambda_v = 787,736.152032$ respectively. The initial values of the host-vector model (2.1) and SIR model (3.1) are given by $(S(0), I(0), R(0), V_S(0), V_I(0)) = (35054487;706;700;556656;33)$ and $(S(0), I(0), R(0)) = (35054487;706;700)$ respectively. We obtained Figure 4 with subetaaphs (a) and (b) by using the real data and its comparison with the fractional order parameter $\alpha$. We make a little change in the parameter value of $\beta_h = 1.361353$ by $\beta_h = 1.461353$ and then checked for different value of $\alpha$ and observed that fractional order provide good fitting to the real data. The graphical results for fractional SIR model (6.1) and vector host model (6.2) are shown respectively in Figures 5–8 and 9–12.

8.   Conclusion

In this study, the particle swarm optimization is implemented to estimate the parameters of the dengue transmission model using host-vector and SIR-type. The parameters of both models are estimated using the cumulative data of dengue fever patient in East Java province, Indonesia. The mean absolute percentage error of the host-vector and SIR model are around 7.77 % and 11.61 % respectively. From the parameter values, the basic reproduction number of both models are greater than one. This prediction result is consistent with the real situation, where dengue infection has not removed yet from the East Java province. Further, we formulated the host-vector and SIR models in Atangana-Baleanu fractional derivative due to the reason of non-singular and non-local kernel which effectively handled the problems of real life phenomenon which is not handled by the ordinary order derivative. Therefore, to use the advantage of the Atanagana-Baleanu fractional order derivative, the results are obtained and discussed with the newly proposed procedure which is applied effectively to many real life problems. The graphical results obtained for both the fractional dengue model (6.1) and (6.2) and considered the value of $\alpha = 1, 0.95, 0.9, 0.85$. The graphical results show that at each instant of time level as well as at arbitrary derivative, one can obtain resealable information which is not possible in the case of integer order derivative. We proven that the vector host model provide better result for the dengue data rather than SIR model.

Acknowledgments

Part of this research is funded by Unversitas Airlangga through Hibah Riset Mandat (No. 342/UN3.14/LT/2019).

Conflict of interest

The authors declare that they have no competing interests.