Designing the nationwide emission trading scheme in China

1.   Introduction

Many countries already recognize measles as a highly contagious and dangerous disease for infants and children, which can cause death in the pre-immunization period ^[1,2]. World Health Organization (WHO) has estimated that there were around 535,000 human deaths due to measles in 2000, and most of them occurred in tropical and sub-tropical such as Indonesia, China, Malaysia, Colombia, and Myanmar. Since then, several measles immunization programs have been launched in high-risk countries which have an impact on reducing global measles deaths by 78% in 2012 ^[3]. Dramatic decreases in the number of cases have also been reported in China since 1986 and can reach around 0.5 per million in the 1995–2005 period. It is an impact of the implementation of the two-dose measles vaccination program, which has been carried out during that period ^[4,5].

Meanwhile, a different condition for measles cases was reported in Australia. The outbreaks occur even with a reasonably large reproduction number ($R_0$) of around 15 to 20, whereas efforts to eliminate measles have been launched ^[6]. Consequently, the coverage of the immunization level was increased to 92–95% in order to reduce the number of measles cases ^[7]. Measles outbreaks have also been reported in Italy, with 4347 confirmed cases and 5404 notified cases in 2017 ^[8].

Measles are caused by enveloped RNA virus with the genus Morbillivirus in the family paramyxovirus. Transmission usually occurs through direct contact with sufferers or via air such as mucus contact, coughing, and sneezing ^[9]. Sharing tableware with an infected person may also have a high potential for measles. Other names for this disease are morbilli or rubeola. Measles transmission starts from infected individuals coughing or sneezing who can make the virus spread through the air and can survive up to several hours.

The virus infects the respiratory tract and then spreads throughout the body. The period of measles transmission begins four days before and four days after the rash appears, the peak of transmission at the time of the initial symptoms (prodromal phase), namely in the first 1–3 days of illness ^[11]. Complications that often occur are diarrhea and bronchopneumonia ^[12]. Specific symptoms of this disease are reddish skin rashes that appear within the incubation period (7–18 days) and can last for 5–6 days ^[10].

Symptoms of measles are characterized by fever and are accompanied by one or more symptoms of cough, runny nose, red eyes, or watery eyes and then followed by reddish spots or rash that usually appears starting from behind the ear. Furthermore, symptoms appear on the body in the form of maculopapular shape and found grayish-white patches on the inner cheeks. In children, this disease can cause deadly severe complications if not handled properly. Therefore, early treatment, such as consultation with a doctor or the nearest health care provider, is needed to avoid complications.

The disease can be prevented by immunization using two doses of MMR (measles, mumps, and rubella) vaccine, which is highly effective and safe ^[13]. However, the implementation of immunization strategy has failed to cover rural populations in developing countries ^[14]. One in five children in the world (around 24 million children) did not receive the first dose of measles vaccination in a routine immunization program in 2008 ^[2]. Efforts to overcome the spread of measles globally by providing safe and cost-effective vaccines have been improving rapidly to reduce the mortality of children from measles. The impact was a dramatic decrease in the number of measles cases ^[10]. World Health Assembly endorsed Global Vaccine Action Plan (GVAP) for eliminating measles, congenital rubella syndrome (CRS), and rubella in 2012 ^[15], followed by all WHO regions who agreed to eradicate measles ^[10]. In Indonesia, there were around 23,164 cases of measles and 30,463 cases of rubella from 2010 to 2015. The number of these cases was allegedly lower than the reality, given the number of cases that were not reported ^[16]. In 2016, there were around 3112 measles cases reported in Jakarta ^[17].

Modeling of measles with vaccination provides a better understanding and interpreting collected data and dynamics of measles to produce effective and efficient control strategies. A spatial pattern of misaligned disease data was investigated to estimate and map the risk of measles by utilizing ecological Bayesian regression models in Namibia ^[18]. A discrete-time seasonal model of measles transmission was proposed by Zhenguo Bai and Liu in 2015. They described the seasonal fluctuation of measles cases in China ^[19]. Muhammad Farman et al. have presented analysis and numerical simulations of an SEIR model represented the measles transmission with non-integer time-fractional derivatives in 2017 and considered the measles control ^[20]. An SEIR epidemic model was performed to investigate the impact of exposed individuals in the dynamics of measles transmission ^[21]. The optimal policy for controlling measles was examined through an SIRV model using supplemental immunization activities with separation of vaccinated compartment ^[22]. In 2018, O.J. Peter et al. proposed an SVEIR model of measles transmission by considering natural recovery and recovery due to vaccination ^[23].

