Estimating the Time Reproduction Number in Kupang City Indonesia, 2016–2020, and Assessing the Effects of Vaccination and Different Wolbachia Strains on Dengue Transmission Dynamics

Abstract

The use of a vaccine and Wolbachia bacterium have been proposed as new strategies against dengue. However, the performance of Wolbachia in reducing dengue incidence may depend on the Wolbachia strains. Therefore, in this paper, the performance of two Wolbachia strains which are WMel and WAu, in combination with the vaccine, has been assessed by using an age-dependent mathematical model. An effective reproduction number has been calculated using the Extended Kalman Filter (EKF) algorithm. The results revealed that the time reproduction number varies overtime with the highest one being around 2.75. Moreover, it has also found that use of the vaccine and Wolbachia possibly leads to dengue elimination. Furthermore, vaccination on one group only reduces dengue incidence in that group but dengue infection in the other group is still high. Furthermore, the performance of the WAu strain is better than the WMel strain in reducing dengue incidence. However, both strains can still be used for dengue elimination strategies depending on the level of loss of Wolbachia infections in both strains.

1. Introduction

The risk of being infected by dengue is higher in tropical and sub-tropical areas where dengue is endemic. It has been estimated that around 105 million infections happen annually with 51 million febrile disease cases [1]. Dengue is caused by four distinct serotypes where being infected by one of the serotypes provides longlife immunity to that serotype but gives a higher chance to get more severe forms of dengue in secondary infection. As there are no effective strategies against all dengue serotypes, the risk for obtaining the severe dengue is possible [2,3,4]. Furthermore, the importation of dengue makes it possible for dengue to spread worldwide easily [5].

Although a number of strategies have been implemented, they were found to be less effective. The current proposed strategies are the use of a vaccine and Wolbachia [6]. The development of a dengue vaccine is underway and its efficacy ranges between 64% and 80% depending on serotypes [7,8,9]. A higher reduction in dengue incidence can be obtained when vaccinating individuals aged 9–45 years [2,4]. Moreover, the risk for secondary infections is higher when vaccinating sero-negative individuals [2,4]. Furthermore, the use of Wolbachia bacterium is a promising strategy [10]. Research found that Wolbachia can reduce the dengue incidence by up to 80% particularly in areas with low to moderate transmission levels [11,12]. A Wolbachia is a bacterium that can reduce the level of virus in the mosquito’s body [13], which reduces the possibility for transmitting the virus when they bite susceptible humans. Furthermore, Wolbachia reduces blood-feeding success in Aedes aegypti [14] and hence minimizes the probability of a successful bite. On the other hand, the use of a dengue vaccine can reduce the dengue transmission in higher transmission regions. The combination of both a vaccine and Wolbachia can significantly reduce the number of dengue infections.

To understand the performance of Wolbachia and the vaccine in reducing dengue infections, various mathematical models have been formulated and studied [15,16,17,18]. A review on dengue modeling for the last 10 years can be found in [19]. Research mostly studied the effects of Wolbachia and the vaccine on dengue transmission dynamics independently and little research has been conducted to assess the effects of the implementation of both strategies simultaneously. Salgado et al. [16] used the dynamic optimization approach to understand the effects of Wolbachia on dengue transmission dynamics and found that the release of Wolbachia can significantly reduce the number of dengue incidence. Ogunlade et al. [18] formulated a mathematical model to study the performance of WAu on dengue transmission dynamics. The results showed that the use of WAu is good at reducing dengue incidence. Shim [20] studied the effects of vaccination on dengue transmission dynamics and found that optimal vaccination rates potentially increase with a higher proportion of seropositive individuals which leads to a higher impact of vaccination. Ndii et al. [4] also studied the effects of vaccination and found vaccinating seropositive individuals provides better reduction in dengue transmission. Not many mathematical models have been formulated to study the effects of both strategies simultaneously. Dorrigati et al. [6] suggested that the use of vaccination and Wolbachia potentially results in a higher reduction in dengue incidence. Several mathematical models have been formulated and studied the use of both strategies simultaneously. Ndii et al. [15] studied the effects of using vaccination and Wolbachia simultaneously and found that the use of the vaccine can potentially reduce the number of dengue infections if the vaccine efficacy is high, otherwise the use of Wolbachia is sufficient [15]. Junsawang [21] performed numerical simulations of a combination of the use of Wolbachia and the vaccine, where the focus was on numerical properties of the model.

The findings have suggested that different Wolbachia strains provide distinct biological characteristics, which potentially leads to their performance in reducing dengue transmission. Several Wolbachia strains that have been used in the field trials are WAu and WMel. These strains give different biological effects on mosquitoes. The characteristics of these two Wolbachia strains are the following [22,23,24,25]. First, WAu has high viral blockage and WMel has only medium. Second, maternal transmission is high in both strains. Third, loss of Wolbachia infections is low in WAu strains and high in WMel. Furthermore, the fitness cost is medium in both strains. WAu has no cytoplasmic incompatibility effects while WMel has cytoplasmic incompatibility. Different biological effects of these strains would provide different performance in reducing dengue incidence. Hence, investigating the effects of different Wolbachia strains in combination with vaccination is of great importance. To the best of our knowledge, no mathematical model has been formulated to study the effects of different Wolbachia strains in combination with vaccination, which is the focus of this work.

In this paper, we formulate a deterministic mathematical model to assess the effects of vaccination and different Wolbachia strains on dengue transmission dynamics. We estimate the effective reproduction number using the Extended Kalman Filter (EKF) algorithm against data of dengue incidence from Kupang-city Indonesia. We employ an optimal control theory to assess the optimal control strategies to result in a higher reduction in dengue incidence, the influential parameters are also determined to gain insights on factors governing the dengue transmission dynamics.

2. Formulation of Mathematical Model

2.1. Modeling Framework for Mosquito Population

For the mosquito population dynamics, we modify the model by Ndii et al. [11] to include the possibility of loss of Wolbachia infections and no effects of cytoplasmic incompatibility. The mosquito population consists of non-Wolbachia and Wolbachia-carrying mosquitoes. The eggs, larvae and pupae are grouped into one compartment called aquatic. The population of adult mosquitoes is divided into male and female mosquitoes. $A_n$, $M_n$, $F_n$, $A_w$, $M_w$, $F_w$ are the aquatic, male and female mosquitoes where the subscripts n and w are to differentiate between non-Wolbachia and Wolbachia-carrying mosquitoes. The model is then governed by the following system of differential equations