The progress in measles transmission modeling has given challenges for researchers in shaping up the health policy decisions. This paper aims to construct two measles transmission models involving vaccination and to investigate the most influential parameters for proposing a control strategy. The proposed model introduces a hospitalized compartment wherein infected individuals entering hospitals are considered isolated and cannot transmit the infection to the susceptible group. First, we presented the weekly data and the measles progress in the data section. We then extracted measles incidence data to be the number of incidences per total population. Second, two mathematical models were constructed, and then afterward, the equilibrium points, as well as their stability conditions, were analyzed. Further, there exist oscillatory solutions if the parameters are satisfying a certain condition. Third, the unobserved parameters and initial conditions were estimated from the data along their feasible regions. Fourth, we did a sensitivity analysis to analyze the effect of changing parameter values and the most influential parameter towards the dynamics. We then conducted several simulations of control strategies to minimize the number of infected by taking different vaccinated proportions in the last section of this paper.

2.   Epidemiological data

2.1.   Measles progress in Jakarta

Following the policy of the Ministry of Health in achieving the elimination target for measles in 2020. the Jakarta Provincial Health Office has targeted for no endemic area of measles for more than 12 months with the implementation of adequate measles surveillance. The indicator is to conduct measles surveillance based on individual cases or CBMS (Case-Based Measles Surveillance), which has been carried out gradually since 2014. In this surveillance, for every case found with symptoms of fever, rash on the body, accompanied by one or more symptoms of cough, runny nose, and red-eye, blood specimens (serum) from the individual are examined at a national reference laboratory to confirm the diagnosis of Measles. Another indicator of the surveillance is discarded ratio (misclassification of measles) obtained from laboratory tests ^[24].

Based on a leaflet from Jakarta Provincial Health Office, hospitals are instructed to identify people groups who are vulnerable to measles by considering the mapping of cases and immunization status at the rural level. It is implemented as a global target for measles elimination and Rubella's control in 2020. Besides, several attempts were made to follow up on each reported case in a health service facility. They are examining, treating, and verifying alerts by carrying out epidemiological investigations in less than 24 hours. Then they are taking and sending specimens for laboratory examinations, reporting each case of measles to Jakarta Health Office surveillance website by submitting individual data (form C1), and monitoring measles cases in the local area or various sources of infection, including cases reported at the local hospital.

Achievement of the measles target is carried out by laboratory sampling in 2017 amounted to 77%, 2018 amounted to 75%, and until July 2019 amounted to 70% based on the evaluation results of measles program holder in Jakarta. Even though the proportion of measles until July 2019 had only reached 70%, the achievement of discarded measles in Jakarta had exceeded 2 per 100,000 population, equal to 4.17 per 100,000 population ^[25].

2.2.   Measles data

In this study, we validated our model with actual data of measles in Jakarta. The data recorded here is the number of hospitalized individuals. We chose Jakarta as our observed area due to the most available of completed data collection compared to other cities in Indonesia. Jakarta Health Office has collected measles data from 44 sub-district health centers, 283 village health centers, and 152 hospitals in Jakarta, which has verified and passed the validation process. Here, the data used is a weekly basis of incidence data started from the first to the last week in 2017, see Figure 1. The left figure is a plot of time series measles data in all districts of Jakarta. The right one is the weekly total incidence of measles in Jakarta, produced by adding up the data in the left data plot.

3.   Methods

We developed a standard SIR model (Kermack–McKendrick epidemic model^[26]) for measles transmission involving vaccination and the addition of a hospitalized compartment to analyze the dynamics of measles transmission. Hereafter a stochastic model was constructed from a theory of stochastic differential equations containing a Wiener process.

3.1.   Deterministic model

We formulated a deterministic model (SIHR) to describe the measles transmission dynamic. The total population is assumed to be constant in this model. To illustrate the model equations, we divide the total population into four compartments. Let $\bar{S}, \bar{I}, \bar{H},$ and $\bar{R}$ are the number of susceptible, infected, hospitalized, and recovered individuals including individuals who were successfully vaccinated, respectively.