$$\begin{array}{cc}\hfill \frac{d{A}_n}{dt}& ={\rho }_n\frac{F_n(M_n+\varphi M_w)}{P}\left(1-\frac{A_n+A_w}{K}\right)-({\gamma }_n+{\mu }_{na})A_n,\hfill \\ \hfill \frac{d{M}_n}{dt}& ={ϵ}_n{\gamma }_n{A}_n+{ϵ}_{nw}(1-\alpha ){\gamma }_w{A}_w+{\gamma }_w{M}_w-{\mu }_n{M}_n,\hfill \\ \hfill \frac{d{F}_n}{dt}& =(1-{ϵ}_n){\gamma }_n{A}_n+(1-{ϵ}_{nw})(1-\alpha ){\tau }_w{A}_w+{\gamma }_w{F}_w-{\mu }_w{M}_w,\hfill \\ \hfill \frac{d{A}_w}{dt}& ={\rho }_w\frac{F_w(M_n+M_w)}{P}\left(1-\frac{A_n+A_w}{K}\right)-({\tau }_w+{\mu }_{wa})A_w,\hfill \\ \hfill \frac{d{M}_w}{dt}& ={ϵ}_w\alpha {\tau }_w{A}_w-{\gamma }_w{M}_w-{\mu }_w{M}_w,\hfill \\ \hfill \frac{d{F}_w}{dt}& =(1-{ϵ}_w)\alpha {\tau }_w{A}_w-{\gamma }_w{F}_w-{\mu }_w{F}_w.\hfill \end{array}$$

(1)

where the $P=M_n+F_n+M_w+F_w$ is the total mosquito population, ${\gamma }_w$ is the rate of loss of Wolbachia infections and $\varphi$ is to denote the effects of Cytoplasmic Incompatibility (CI) where $\varphi =1$ means there is an effect of CI and $\varphi =0$ means no CI effect. The ratio between male and female mosquitoes is approximately 1:1 and hence we set ${ϵ}_n={ϵ}_w={ϵ}_{nw}=1/2$. We then obtain the following Equation:

$$\begin{array}{cc}\hfill \frac{d{A}_n}{dt}& ={\rho }_n\frac{F_n(F_n+\varphi F_w)}{2(F_n+F_w)}\left(1-\frac{(A_n+A_w)}{K}\right)-({\gamma }_n+{\mu }_{na})A_n,\hfill \\ \hfill \frac{d{F}_n}{dt}& =\frac{{\gamma }_n}{2}{A}_n+\frac{(1-\alpha ){\tau }_w}{2}{A}_w+{\gamma }_w{F}_w-{\mu }_n{F}_n,\hfill \\ \hfill \frac{d{A}_w}{dt}& ={\rho }_w\frac{F_w}{2}\left(1-\frac{(A_n+A_w)}{K}\right)-({\tau }_w+{\mu }_{wa})A_w,\hfill \\ \hfill \frac{d{F}_w}{dt}& =\frac{\alpha {\tau }_w}{2}{A}_w-{\gamma }_w{F}_w-{\mu }_w{F}_w.\hfill \end{array}$$

(2)

In the absence of loss of Wolbachia infections and the existence of cytoplasmic incompatibility, the dynamics of mosquito population can be found in Ndii et al. [26].

2.2. Modeling Framework for Host-Vector Involving Wolbachia and Vaccination

In this section, we formulate a deterministic host-vector mathematical model involving Wolbachia and the vaccine, where the vector population is governed by Equation (2). We divided human population into disjoint compartments. The human population comprises susceptible child and adult individuals ($S_c$ and $S_a$, respectively), infected child and adult individuals ($I_c$ and $I_a$, respectively), and recovered child and adult individuals ($R_c$ and $R_a$, respectively). The mosquito population is divided into aquatic mosquitoes ($A_n$ and $A_w$), susceptible ($S_n$ and $S_w$) and infected ($I_n$ and $I_w$) groups, where subscript n and w is to represent non-Wolbachia and Wolbachia-carrying mosquitoes.

Let ${\alpha }_h$ be the progression rate from child to adult and the value of $\alpha =1/T$, where T is the age at which individuals in the child class move to the adult class. The parameters ${ϵ}_1$ and ${ϵ}_2$ are vaccine efficacy on child and adult individuals. The parameters $v_1$ and $v_2$ are the vaccination rate, the $\gamma$ is the recovery rate. $\mathsf{\Lambda }$ is the recruitment rates of humans. The susceptible individuals are infected after being bitten by infectious non-Wolbachia and Wolbachia carrying mosquitoes at a rate of ${\beta }_h^n$ and ${\beta }_h^w$, respectively. Furthermore, the vaccinated individuals move to the recovered class. After a certain period, the recovered individuals move to the susceptible class. The model is given by the following system of differential equation.

$$\begin{array}{cc}\hfill \frac{d{S}_c}{dt}& =\mathsf{\Lambda }-{\alpha }_h{S}_c-\frac{{\beta }_h^n{I}_n}{N_h}{S}_c-\frac{{\beta }_h^w{I}_w}{N_h}{S}_c-{\mu }_h{S}_c-{ϵ}_1{v}_1{S}_c+q_1{R}_c,\hfill \\ \hfill \frac{d{S}_a}{dt}& ={\alpha }_h{S}_c-\frac{{\beta }_h^n{I}_n}{N_h}{S}_a-\frac{{\beta }_h^w{I}_w}{N_h}{S}_a-{ϵ}_2{v}_2{S}_a-{\mu }_h{S}_a+q_2{R}_a,\hfill \\ \hfill \frac{d{I}_c}{dt}& =\frac{{\beta }_h^n{I}_n}{N_h}{S}_c+\frac{{\beta }_h^w{I}_w}{N_h}{S}_c-(\gamma +{\mu }_h)I_c,\hfill \\ \hfill \frac{d{I}_a}{dt}& =\frac{{\beta }_h^n{I}_n}{N_h}{S}_a+\frac{{\beta }_h^w{I}_w}{N_h}{S}_a-(\gamma +{\mu }_h)I_a,\hfill \\ \hfill \frac{d{R}_c}{dt}& =\gamma I_c+{ϵ}_1{v}_1{S}_c-({\mu }_h+q_1)R_c,\hfill \\ \hfill \frac{d{R}_a}{dt}& =\gamma I_a+{ϵ}_2{v}_2{S}_a-({\mu }_h+q_2)R_a,\hfill \\ \hfill \frac{d{A}_n}{dt}& ={\rho }_n\frac{(F_n^2+\varphi F_n{F}_w)}{2(F_n+F_w)}\left(1-\frac{(A_n+A_w)}{K}\right)-({\tau }_n+{\mu }_{na})A_n,\hfill \\ \hfill \frac{d{S}_n}{dt}& =\frac{{\tau }_n{A}_n}{2}+(1-\alpha ){\tau }_w\frac{A_w}{2}+{\gamma }_w{S}_w-\frac{{\beta }_n(I_c+I_a)}{N_h}{S}_n-{\mu }_n{S}_n,\hfill \\ \hfill \frac{d{I}_n}{dt}& =\frac{{\beta }_n(I_c+I_a)}{N_h}{S}_n-{\mu }_n{I}_n,\hfill \\ \hfill \frac{d{A}_w}{dt}& ={\rho }_w\frac{F_w}{2}\left(1-\frac{(A_n+A_w)}{K}\right)-({\tau }_w+{\mu }_{wa})A_w,\hfill \\ \hfill \frac{d{S}_w}{dt}& =\frac{{\tau }_w}{2}\alpha A_w-{\gamma }_w{S}_w-\frac{{\beta }_w(I_c+I_a)}{N_h}{S}_w-{\mu }_w{S}_w,\hfill \\ \hfill \frac{d{I}_w}{dt}& =\frac{{\beta }_w(I_c+I_a)}{N_h}{S}_w-{\mu }_w{I}_w.\hfill \end{array}$$