The susceptible compartment $\bar{S}$ increased by birth at a rate of $\Lambda$. The individuals probably got their infection by contact with infected individuals at rate $\beta$. Vaccination process is given against the susceptible population $\bar{S}$ at rate $rp$, where $r$ and $p$ are a probability of success in preventing measles and the proportion of vaccinated individuals, respectively. Consequently, the susceptible individuals might be moved to the infected compartment $\bar{I}$ at rate $\beta(1-p)+\beta(1-r)p = \beta(1-rp)$. Then, the infected individuals tend to go to the hospital at a rate of $\delta$ for getting treatment. The infected and hospitalized individuals recovered naturally at a rate of $\gamma_1$ and due to treatment at a rate of $\gamma_2$, respectively. Logically, it should be $\gamma_2 > \gamma_1$. All compartments decreased by a natural death at a rate of $\mu$.

The transmission dynamics of measles can be expressed by a system of differential equations as follows.

where $\dfrac{d\bar{N}(t)}{dt} = \Lambda-\mu N(t)$. We assumed that the population is constant then $\Lambda = \mu N, N\in \mathbb{N}$. We then normalized the variables by substituting $S(t) = \dfrac{\bar{S}(t)}{N}, I(t) = \dfrac{\bar{I}(t)}{N}, H(t) = \dfrac{\bar{H}(t)}{N}, R(t) = \dfrac{\bar{R}(t)}{N}$. Consequently, it reduced the system (3.1) into three-dimensional space in (3.2) by ignoring the recovered compartment due to the fact that recovered and successfully vaccinated individuals did not contribute to the transmission.

Since, the parameters and state variables are non negative ie $t \geq 0$, we investigated model (3.2) in the following feasible region:

All of the parameters used here are positive and described in Table 1. The non-negative octant $\Omega$ is positively invariant of model (3.2). This follows immediately evaluating normal vector at each plane $S = 0, I = 0$, and $H = 0$, respectively. The model has two equilibria; we called it as a disease-free equilibrium (DFE) $E_0 = (S^0, I^0, H^0) = \left(\dfrac{\mu}{rp+\mu}, 0, 0\right)$ and a endemic equilibrium $E_1$ that can be expressed as follows.

where $R_0 = \dfrac{\beta (1-rp) \mu}{(\mu+rp)(\delta+\gamma_1+\mu)}$ is a basic reproduction number for the model. It can be defined as the spectral radius of matrix $\rho(FV^{-1})$ which is

Further, we give the proof of stability of the equilibria both locally and globally, in Theorem 3.1.1 and 3.1.2. The existence of oscillatory solutions is also presented in Proposition 3.1.1.. Although it is trivial, we stated it explicitly to support simulations and to provide easy understanding.

Theorem 3.1.1 Local stability of the equilibria.

1. If $R_0 < 1$, then the disease-free equilibrium $E_0$ is locally asymptotically stable in $\Omega$.

2. If $R_0 > 1$, then the endemic equilibrium $E_1$ is locally asymptotically stable in $\Omega$.

Proof.

1. The Jacobian matrix of model (3.2) at $E_0$ is

The eigen values of the matrix are $-(\mu+\gamma_2), -(rp+\mu),$ and $(\delta+\gamma_1+\mu)(R_0-1)$. If $R_0 < 1$ then all eigen values are real negative, implies the disease-free equilibrium $E_0$ is locally asymptotically stable in $\Omega$.

2. The Jacobian matrix of model (3.2) at $E_0$ is

The characteristic polynomial of $J(E_1)$ is given by $f(\lambda) = \lambda^3+A\lambda^2+B\lambda+C$ where $A = (pr+\mu)R_0+(\mu+\gamma_2), B = (R_0-1)(pr+\mu)(\delta+\gamma_1+\mu)+R_0(pr+\mu)(\mu+\gamma_2),$ and $C = (\gamma_2+\mu)(R_0-1)(pr+\mu)(\delta+\gamma_1+\mu)$. If $R_0 > 1$ then coefficients of the polynomial are positive and $AB > C$. It indicates that Routh-Hurwitz conditions for the polynomial are satisfied. Therefore the endemic equilibrium $E_1$ is locally asymptotically stable in $\Omega.$

Theorem 3.1.2 Global stability of the equilibria.

1. If $R_0 \leq 1$, then the disease-free equilibrium $E_0$ is globally asymptotically stable in $\Omega$.