(3)

with non-negative initial conditions $S_c\left(0\right)\ge 0,S_a\left(0\right)\ge 0,I_c\left(0\right)\ge 0,I_a\left(0\right)\ge 0,R_c\left(0\right)\ge 0,R_a\left(0\right)\ge 0,A_n\ge 0,S_n\left(0\right)\ge 0,I_n\left(0\right)\ge 0,A_w\ge 0,S_w\left(0\right)\ge 0,I_w\left(0\right)\ge 0$ and $N_h=S_c+S_a+I_c+I_a+R_c+R_a$. We can verify that the solutions of the Model (3) with non-negative initial conditions remain non-negative. $F_n=S_n+I_n$ and $F_w=S_w+I_w$.

3. Basic Reproduction Number

The basic reproduction number is generated using the concept of the next generation matrix [27]. We construct the transmission and transition matrices. The transmission matrix, $\mathbb{T}$ and the transition matrix, $\mathsf{\Sigma }$ are

$$\mathbb{T}=\left(\begin{array}{cccc}0& 0& \frac{{\beta }_h^n{S}_c^{*}}{N_h}& \frac{{\beta }_h^w{S}_c^{*}}{N_h}\\ 0& 0& \frac{{\beta }_h^n{S}_a^{*}}{N_h}& \frac{{\beta }_h^w{S}_a^{*}}{N_h}\\ \frac{{\beta }_n{S}_n^{*}}{N_h}& \frac{{\beta }_n{S}_n^{*}}{N_h}& 0& 0\\ \frac{{\beta }_w{S}_w^{*}}{N_h}& \frac{{\beta }_w{S}_w^{*}}{N_h}& 0& 0\end{array}\right)\mathsf{\Sigma }=\left(\begin{array}{cccc}-(\gamma +{\mu }_h)& 0& 0& 0\\ 0& -(\gamma +{\mu }_h)& 0& 0\\ 0& 0& -{\mu }_n& 0\\ 0& 0& 0& -{\mu }_w\end{array}\right),$$

(4)

where $S_c^{*}$ and $S_a^{*}$ are the number of child and adult susceptible individuals at disease free equilibrium.

We then take the inverse of the transition matrix ${\mathsf{\Sigma }}^{-1}$ and obtain

$${\mathsf{\Sigma }}^{-1}=\left(\begin{array}{cccc}-\frac{1}{(\gamma +{\mu }_h)}& 0& 0& 0\\ 0& -\frac{1}{(\gamma +{\mu }_h)}& 0\\ 0& 0& -\frac{1}{{\mu }_n}& 0\\ 0& 0& 0& -\frac{1}{{\mu }_w}\end{array}\right).$$

(5)

The next generation matrix is obtained by $-\mathbb{T}{\mathsf{\Sigma }}^{-1}$ and hence

$$-\mathbb{T}{\mathsf{\Sigma }}^{-1}=\left(\begin{array}{cccc}0& 0& \frac{{\beta }_h^n{S}_c^{*}}{{\mu }_h{N}_h}& \frac{{\beta }_h^w{S}_c^{*}}{{\mu }_w{N}_h}\\ 0& 0& \frac{{\beta }_h^n{S}_a^{*}}{{\mu }_h{N}_h}& \frac{{\beta }_h^w{S}_a^{*}}{{\mu }_w{N}_h}\\ \frac{{\beta }_n{S}_n^{*}}{N_h(\gamma +{\mu }_h)}& \frac{{\beta }_n{S}_n^{*}}{N_h(\gamma +{\mu }_h)}& 0& 0\\ \frac{{\beta }_w{S}_w^{*}}{N_h(\gamma +{\mu }_h)}& \frac{{\beta }_w{S}_w^{*}}{N_h(\gamma +{\mu }_h)}& 0& 0\end{array}\right).$$

(6)

The reproduction number is the spectral radius of the next generation matrix. We obtain the reproduction number, $R_0$, as

$$R_0^{vw}=\sqrt{\frac{(S_c^{*}+S_a^{*})(S_n^{*}{\beta }_h^n{\beta }_n{\mu }_w+S_w^{*}{\beta }_h^w{\beta }_w{\mu }_n)}{{\mu }_n{\mu }_w(\gamma +{\mu }_h)N_h^2}},$$

(7)

It is clear that in the absence of vaccination, $S_c^{*}+S_a^{*}=N_h$ and therefore the basic reproduction number in the absence of vaccination is

$$R_0^w=\sqrt{\frac{(S_n^{*}{\beta }_h^n{\beta }_n{\mu }_w+S_w^{*}{\beta }_h^w{\beta }_w{\mu }_n)}{{\mu }_n{\mu }_w(\gamma +{\mu }_h)N_h}}.$$

In the presence of vaccination, the $(S_c^{*}+S_a^{*})<N_h$, $R_0^{vw}<R_0^w$.

4. Parameter Estimation

4.1. Data

In this study, weekly dengue data from Kupang city since January 2016 to December 2020 has been used to estimate the reproduction number and the transmission rate. The data is given in Figure 1. This is secondary data that has been obtained from the Health Office of Kupang City, Indonesia. The Health Office has collected the data from all Public Health Centers in Kupang city. It can be seen that the highest incidence occurred in 2019 and dengue epidemic has happened annually. Furthermore, the number of infected mosquitoes is set to be three times the number of infected humans.

4.2. Mathematical Model for Parameter Estimation

For this purpose, as the dengue data is for situation in the absence of Wolbachia and vaccination, we estimate using the model in the absence of Wolbachia and vaccination. Hence the Model (3) has been reduced to

$$\begin{array}{cc}\hfill \frac{d{S}_h}{dt}& =\mathsf{\Lambda }-\frac{{\beta }_h^n{I}_n}{N_h}{S}_h-{\mu }_h{S}_h,\hfill \\ \hfill \frac{d{I}_h}{dt}& =\frac{{\beta }_h^n{I}_n}{N_h}{S}_h-(\gamma +{\mu }_h)I_h,\hfill \\ \hfill \frac{d{R}_h}{dt}& =\gamma I_h-{\mu }_h{R}_h,\hfill \\ \hfill \frac{d{A}_n}{dt}& ={\rho }_n\frac{F_n}{2}\left(1-\frac{A_n}{K}\right)-({\tau }_n+{\mu }_{na})A_n,\hfill \\ \hfill \frac{d{S}_n}{dt}& =\frac{{\tau }_n{A}_n}{2}-\frac{{\beta }_n{I}_h}{N_h}{S}_n-{\mu }_n{S}_n,\hfill \\ \hfill \frac{d{I}_n}{dt}& =\frac{{\beta }_n{I}_h}{N_h}{S}_n-{\mu }_n{I}_n.\hfill \end{array}$$

(8)

Using the Next Generation Approach, we had the basic reproduction number of Model (8) is

$$R_0=\sqrt{\frac{{\beta }_h^n{\beta }_n{S}_n^{*}}{N_h{\mu }_n(\gamma +{\mu }_h)}}.$$

Model (8) is discretized using the forward Euler method and, together with Extended Kalman Filter, is used to estimate the transmission rate and the time-varying effective reproduction number for dengue in Kupang city, Indonesia.

4.3. Discrete Time Stochastic Augmented Model

A discrete time stochastic augmented model is constructed. We set the transmission rate ${\beta }_h^n={\beta }_n$ and denoted by $\beta$. We discretize Model (8) using forward Euler and augmenting the transmission rate $\beta$ as a state variable and hence we obtained the following discrete-time stochastic augmented model:

$$\begin{array}{cc}\hfill S_h(k+1)& =S_h\left(k\right)+\mathsf{\Lambda }-\frac{\beta \left(k\right)I_n\left(k\right)}{N_h}{S}_h\left(k\right)\Delta t-{\mu }_h{S}_h\left(k\right)\Delta t+{\omega }_1\left(k\right),\hfill \\ \hfill I_h(k+1)& =I_h\left(k\right)+\frac{\beta \left(k\right)I_n\left(k\right)}{N_h}{S}_h\left(k\right)\Delta t-\gamma I_h\left(k\right)\Delta t-{\mu }_h{I}_h\left(k\right)\Delta t+{\omega }_2\left(k\right),\hfill \\ \hfill R_h(k+1)& =R_h\left(k\right)+\gamma R_h\left(k\right)\Delta t-{\mu }_h{R}_h\left(k\right)\Delta t+{\omega }_3\left(k\right),\hfill \\ \hfill A_n(k+1)& =A_n\left(k\right)+{\rho }_n{\displaystyle \frac{F_n\left(k\right)}{2}}\left(1-{\displaystyle \frac{A_n\left(k\right)}{K}}\right)\Delta t-\left({\tau }_n+{\mu }_{na}\right)A_n\left(k\right)\Delta t+{\omega }_4\left(k\right),\hfill \\ \hfill S_n(k+1)& =S_n\left(k\right)+\frac{{\tau }_n}{2}{A}_n\left(k\right)\Delta t-\beta \left(k\right)\frac{S_n\left(k\right)I_h\left(k\right)}{N_h}\Delta t-{\mu }_n{S}_n\left(k\right)\Delta t+{\omega }_5\left(k\right),\hfill \\ \hfill I_n(k+1)& =I_n\left(k\right)+\beta \left(k\right)\frac{S_n\left(k\right)I_h\left(k\right)}{N_h}\Delta t-{\mu }_n{I}_n\left(k\right)\Delta t+{\omega }_6\left(k\right),\hfill \\ \hfill \beta (k+1)& =\beta \left(k\right)\Delta t+{\omega }_7\left(k\right)\hfill \end{array}$$

(9)

where $\Delta t$ is the time step and $F_n\left(k\right)=S_n\left(k\right)+I_n\left(k\right)$. We add noise

$$\mathit{\omega }\left(k\right)=\left({\omega }_1\left(k\right),{\omega }_2\left(k\right),{\omega }_3\left(k\right),{\omega }_4\left(k\right),{\omega }_5\left(k\right),{\omega }_6\left(k\right),{\omega }_7\left(k\right)\right)$$

to model uncertainty.

4.4. Extended Kalman Filter

The Extended Kalman Filter has been used to estimate the parameter values and the reproduction number [28,29,30]. We estimate the reproduction number by applying the Extended Kalman Filter (EKF) to the discrete-time stochastic augmented compartmental as given in Model (9). Let us define

$$\mathit{x}\left(k\right)={\left(S_h\left(k\right),I_h\left(k\right),R_h\left(k\right),A_n\left(k\right),S_n\left(k\right),I_n\left(k\right)\right)}^T.$$

Model (9) can be written as

$$\mathit{x}(k+1)=\mathit{f}\left(x\right(k\left)\right)+\mathit{\omega }\left(k\right)$$

(10)

where $\mathit{f}$ is the right hand side of Model (9). Denote $\widehat{\mathit{x}}$ as the estimate of $\mathit{x}\left(k\right)$ from the EKF. We set the Jacobian of $\mathit{f}$ at the estimate of $\mathit{x}\left(k\right)$ which is given as