2. If $R_0 > 1$, then the endemic equilibrium $E_1$ is globally asymptotically stable in the sub-region $\Omega_2 \subseteq \Omega$ that satisfy $H\geq \max \{H^*, \dfrac{\delta I}{\hat{\gamma}}\}$ or $H\leq \min \{H^*, \dfrac{\delta I}{\hat{\gamma}}\}$.

Proof.

1. We first considered the Lyapunov function $V(S, I, H) = \dfrac{S^0}{(\delta+\gamma_1+\mu)}I+\dfrac{1}{\delta}H$ and got

which is less than zero in $\Omega$ if $R_0\leq 1$. A dot represents differentiation with respect to time. Furthermore, $\dot{V} = 0$ if $I = I^0$ and $H = H^0$. Hence, the largest compact invariant set contained in the set $\Psi = \{(S, I, H)\in \Omega, \dot{V}(S, I, H) = 0\}$ when $R_0 < 1$ is the singleton $E_0$. Therefore, the DFE is globally asymptotically stable in $\Omega$ based on the LaSalle invariance principle theorem in ^[29].

2. For simplifying this proof, we assumed that $\hat{\beta} = \beta(1-rp), \hat{\alpha} = rp+\mu, \hat{\delta} = \delta+\gamma_1+\mu$ and $\hat{\gamma} = \gamma_2+\mu$. Consequently, the system (3.2) can rewrite as follows.

where $R_0 = \dfrac{\hat{\beta}\mu}{\hat{\alpha}\hat{\delta}}, S^* = \dfrac{\mu}{R_0 \hat{\alpha}}, I^* = \dfrac{(R_0-1)\mu}{R_0 \hat{\delta}},$ and $H^* = \dfrac{\delta(R_0-1)\mu}{\hat{\gamma}R_0\hat{\delta}}$. We then considered a Lyapunov function

Clearly $U:\Omega \rightarrow \mathbb{R}$ and continuously differentiable. In addition, $U(S^*, I^*, H^*) = 0$ while $U(S, I, H) > 0, \forall (S, I, H) \neq (S^*, I^*, H^*)$. The derivative with respect to the time can be derived as follows.

If $R_0 > 1$ and $H\geq \max \{H^*, \dfrac{\delta I}{\hat{\gamma}}\}$ or $H\leq \min \{H^*, \dfrac{\delta I}{\hat{\gamma}}\}$ are satisfied then $\dot{U}\leq 0$. Therefore, $E_1$ is a globally asymptotically stable within $\Omega_2.$

Proposition 3.1.1 There exists oscillatory solutions of system (3.2) to the endemic equilibrium if $R_0 > 1$ and $I^* > \dfrac{\mu(rp+\mu)R_0}{4(\delta+\gamma_1+\mu)^2}$.

Proof. The characteristic polynomial of the Jacobian matrix at $E_1$ can be used as a sign of the existence of oscillatory solutions around the equilibrium of $E_1$. Remember that we can check it through its discriminant. The eigenvalue of the matrix has non-zero imaginary parts if the discriminant $\varDelta = \dfrac{A^3C}{27}-\dfrac{A^2B^2}{108}-\dfrac{ABC}{6}+\dfrac{B^3}{27}+\dfrac{C^2}{4}$ is bigger than zero ^[30,31]. Therefore, we have $\varDelta = -\dfrac{1}{108}(pr+\mu)Q_1^2 Q_2$ where $Q_1 = R_0(pr+\mu)(\delta+\gamma_1-\gamma_2)+\mu(2\gamma_2-pr-\delta-\gamma_1)-pr(\delta+\gamma_1)+\gamma_2^2 \geq 0$ and $Q_2 = R_0^2(pr+\mu)+(1-R_0)4(\mu+\delta+\gamma_1)$.

If $R_0 > 1$ and $I^* > \dfrac{\mu(pr+\mu)R_0}{4(\delta+\gamma_1+\mu)^2}$ then $Q_2 < 0$, implies $\varDelta > 0.$

3.2.   Stochastic model