$$J\left(\widehat{\mathit{x}}\left(k\right)\right)=\left(\begin{array}{ccccccc}{J}_{11}\left(\widehat{\mathit{x}}\left(k\right)\right)& 0& 0& 0& 0& J_{16}\left(\widehat{\mathit{x}}\left(k\right)\right)& J_{17}\left(\widehat{\mathit{x}}\left(k\right)\right)\\ J_{21}\left(\widehat{\mathit{x}}\left(k\right)\right)& J_{22}\left(\widehat{\mathit{x}}\left(k\right)\right)& 0& 0& 0& J_{26}\left(\widehat{\mathit{x}}\left(k\right)\right)& J_{27}\left(\widehat{\mathit{x}}\left(k\right)\right)\\ 0& 0& J_{33}\left(\widehat{\mathit{x}}\left(k\right)\right)& 0& 0& 0& 0\\ 0& 0& 0& J_{44}\left(\widehat{\mathit{x}}\left(k\right)\right)& J_{45}\left(\widehat{\mathit{x}}\left(k\right)\right)& J_{46}\left(\widehat{\mathit{x}}\left(k\right)\right)& 0\\ 0& J_{52}\left(\widehat{\mathit{x}}\left(k\right)\right)& 0& J_{54}\left(\widehat{\mathit{x}}\left(k\right)\right)& J_{55}\left(\widehat{\mathit{x}}\left(k\right)\right)& 0& J_{57}\left(\widehat{\mathit{x}}\left(k\right)\right)\\ 0& J_{62}\left(\widehat{\mathit{x}}\left(k\right)\right)& 0& 0& J_{65}\left(\widehat{\mathit{x}}\left(k\right)\right)& J_{66}\left(\widehat{\mathit{x}}\left(k\right)\right)& J_{67}\left(\widehat{\mathit{x}}\left(k\right)\right)\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right)$$

where

$$\begin{array}{cc}\hfill J_{11}\left(\widehat{\mathit{x}}\left(k\right)\right)& =1-\frac{\beta \left(k\right)I_n\left(k\right)}{N_h}\Delta t-{\mu }_h\Delta t,J_{16}\left(\widehat{\mathit{x}}\left(k\right)\right)=-\frac{\beta \left(k\right)S_h\left(k\right)}{N_h}\Delta t,J_{17}\left(\widehat{\mathit{x}}\left(k\right)\right)=-\frac{I_n\left(k\right)S_h\left(k\right)}{N_h},\hfill \\ \hfill J_{21}\left(\widehat{\mathit{x}}\left(k\right)\right)& =\frac{\beta \left(k\right)I_n\left(k\right)}{N_h}\Delta t,J_{22}\left(\widehat{\mathit{x}}\left(k\right)\right)=1-(\gamma +{\mu }_h)\Delta t,J_{26}\left(\widehat{\mathit{x}}\left(k\right)\right)=\frac{\beta \left(k\right)S_h\left(k\right)}{N_h}\Delta t,\hfill \\ \hfill J_{27}\left(\widehat{\mathit{x}}\left(k\right)\right)& =\frac{I_n\left(k\right)S_h\left(k\right)}{N_h}\Delta t,J_{33}\left(\widehat{\mathit{x}}\left(k\right)\right)=1+\gamma \Delta t-{\mu }_h\Delta t,\hfill \\ \hfill J_{44}\left(\widehat{\mathit{x}}\left(k\right)\right)& =1-\frac{{\rho }_n(S_n\left(k\right)+I_n\left(k\right))}{2K}\Delta t-({\tau }_n+{\mu }_{na})\Delta t,J_{45}\left(\widehat{\mathit{x}}\left(k\right)\right)=\frac{{\rho }_n}{2}\Delta t-\frac{{\rho }_n{A}_n\left(k\right)}{2K}\Delta t,\hfill \\ \hfill J_{46}\left(\widehat{\mathit{x}}\left(k\right)\right)& =\frac{{\rho }_n}{2}\Delta t-\frac{{\rho }_n{A}_n\left(k\right)}{2K}\Delta t,J_{52}\left(\widehat{\mathit{x}}\left(k\right)\right)=-\frac{\beta \left(k\right)S_n\left(k\right)}{N_h}\Delta t,\hfill \\ \hfill J_{54}\left(\widehat{\mathit{x}}\left(k\right)\right)& =\frac{{\tau }_n}{2}\Delta t,J_{55}\left(\widehat{\mathit{x}}\left(k\right)\right)=1-\frac{\beta \left(k\right)I_h\left(k\right)}{N_h}\Delta t-{\mu }_n\Delta t,J_{57}\left(\widehat{\mathit{x}}\left(k\right)\right)=-\frac{S_n\left(k\right)I_h\left(k\right)}{N_h}\Delta t,\hfill \\ \hfill J_{62}\left(\widehat{\mathit{x}}\left(k\right)\right)& =\frac{\beta \left(k\right)S_n\left(k\right)}{N_h}\Delta t,J_{65}\left(\widehat{\mathit{x}}\left(k\right)\right)=\frac{\beta \left(k\right)I_h\left(k\right)}{N_h}\Delta t,\hfill \\ \hfill J_{66}\left(\widehat{\mathit{x}}\left(k\right)\right)& =1-{\mu }_n\Delta t,J_{67}\left(\widehat{\mathit{x}}\left(k\right)\right)=\frac{S_n\left(k\right)I_h\left(k\right)}{N_h}\Delta t.\hfill \end{array}$$

In the EKF algorithm, there are two tuning parameters which are the process covariance matrix ${\mathit{Q}}_{\mathit{f}}$ and observation covariance matrix ${\mathit{R}}_{\mathit{f}}$. Note that the tuning parameters are chosen such that the Relative Root Mean Square Error (RRMSE) between the data and the estimated data is sufficiently small. The RRMSE for each variable is defined as

$$RRMSE=\frac{1}{N_w}\sum _{j=1}^{N_w}\frac{\parallel X_j-\widehat{X_j}{\parallel }_2^2}{\parallel X_j{\parallel }_2^2},$$

where $N_w$ is the number of weeks observed and $X_i\in \left(S_h,I_h,I_n\right)$ and ${\widehat{X}}_i\in \left(\widehat{S_h},\widehat{I_h},\widehat{I_n}\right)$.

4.5. Estimation of Reproduction Number

In the estimation of time reproduction number, we have estimated the transmission rate and used its value to calculate the time reproduction number. The other parameters values are obtained from the literature and given in

Table 1. Parameter values and descriptions.

Parameter	Description	Value	Unit
$\mathsf{\Lambda }$	Recruitment rate of human	$\frac{402,286}{65\times 52}$	week ${}^{-1}$
$\mu$	Natural Death Rate	$\frac{1}{65\times 52}$	week ${}^{-1}$
${\alpha }_h$	Progression rate from child to adult	$1/26$	week ${}^{-1}$
${\beta }_h^n$	Transmission rate from Non-W to human	$0.4417$	week ${}^{-1}$
${\beta }_h^w$	Transmission rate from W to human	0.2098	week ${}^{-1}$
$q_1$	Progression rate from recovered child to susceptible child	$1/(6\times 52)$	week ${}^{-1}$
$q_2$	Progression rate from recovered adult to susceptible adult	$1/(12\times 52)$	week ${}^{-1}$
$\gamma$	Recovery rate	7/5	week ${}^{-1}$
${\rho }_N$	Reproduction rate of non-W mosquitoes	8.75	week ${}^{-1}$
${\tau }_N$	Maturation rate of non-W mosquitoes	2	week ${}^{-1}$
${\mu }_{na}$	Death rate of aquatic non-W mosquitoes	1/2	week ${}^{-1}$
$\alpha$	Maternal transmission	0.9	n/a
${\beta }_n$	Transmission rate from human to Non-W mosquitoes	$0.4417$	week ${}^{-1}$
${\rho }_w$	Reproductive rate of W-mosquitoes	$15.75$	week ${}^{-1}$
${\tau }_w$	Maturation rate of W-mosquitoes	2	week ${}^{-1}$
${\gamma }_w$	Loss of Wolbachia infections	$0.28$	week ${}^{-1}$
${\mu }_n$	Death rate of adult non-W mosquitoes	1/2	week ${}^{-1}$
${\mu }_w$	Death rate of adult W mosquitoes	0.45	week ${}^{-1}$

. The plot of the reported data and the simulations is given in Figure 2. It showed that the EKF algorithm estimates the reported data well. Furthermore, the transmission rate and the effective reproduction number are given in Figure 3 where they vary overtime depending on reported data. The highest effective reproduction number has been occurred in 2019 where the $R_t$ is higher than 2.5. It is consistent with the estimate of basic reproduction number by Ndii et al. [31] which estimated using the early growth rate method. It indicates that a single infectious individual may generate around two newly infected individuals. Furthermore, the transmission rate fluctuates between 0 and 2.3 week ${}^{-1}$ with average transmission rate is around 0.4417 week ${}^{-1}$.

5. Sensitivity Analysis

To understand parameters governing the dynamics of dengue transmission in the presence of vaccination and Wolbachia bacterium, a global sensitivity analysis has been performed. We use Latin Hypercube Sampling (LHS) in conjuction with Partial Rank Correlation Coefficient (PRCC) Multivariate analysis to determine the influential parameters. The output of interest is the reproduction number.

Figure 4 showed the results of sensitivity analysis to determine the most influential parameters on the reproduction numbers. It showed that the ${\beta }_h^n$, $S_c$,$S_n$,${\beta }_n$ are the influential parameters and have a positive relationship, which imply the possibility of an increase in the reproduction number when these parameter values increase. Note that the susceptible non-Wolbachia mosquitoes ($S_n$) also govern the dynamics, and hence, we then explore the effects of parameters on the susceptible non-Wolbachia mosquitoes. For this, we measure against an increasing number of susceptible non-Wolbachia mosquitoes which is the solution of

$$\frac{d{C}_{S_n}}{dt}=\frac{{\tau }_n{A}_n}{2}+(1-\alpha ){\tau }_w\frac{A_w}{2}+{\gamma }_w{S}_w.$$

(11)

We found that the population of susceptible mosquitoes is affected by the parameters ${\tau }_n$, ${\mu }_n$ and ${\mu }_{na}$ as given in Figure 5. The loss of Wolbachia infections (${\gamma }_w$) and the maturation rate of Wolbachia (${\tau }_w$) affect in the early period only.

6. Numerical Simulations in the Absence and Presence of Wolbachia

6.1. Numerical Simulations in the Absence of Vaccination

In the numerical simulations, we perform two vaccinations in combination with two different strains of Wolbachia which are WMel and WAlb. We simulate the solutions of the model in the presence of vaccination with different Wolbachia strains. Furthermore, we simulate the model with constant and time-dependent controls. In our simulation, we use the values as given in Table 1 and the transmission rate used is the mean of the estimated values given in Figure 3.

Figure 6 shows that if the loss of infection is high, the number of dengue infection is similar to that in the absence of Wolbachia. This indicates that the use of WMel strain is possible if the loss of infection is not high. This means that a loss of Wolbachia infections affects the performance of Wolbachia in reducing dengue transmission. Hence, the use of WMel may be less effective in comparison to WAu strains.

Figure 7 shows that the performace of the WAu strain is better than WMel. Interestingly, with the same values of loss of Wolbachia infection, WAu still performs better than WMel. This is caused by the effects of CI on WMel. Furthermore, an increase in the loss of Wolbachia infections results in a higher number of dengue incidence. Results suggest that the use of WAu is better in reducing the number of dengue incidence.

6.2. Numerical Simulation in the Presence of Vaccination

In this section, we simulate the presence of vaccination for two different Wolbachia strains. As in the previous section, the performance of WAu is better than WMel. Here, we investigate the effects of vaccination in the presence of Wolbachia. We use an optimal control approach for this.

6.3. Numerical Solution for Dengue Transmission Dynamics with Vaccination with Wolbachia

In this section, numerical solutions of the model in the presence of vaccination and Wolbachia are presented. The results showed that WAu performs better than WMel in reducing dengue incidence and hence in the presence of vaccination, we simulate using the WMel parameter values as used in simulating Figure 7. We consider different vaccination rates.

Figure 8 shows that if the vaccination rate is low, the number of dengue infections is high. To reduce dengue incidence, a higher vaccination rate should be applied. This implies that if a vaccine efficacy and the rate of vaccination is higher in combination with Wolbachia bacterium, dengue elimination is possible.

7. Optimal Control Approach

Numerical simulation showed that the loss of Wolbachia infection governs the dynamics of dengue transmission. When the values of loss Wolbachia infection increases, the number of dengue incidence increases (see Figure 6). Furthermore, the loss of Wolbachia infection is influenced by temperature particularly for WMel strain. To account for this phenomena, in the optimal control analysis, we consider a seasonal loss of Wolbachia infections, and hence, we set the loss of Wolbachia infection rate by the following sinusoidal function

$${\gamma }_w\left(t\right)={\gamma }_{w0}\left(1+\eta sin\left(\omega t\right)\right)$$

where the ${\gamma }_{w0}$ is the average loss of Wolbachia infections, $\eta$ the strength of seasonality. When analysing the effects of WAu strain, the values of ${\gamma }_w=0.28$ and $\varphi =0$. Furthermore, the vaccination acts as a control variable. Let ${ϵ}_1{v}_1$ and ${ϵ}_2{v}_2$ be $u_1$ and $u_2$, respectively.