In this part, we constructed a system of It$\hat{o}$ SDEs (stochastic differential equations) with continuous-valued random variables by utilizing the probability transition rates in Table 2 and assuming $\triangle t$ (time step) approach zero ^[32,33]. A diffusion approximation is a designation of the process resulting from the SDE model by Kurtz in 1978 ^[34]. Let $y(t)$ denote the number of compartments in the population at generation $t$, $x(t)$ denote the proportion of compartments in a population, and $x(t) = y(t)/N$. The random vectors $x(t) = \left(\begin{matrix} s(t) & i(t) & h(t) \end{matrix}\right)^T$ and $y(t) = \left(\begin{matrix} \bar{s}(t) & \bar{i}(t) & \bar{h}(t) \end{matrix}\right)^T$ are vectors of continuous random variables with state space $s(t), i(t), h(t) \in [0, 1]$ and $\bar{s}(t), \bar{i}(t), \bar{h}(t) \in [0, N], t\in [0, \infty)$.

Based on the transitions and rates in Table 2, the corresponding SDE model can be written as the form

where $G$ is a diffusion matrix and $\bar{W}(t) = (\bar{W}_1(t), \bar{W}_2(t), \cdots, \bar{W}_9(t))^T$ is a $9 \times 1$ matrix (vector) of independent Wiener processes. To order $\triangle t$,

where $\bar{f}$ is drift vector and

Next, we computed the expectation and the covariance matrix based on Table 2. For a sufficiently small $\triangle t$, the expectation

The covariance matrix can be written as

where $m_{11} = \mu+\beta(1-rp)s(t)i(t)+(rp+\mu)s(t),$ $m_{22} = \beta(1-rp)s(t)i(t)+(\delta+\gamma_1+\mu)i(t),$ $m_{33} = \delta i(t) +(\gamma_2+\mu)h(t)$, $m_{12} = -\beta(1-rp)s(t)i(t)$ and $m_{23} = -\delta i(t)$. The diffusion matrix $G(x(t))$ is given by

such that $GG^T = V$.

3.3.   Parameter estimation

We defined an optimization problem for parameter estimation in model (3.2). For simplification of the optimization problem, the model can be rewritten

where $X(t) = \left(\begin{matrix}S(t) & I(t) & H(t)\end{matrix}\right)^T$ and $\vartheta = \{\beta, \delta, I(0), H(0)\}$ is sets of estimated parameters and initial conditions. In this case, we have two unobserved parameters, $\beta$ and $\delta$, and two unknown initial conditions, $I(0)$ and $H(0)$. $p$ is vaccination related parameter, which is not included into consideration for parameter estimation. This parameter is then used to compare the impact of vaccination against the measles epidemic. Hence, the optimization problem can be written as follows.

$\Psi$ is squared residual norm (resnorm). The main idea of parameter estimation is how to find the best combination of $\vartheta$ such that minimizing the resnorm. In this case, we used nonlinear least squares method to look for the best combination of $\vartheta$.

3.4.   Sensitivity analysis

In this part, we discuss a local sensitivity analysis. It gives us information on which parameters are the most influential to the dynamics of the system, with the initial value of the given parameter. We investigate the sensitivity of parameters that are meaningful in the biological interpretation of measles transmission. The idea of this sensitivity analysis is to determine the effect of changing parameter values on the solution of the system $\dot{X} = f(t, X, \theta)$ in Eq (3.8) where $X(t, \theta)\in \mathbb{R}^3$ depend on time $t$ and $\theta$ is the set of parameters, $\theta = \left\lbrace \beta, p, \delta, \gamma_1, \gamma_2 \right\rbrace \in \mathbb{R}^5$. We first define a sensitivity function of $X(t, \theta)$ as follows.

where $i = 1, ...5, k = 1, 2, 3,$ and $\zeta$ is time dependent $t$ ^[35]. By differentiating $\zeta$ with respect to time $t$, then the first order solution sensitivity is a solution of differential equation,

with $\zeta(t = 0) = 0$. $J(t, X, \theta)$ and $Q(t, X, \theta)$ are 3-by-3 Jacobian matrix and 5-by-3 matrix, respectively. The entries of matrix $Q$ are derivatives of the right side of Eq (3.10) with respect to the parameters.

3.5.   Control strategy

In the past decade, some countries in the European Union (EU) are burdened with cases of high measles due to low immunization coverage ^[36]. Several factors potentially reduced coverage of the measles vaccine, such as low family income ^[37], parent education, lack of access to health care networks, migration ^{[37,38,39,40]}, and beliefs about the risk of vaccination ^[41]. Therefore, high vaccine coverage is needed to prevent outbreaks in the general population as well as in the poorly pooled population ^[42]. The effectiveness of immunization also depends on reducing disparities associated with socioeconomic levels and access to health care ^[43,44].