We aim to minimize the number of human infections with minimal cost. The objective functional is defined as

$$J(u_1,u_2)={\int }_0^{t_f}(A_1(I_c+I_a)+A_2(u_1^2+u_2^2))dt$$

(12)

where $A_1$, $A_2$ are the balancing coefficient and $T_f$ is the final time of interest. We use the quadratic objective as can be found in [32,33,34]. The Hamiltonian functional is defined as

$$\mathcal{H}=(A_1(I_c+I_a)+A_2(u_1^2+u_2^2))+{\lambda }_i\frac{d\mathbf{X}}{dt}$$

(13)

where $\mathbf{X}=\left(S_c,S_a,I_a,I_c,R_a,R_c,A_n,S_n,I_n,A_w,S_w,I_w\right)$ and $i=1,\dots ,12.$

7.1. The Existence and Characterization of Optimal Control

In this section, we present the existence and characterization of the optimal control.

Theorem 1.

There exists an optimal control $\mathit{u}=(u_1,u_2)$ with a corresponding state solutions such that

$$\underset{u_1,u_2\in \mathcal{U}}{min}J(u_1,u_2)=J(u_1^{*},u_2^{*})$$

Proof. 

This theorem is proved based on the results in Flemming and Risher [35] and Lukes [36]. We state the following properties A1–A4:

A1.

The solution of the model with non-negative initial conditions and the associated control function in $\mathcal{U}$ is non-empty.

A2.

The control set $\mathcal{U}$ is convex and closed.

A3.

The dengue model can be expressed as a linear function of $u_1$, $u_2$ with time and state dependent coefficients.

A4.

There exists constant $m_1>0$, $m_2>0$ and $p>1$ such that the integrand in (12) is convex and satisfy

$$H(S_a,I_a,I_c,u_1,u_2)\ge m_1(\sum _{i=1}^2|u_i{{|}^2)}^{\frac{p}{2}}-m_2.$$

The condition A1 and A2 is fulfilled as the state variables and control variables are non-empty and bounded. Condition 3 is met due to the linear dependence of the state system on controls $u_1$ and $u_2$ The algorithm to show this can be seen in [37]. Condition A4 can be verified by writing the following

$$\begin{array}{cc}\hfill (A_1(I_c+I_a)+A_2(u_1^2+u_2^2))& \ge A_2(u_1^2+u_2^2)\hfill \\ \hfill & \ge A_2(u_1^2+u_2^2)-m_2,\hfill \\ \hfill & \ge m_1{(u_1^2+u_2^2)}^{p/2}-m_2\hfill \end{array}$$

Hence, A4 is verified. Therefore, there exists an optimal control $u_1$, $u_2$ which minimizes $J(u_1,u_2)$. □

Theorem 2.

Given the optimal control $(u_1^{*},u_2^{*})$ and the corresponding state trajectories, there exists an adjoit vector function that satisfy

$$\frac{d{\lambda }_i}{dt}=-\frac{\partial H}{\partial x_i}$$

where $i=1,2,\dots ,12$ and $x_i=(S_c,S_a,I_c,I_a,R_c,R_a,A_n,S_n,I_n,A_w,S_w,I_w)$. Furthermore, $u_1^{*}$ and $u_2^{*}$ are characterized by

$$\begin{array}{cc}\hfill u_1^{*}\left(t\right)& =min\left\{max\left\{\frac{S_c({\lambda }_1-{\lambda }_5)}{2A_2},0\right\},u_1^{max}\right\},\hfill \\ \hfill u_2^{*}\left(t\right)& =min\left\{max\left\{\frac{S_a({\lambda }_2-{\lambda }_6)}{2A_2},0\right\},u_2^{max}\right\},\hfill \end{array}$$

(14)

Proof. 

By taking the partial derivative of the Hamiltonian functional with respect to state variables, we then obtain the following adjoint variables.

$$\begin{array}{cc}\hfill \frac{d{\lambda }_1}{dt}& =-{\lambda }_1\left(-{\alpha }_h-\frac{{\beta }_h^n{I}_n}{N_h}-\frac{{\beta }_h^w{I}_w}{N_h}-{\mu }_h-u_1\right)-{\lambda }_2{\alpha }_h-{\lambda }_3\left(\frac{{\beta }_h^n{I}_n}{N_h}+\frac{{\beta }_h^w{I}_w}{N_h}\right)-{\lambda }_5{u}_1,\hfill \\ \hfill \frac{d{\lambda }_2}{dt}& =-{\lambda }_2\left(-\frac{{\beta }_h^n{I}_n}{N_h}-\frac{{\beta }_h^w{I}_w}{N_h}-{\mu }_h-u_2\right)-{\lambda }_4\left(\frac{{\beta }_h^n{I}_n}{N_h}+\frac{{\beta }_h^w{I}_w}{N_h}\right)-{\lambda }_6{u}_2,\hfill \\ \hfill \frac{d{\lambda }_3}{dt}& =-A_1-{\lambda }_3(-\gamma -{\mu }_h)-{\lambda }_5\gamma +\frac{{\beta }_n{S}_n}{N_h}({\lambda }_8-{\lambda }_9)+\frac{{\beta }_w{S}_w}{N_h}({\lambda }_{11}-{\lambda }_{12}),\hfill \\ \hfill \frac{d{\lambda }_4}{dt}& =-A_1-{\lambda }_4(-\gamma -{\mu }_h)-{\lambda }_6\gamma +\frac{{\beta }_n{S}_n}{N_h}({\lambda }_8-{\lambda }_9)+\frac{{\beta }_w{S}_w}{N_h}({\lambda }_{11}-{\lambda }_{12}),\hfill \\ \hfill \frac{d{\lambda }_5}{dt}& =-{\lambda }_1{q}_1-{\lambda }_5(-{\mu }_h-q_1),\hfill \\ \hfill \frac{d{\lambda }_6}{dt}& =-{\lambda }_2{q}_2-{\lambda }_6(-{\mu }_h-q_2),\hfill \\ \hfill \frac{d{\lambda }_7}{dt}& =-{\lambda }_7\left(-{\rho }_n\frac{F_n^2+\varphi F_n{F}_w}{2K(F_n+F_w)}-{\tau }_n-{\mu }_{na}\right)-{\lambda }_8\frac{{\tau }_n}{2}+{\lambda }_{10}\frac{{\rho }_w{F}_w}{2K},\hfill \\ \hfill \frac{d{\lambda }_8}{dt}& =-{\lambda }_8\left(-\frac{{\beta }_n(I_c+I_a)}{N_h}-{\mu }_n\right)-{\lambda }_9\frac{{\beta }_n(I_c+I_a)}{N_h},\hfill \\ \hfill \frac{d{\lambda }_9}{dt}& ={\lambda }_1\frac{{\beta }_h^n{S}_c}{N_h}+{\lambda }_2\frac{{\beta }_h^n{S}_a}{N_h}-{\lambda }_3\frac{{\beta }_h^n{S}_c}{N_h}-{\lambda }_4\frac{{\beta }_h^n{S}_a}{N_h}+{\lambda }_9{\mu }_n,\hfill \\ \hfill \frac{d{\lambda }_{10}}{dt}& ={\lambda }_7{\rho }_n\frac{F_n^2+\varphi F_n{F}_w}{2K(F_n+F_w)}-{\lambda }_8\frac{(1-\alpha ){\tau }_w}{2}-{\lambda }_{10}\left(-\frac{{\rho }_w{F}_w}{2K}-{\tau }_w-{\mu }_{wa}\right)-{\lambda }_{11}\frac{{\tau }_w\alpha }{2},\hfill \\ \hfill \frac{d{\lambda }_{11}}{dt}& =-{\lambda }_8{\gamma }_w\left(t\right)-{\lambda }_{11}\left(-\frac{{\beta }_w(I_c+I_a)}{N_h}-{\mu }_w-{\gamma }_w\left(t\right)\right)-{\lambda }_{12}\left(\frac{{\beta }_w(I_c+I_a)}{N_h}\right),\hfill \\ \hfill \frac{d{\lambda }_{12}}{dt}& =\frac{{\beta }_h^w{S}_c}{N_h}\left({\lambda }_1-{\lambda }_3\right)+\frac{{\beta }_h^w{S}_a}{N_h}\left({\lambda }_2-{\lambda }_4\right)+{\lambda }_{12}{\mu }_w.\hfill \end{array}$$