Measles vaccine is available and circulating in society. Jakarta provincial office also supports measles vaccination programs to minimize the cases as much as possible by implementing several government programs that have been launched. There are three stages of vaccine implementation in Jakarta. The stages of vaccination are given to susceptible individuals aged 9–11 months, 18–24 months, and 7–15 years^[45]. However, the efficiency of this intervention program has not been fully assessed. Hence it is necessary to analyze it to get the best strategy that is efficient and adjusted for the budget. Here, we discussed the optimal way to control measles transmission via model (3.2) by a term of vaccination $p$. The parameter values obtained from the estimation parameter are used in this simulation by varying the value of $p$ to observe the impact of vaccination.

4.   Results and discussion

We confirmed the constructed model based on the above analysis results and the actual data in this part. The numerical method used in this simulation is Runge-Kutta method. First, the oscillatory solution was performed by taking parameter values from Table 1, $\beta = 2, p = 0,$ and $\delta = 0.5$. The other values are taken according to either references or assumptions, see Table 1. Second, the solution oscillates and tends to the equilibrium point, see Figure 2. It is a numerical simulation of Proposition 3.1.1 about the existence of oscillatory solutions around $E_1$.Third, we estimated parameters in the model satisfying the minimum error between data and hospitalized solution and used the parameters for stochastic simulation. There are two unobserved parameters in model (3.2), they are $\beta$ and $\delta$. For the initial condition, we assumed $N = 9.5\cdot 10^6$ as a total population in Jakarta, $S(0) = 1-(I(0)+H(0)), R(0) = 0$. The unknown initial conditions, $H(0)$ and $I(0)$ are estimated together with the unobserved parameters. Fourth, we presented a study of sensitivity analysis and then performed control a strategy which is simulated with several levels of immunization coverage.

The result of parameter estimation shows that $\vartheta$ converges to $\beta = 0.67, \delta = 0.3583, I(0) = 2.0157\cdot10^{-5}, H(0) = 2 \cdot10^{-6}$ with a resnorm of $1.0486\cdot10^{-10}$, see Figure 4(a). The residual plot and fitting the distribution of the residual are shown in Figure 4(b) and Figure 4(c), respectively. The results show that the absolute errors are mostly below $2\cdot 10^{-6}$, indicated in the yellow area, which can be interpreted that the average estimated error is under two people per million population. Further, the residual follows a normal distribution with mean equal to zero and $1.4242\cdot 10^{-6}$ variance by using a one-sample $t-$test at significance level 99 %. Afterward, we used the estimated parameter for stochastic simulation with $\triangle t = 10^{-4}$. The results are shown in Figure 5. The simulation results show that the data are within the 95 % confidence interval of the stochastic model and the average of the stochastic solutions across a thousand runs is relatively close to the deterministic solution.

In sensitivity analysis, we are interested in observing the parameter sensitivity under endemic conditions $(R_0 > 1)$ such that the sensitivity analysis was conducted locally in neighborhood of $\theta^\star$ where,

The results of this sensitivity analysis are shown in Figure 6 by taking fixed biological parameters $\mu = 1/(52 \cdot 70)$ and $r = 0.9$. We then compared which parameters are dominantly contributing to the change of the solutions. The parameters discussed here for sensitivity analysis are $\beta, p, \delta, \gamma_1,$ and $\gamma_2$, respectively. Recalling the Eq (3.10), if $\zeta$ is large enough, it means that a small change in parameter $\theta_i$ makes a big contribution to $X_k$. Therefore, significant parameters have large values of $\zeta$. A positive or negative value of $\zeta$ indicates the effect which is directly or inversely proportional, respectively.

In Figure 6(a)–(c), the most sensitive parameter to the susceptible population is $p$. The second and third are $\delta$ and $\gamma_1$. They have a positive effect on the susceptible population even though changes to these parameters do not greatly affect the susceptible population. It indicates that when the parameters increase, then the susceptible population also increases. Unlike the $\beta$ parameter, this parameter is the fourth sensitive parameter after $\delta$ and $\gamma_1$. But the effect of changing the parameter $\beta$ and $p$ have an opposite effect on the population. It means that by increasing infection rates or the proportion of vaccinated individuals, it decreases the susceptible population. Meanwhile, the parameter $\gamma_2$, the effect of the changes is not significant on the susceptible population.