with the transversality condition ${\lambda }_i\left(t_f\right)=0$, where $i=1,2,3,\dots ,12$. We then take the derivative with respect to $u_1$ and $u_2$ to obtain

$$u_1=\frac{S_c({\lambda }_1-{\lambda }_5)}{2A_2}\mathrm{a}\mathrm{n}\mathrm{d}{u}_2=\frac{S_a({\lambda }_2-{\lambda }_6)}{2A_2}$$

Using the bounds we obtain the characterization of control as given in Equation (14). □

7.2. Numerical Simulations of Optimal Control

In this section, we present numerical simulations of the optimal control. The parameter values used are given in Table 1 and the balancing coefficients are $A_1=A_2=1$. We consider two cases which are vaccination of child individuals only and vaccination of adult individuals only, vaccination both child and adult individuals.

Figure 9 shows the infected individuals with the implementation of vaccination on children only ($u_1\ne 0$ and $u_2=0$). It reveals that the implementation of $u_1$ only can reduce the number of infected children. Dengue infections on adult individuals are still high. Furthermore, the implementation of $u_1$ should be at highest level for the entire period.

Figure 10 shows that the number of infected adult individuals is reduced in the implementation of $u_2$ control only and the number of infected children remains high. The control on $u_2$ is not at the highest level in the early period but then it goes to highest level after around 100 days.

Figure 11 shows that the implementation of vaccination on both child and adult individuals in combination with Wolbachia could possibly eliminate the dengue infections. Furthermore, a highest vaccination rate on children and adult should be implemented to obtain optimal results.

The implementation of vaccination on one group of individuals only cannot eradicate the diseases. It can reduce dengue infections for that group only. Overall, to obtain the optimal results, the vaccination should be implemented in both groups and in combination with the use of Wolbachia.

8. Discussion and Conclusions

In this paper, a mathematical model in the presence of Wolbachia and vaccination has been formulated. The model is studied to examine the performance of both strategies in reducing dengue transmission. We estimate the effective reproduction number using the Extended Kalman Filter algorithm against data of dengue infections from Kupang city, Indonesia. A global sensitivity analysis has been performed to determine the influential parameters. An optimal control approach has been performed to assess the optimal reduction in dengue with minimal cost. The contributions of this paper are the following. The first contribution is the estimation of the time reproduction number in Kupang-city Indonesia. To the best of our knowledge, this is the first work to estimate the time reproduction number in Kupang-city, Indonesia. The second is scientific knowledge regarding the effects of vaccination and the use of WMel and WAu Wolbachia strain in reducing dengue transmission dynamics.

The results show that the Extended Kalman Filter algorithm estimates the reported data well. It reveals that the effective reproduction number fluctuates over time with the highest $R_t$ of of around 2.6 occurs in 2019. This means that a single infectious individual can generate two to three newly infectious individuals. The estimate is consistent with estimates of $R_0$ using early growth rate methods [31]. The results of the effective reproduction number are similar to the previously estimated $R_0$ for different areas such as Cali, Colombia [38], Pakistan [39], and East Java, Indonesia [40]. Furthermore, non-Wolbachia mosquito-related and the transmission parameters are the more influential and govern the dynamics of dengue transmission. This implies that these factors need to be carefully managed to design better actions for reducing dengue incidence. The performance of WAu in reducing dengue is better than WMel. Furthermore, when the rate of loss of Wolbachia infections is high, the implementation of Wolbachia cannot perform well in reducing dengue incidence. The results from the optimal control approach suggest the implementation of vaccination on both groups in the presence of Wolbachia could potentially eliminate the dengue. The implementation of vaccination in certain groups only can reduce the dengue infections in that group only.

Given that the loss of Wolbachia infections affects the performance of Wolbachia in reducing the dengue transmission dynamics, an appropriate use of Wolbachia strains for the dengue elimination strategy is required. For areas with relatively stable temperature over years, the use of either the WMel or WAu strain is appropriate. For areas with strong fluctuations in temperature, the use of the WAu strain is possibly more appropriate. However, both strains can potentially reduce the number of dengue infections. Furthermore, the administration of the vaccine on a single group only, either child or adult, can only reduce the number of infections in that group. However, further research needs to be conducted to understand the effects of temperature and antibody-dependent enhancement on the effects with the implementation of vaccination and different Wolbachia strains. This may add new insights into the use of both strategies. In conclusion, the integrated strategies by using Wolbachia and the vaccine should be implemented in order to possibly reach dengue elimination.