For the sensitivity of the infected population with respect to the parameters, $\delta$ and $\gamma_1$ are the most sensitive parameter and have a opposite effect on the infected population. The parameter $p$ is more sensitive than the parameter $\beta$. But the parameter $p$ has a opposite impact on the infected population, see Figure 6(d)–(f). It can be interpreted that when the rate of treatment to the hospital increases, the infected population decreases. The government (health department) can issue policies such as the health insurance program to increase the rate of treatment at the hospitals to reduce the number of infected populations. Natural healing is more sensitive than vaccination parameters. It indicates that when an outbreak is occurring, providing treatment accesses easier for infected individuals is better than vaccinating.

The most sensitive parameter in the hospitalized population is the recovery rate caused by treatments ($\gamma_2$). The second one is the hospitalized rate ($\delta$). The recovery rates, both natural healing, and healing caused by treatments have an opposite effect on the hospitalized population. It can be seen that the increase in the recovery rate will further reduce the number of hospitalized population. The infection rate ($\beta$) has a positive effect on the hospitalized population even though the sensitivity level is lower than $p$, see Figure 6(g)–(i). Vaccination parameter ($p$) has an opposite impact, which means that giving vaccines with a large proportion will reduce the number of hospitalized population.

Jakarta has about 2.5 million children under 15 years old or around 26% of the total population. They are the most dominant measles vaccination target. Therefore, the variation of vaccination coverage is less than 30% in this simulation. Variation of the level of vaccine coverage is represented by the parameter $p$ and the results are presented in Figure 7. The far-left graphics represent the dynamics infected with respect to the vaccine coverage 0.2% to 1% and 2% to 6% in the right image. Both provide information that increasing vaccine coverage, the number of cases is getting smaller. Nevertheless, this level of vaccine coverage must be adjusted and synchronized with the intended budget to overcome this disease.

5.   Conclusion

We have constructed deterministic and stochastic models for a measles transmission with vaccination. Two equilibria are obtained with their stability, both locally and globally. Oscillatory solutions are also found around the endemic equilibrium when certain conditions are satisfied, are shown in Figures 2 and 3. In the case of parameter estimation, the residuals are mostly less than two people per one million population and follow a normal distribution with a mean of 0 at the significance level of 99%. The stochastic model is shown to capture the randomness of the transmission, represented by independent Wiener process variables and diffusion matrices. It is caused by the influence of other factors in which not all of them can be accommodated into the model. However, the simulation results show that the actual data are within the 95% confidence interval of the stochastic model, which means that the model can reach the random transmission of this disease contained in the data (see Figure 5). The average of the stochastic solutions across a thousand runs is relatively close to the deterministic solution. Moreover, the sensitivity analysis results show that the most influential parameters in the infected population are hospitalized rate and natural recovery rate, then followed by the vaccinated proportion. It can be a valuable information for decision-makers that improving access to health care networks is one of the dominant factors in reducing the number of measles cases. Vaccine coverage is also very influential in reducing the spread of the virus as a preventive measure. We hence proposed a control strategy which is providing treatment accesses easier is better than vaccinating when an outbreak occurs.

There are several limitations to this study. The first problem is the unique solution for parameter estimation. It brings us to the global optimization problem in high dimensional space, which is still a challenge for researchers in the optimization field. The assumption that the total population is almost constant over time and some uniformity assumptions on parameter values taken are also limitations of this study. In the case of sensitivity analysis, it is a local sensitivity in the neighborhood of $\theta^*$, which is time-dependent dynamics, and its sensitivities have adhered to the initial parameters.

Acknowledgment

MF and HF gratefully acknowledge the financial support provided by Indonesia Ministry of Research and Technology through PMDSU Program [grant No. 1511/E4.4/2015]. Parts of this work of ES and NN are funded by ITB Research Grant 2019 and ITB Research, Community Service, and Innovation Program (P3MI-ITB) 2019, respectively [grant No. 1000M/I1.C01/PL/2019].

Author's contributions

MF performed the mathematical analysis, stochastic model, and constructed the manuscript. DS did the mathematical analysis and sensitivity analysis. S collected and interpreted data. HF did the sensitivity analysis. NN reviewed the manuscript, control strategy, and explained the results. ES reviewed the manuscript and analyzed the mathematical method.

Conflict of interest

The authors declare there is no conflicts of interest.