Analysis of yellow fever prevention strategy from the perspective of mathematical model and cost-effectiveness analysis

1.   Introduction

Yellow fever is a hemorrhagic fever disease that is transmitted by the Aedes or Haemagogus aegypti infected mosquito, which contains the yellow fever virus. The virus that causes this disease belongs to the genus Flavivirus, a large group of RNA viruses. It is one of the most dangerous infectious diseases. The mortality rate of this disease ranges from 20 to 50%, but in severe cases, it can exceed 50%. There is no specific treatment for this disease, but it can be reduced through vaccination ^[1]. Although the vaccine for yellow fever is lifelong, a booster is needed to guarantee the efficacy of this vaccine to protect humans from infection.

Various interventions can be performed to suppress the spread of yellow fever. Based on ^[2], fumigation is the most recommended form of intervention to control mosquito populations. Even so, this intervention needs to be carefully considered in its application in the field because it can trigger mosquito resistance to fumigation materials ^[3]. In addition to fumigation, vaccination is highly recommended, especially for tourists entering areas where yellow fever is endemic. Although it is claimed that yellow fever vaccination lasts forever, a vaccine booster must maintain the vaccine's ability to protect the individual from infection at maximum levels ^[4]. Although the yellow fever vaccine is claimed to reach the maximum protection at 4th week after the injection, the implementation of the yellow fever vaccine has some side effects for high-risk groups. These groups include infants aged 6–8 months, individuals more than 60 years old, pregnant, and breastfeeding women ^[3]. Despite the long history of safety, there are some cases of worsening effect of the use of yellow fever vaccine ^[1]. The author of ^[1] expresses more concern about yellow fever vaccine-associated neurotropic disease (YEL-AND) because it is fatal, despite the low number of reported cases.

Mathematical model have been used by many authors to model the dynamic of disease, whether it in human ^[5,6,7,8,9], animals ^[10], or plants ^[11,12].Compared with other vector-borne diseases such as malaria or dengue, as far as we know, there are a few of mathematical models discussing yellow fever. The author of ^[13] introduced a yellow fever transmission model involving two types of hosts, namely, humans and primates. They found that their basic reproduction number consists of the addition of reproduction from human-mosquito interaction and primate (monkey)-mosquito interaction. Additionally, ^[14] introduced a yellow fever model that considers the vaccination of the human population and aquatic phase of the mosquito. Recently, ^[15] modified the model by ^[14] by adding a vertical transmission in the mosquito population to their model. Furthermore, they also involved treated bed nets, larvicides, and adulticides as alternative means of intervention. They showed that their model exhibits backward bifurcation phenomena that make the eradication program for yellow fever more difficult. A mathematical model was used to analyze the impact of vaccination for yellow fever in Angola, as analyzed by ^[16].

Based on the above description, we introduce a new mathematical model for yellow fever transmission, considering vaccination, fumigation, and treatment. We also include the possible side effects of the yellow fever vaccine, which may lead to YEL-AND phenomena. Our model includes the aforementioned strategies as time-dependent variables. We aim to minimize the number of infected individuals keeping the cost of interaction as low as possible. Hence, the most cost-effective possible strategies are then analyzed using the cost-effectiveness method ^[17]. To our knowledge, this work is the first yellow fever model to consider the side effects of the yellow fever vaccine with a time-dependent intervention.

The remainder of this paper is organized as follows. In Section 2, we develop our model using time-dependent interventions. Treating the interventions as a constant parameter, we analyze our special case model regarding its equilibrium points and local stability criteria in Section 3. The characterization of the optimal control problem and its numerical experiments are presented in Section 4. The conclusions of our study are presented in Section 5.

2.   The construction of the yellow fever optimal control model

In this study, we introduced a new mathematical model to describe the spread of yellow fever. Our model includes vaccination $(u_1(t))$, hospitalization $(u_2(t))$, and fumigation $(u_3(t))$ as time-dependent control variables. We aimed to understand the impact of the worsening effect of vaccination intervention on the dynamics of yellow fever and determine the best possible combination of cost-effective strategies in the yellow fever eradication policy effort.

First, we divide the total human population (denoted by $N(t)$) based on their health and vaccination status into eight compartments, namely, susceptible humans $(S(t))$, vaccinated humans $(V(t))$, unvaccinated exposed humans $(E_1(T))$, vaccinated exposed humans $(E_2(t))$, infected humans without vaccination worsening effect $(I(t))$, infected humans with worsening effect of vaccine $(T_1(t))$, infected humans who get worsening effect from the vaccine and were hospitalized $(T_2(t))$, and recovered humans $(R(t))$. Additionally, we divide the mosquito population into susceptible and infected mosquitoes, denoted by $U(t)$ and $W(t)$, respectively. Owing to the short life expectancy of mosquitoes, we exclude the recovered compartments of mosquitoes. This is because once the mosquitoes are infected by the yellow fever virus, they will stay infected forever. Hence, the total human and mosquito population at time $t$ is expressed as follows:

We use the transmission diagram in Figure 2 to construct our yellow fever model, and the description of the parameters is given in Table 1. For written simplification, we use $S, V, E_1, E_2, I, T_1, T_2, R, U$, and $W$ instead of $S(t), V(t), E_1(t), E_2(t), I(t), T_1(t), T_2(t), R(t), U(t)$, and $W(t)$.

In our model, we assume that all recruitment rates enter the susceptible population through a constant rate $\Lambda_h$. Owing to a vaccination rate of $u_1$, the susceptible population will be transferred into compartment $V$. We assume that after a period of $\omega_h^{-1}$, the effect of the vaccine disappears. Hence, the vaccinated individual was sent back to $S$. In our model, we assume that vaccination cannot provide perfect protection against yellow fever infection. Hence, there is an efficacy level of the vaccine, which is denoted by $1-\xi$. Thus, when a susceptible individual is infected with a success rate of $\beta_h$, vaccinated individuals will have a success probability of $\xi \beta_h$.

Susceptible individuals who are already infected will be classified into the exposed compartment $E_1$, which has no chance of experiencing any worsening vaccine effect. However, vaccinated individuals who are already infected and experienced any worsening vaccine effect will be classified into the exposed compartment $E_2$. In our model, we assume that a yellow fever vaccine can have side effects on several types of individuals, including severely immunodeficient individuals, elderly individuals, pregnant women, newborns under six months of age, and individuals who already have thymic diseases, including thymoma and thymectomy ^[18]. These side effects can result in severe symptoms of infection. Because of this worsening effect, after the incubation period of $\gamma_2^{-1}$, an exposed individual $E_2$ will progress into infected stage $I$ with a proportion of $1-p$, where he/she has no worsening effect of the vaccine, and a proportion of $p$ will progress to $T_1$ when he/she has the worsening effect of the vaccine. In contrast, the exposed individual $E_1$ progresses to $I$ after the incubation period of $\gamma_1^{-1}$. Note that $\gamma_1^{-1}\leq \gamma_2^{-1}$.

Infected individuals in $T_1$ are assumed to have more severe symptoms than $I$ individuals. Hence, individuals in $T_1$ must receive more intensive care in hospital. Therefore, we assume that there is a hospitalization rate of $u_2$ that will bring an individual in $T_1$ to a hospitalized individual in $T_2$. The infected individuals in $I$ are assumed to recover from yellow fever after $\delta_1^{-1}$, whereas $T_1$ and $T_2$ are assumed to recover after $\delta_1^{-1}$ and $\delta_2^{-1}$, respectively. We assume that yellow fever may cause death in infected individuals $I, T_1$, and $T_2$ with a constant rate $\alpha_1, \alpha_2,$ and $\alpha_3$, respectively.

However, the mosquito population will increase because of a constant natural birth rate of $\Lambda_v$ and decrease because of the natural death rate $\mu_v$ and fumigation $u_3$. All newborns in the mosquito population were assumed to be susceptible. Additionally, we assume that susceptible mosquitoes may get infected by biting infected individuals who do not go to the hospital, $I$ and $T_1$, with a constant infection rate $\beta_v$.

Our aim is to understand how the spread of yellow fever may be reduced under three types of intervention: vaccination $(u_1)$, hospitalization $(u_2)$, and fumigation $(u_3)$. In so many realities of disease eradication in the field, the limited costs of implementing interventions are always an unavoidable limitation. Therefore, instead of treating the above mentioned interventions as constant parameters, it would be better to make the intervention parameter both time-dependent and adaptive to the number of infected people and environmental conditions at time $t$. One approach is to find the optimal value of $u_1$ for each time $t$. This can be achieved by constructing a model to be discussed using an optimal control approach. Therefore, we discuss the time-dependent value of $u_i$, namely, $u_1(t), u_2(t)$, and $u_3(t)$.

Based on the above descriptions and formulations, the mathematical model of yellow fever under vaccination, hospitalization, and fumigation, as well as the worsening effect of vaccination, is given by the following system of ordinary differential equations:

Note that the yellow fever model in system (2.1) is completed with non-negative initial conditions. Additionally, the parameters in system (2.1) are assumed to be non-negative.

This study aimed to minimize the number of infected individuals as much as possible while keeping the costs of vaccination, hospitalization, and fumigation as low as possible. Because we only want to minimize the number of infected individuals, we set the weight cost for $S, V, R$, and $U$, denoted by $\omega_1, \omega_2, \omega_8,$ and $\omega_9$, respectively, to zero. Hence, we describe our aim in the following cost function:

where $\varphi_i$ for $i = 1, 2, 3$ denotes the weight cost for the control variables $u_1, u_2,$, and $u_3$, respectively. The interval $[0, T]$ represents the duration of the simulation time for the intervention that must be implemented. The total cost for interventions, such as vaccination, hospitalization, and fumigation, is presented as follows:

The quadratic cost function was chosen because of a convex function that ensures the existence of an optimal solution. On the other hand, the quadratic cost function is biologically interpreted as a nonlinear cost for intervention. For example, the cost of fumigation in a small area may only require a cost for buying the fumigation material and hiring an operator. However, when the area that needs to be fumigated becomes much wider, the cost includes not only the purchase of fumigation material and the hiring of operator, but also coordination costs, counseling before implementation, etc. Another quadratic cost function for another type of epidemiological optimal control model can be found in ^[28,29]. On the other hand, all costs associated with the high number of infected humans and mosquitoes, which are not related to vaccination, hospitalization, or fumigation intervention were assumed to be proportional to the number of infected individuals, which is presented as follows:

Our model's aim is to minimize the number of infected individuals in human and mosquito populations while keeping the costs of vaccination, hospitalization, and fumigation as low as possible. Therefore, we seek to optimize the value of $\hat{u_i}$ for $i = 1, 2, 3$ such that

where

is the admissible control set. Note that $u_i^{\min}$ and $u_i^{\max}$ for $i = 1, 2, 3$ represent the lower and upper bounds of each control variable, respectively.

3.   Dynamical analysis of autonomous model

In this section, we analyze the autonomous version of system (2.1) regarding its existence and local stability criteria of equilibrium points, the form of the basic reproduction number, and bifurcation analysis. To accomplish this, let us assume that all control variables remain constant over time. Hence, using $u_1(t) = u_1, u_2(t) = u_2,$ and $u_3(t) = u_3$, the autonomous version of system (2.1) is given by:

Hence, the dynamics of the total human population is given by

and the dynamics of the total mosquito population is given by

3.1.   Preliminary analysis

Before we begin analyzing the model in system (3.1) regarding the existence and stability of its equilibrium points, it is essential to ensure that the model in system (3.1) is well-defined biologically and mathematically. That is, the solution in system (3.1) must always be non-negative for all $t > 0$. Properties related to these are written in the theorems below, where Theorem 1 deals with the positivity of the solution of each variable in system (3.1), and Theorem 2 discusses the upper bound of the total human and mosquito populations.

Theorem 1. Given that the initial values of the system (2.1) are $S(0)$, $V(0)$, $E_1(0)$, $E_2(0)$, $I(0)$, $T_1(0)$, $T_2(0)$, $R(0)$, $U(0)$, and $W(0) \geqslant 0$, the solution of $S(t)$, $V(t)$, $E_1(t)$, $E_2(t)$, $I(t)$, $T_1(t)$, $T_2(t)$, $R(t)$, $U(t)$, and $W(t)$ will always be non-negative for all time $t \geq 0$.

Proof. Please see Appendix A for the proof.

Theorem 2. The total human $(N)$ and mosquito $(M)$ population in system (2.1) is eventually bounded as $t \rightarrow \infty$.

Proof. The proof can be found in Appendix B.

Based on Theorems 1 and 2, we can conclude that each variable in system (3.1) is also bounded.

3.2.   Yellow fever-free equilibrium point and the basic reproduction number

The purpose of this section is to determine the possible final state of the dynamics of system (3.1), which at times tends to infinity and is known as the equilibrium point. The equilibrium of system (3.1) is obtained by setting the right-hand side of system (3.1) to zero and solving this non-linear set of equations with respect to each variable on human and mosquito populations. The first equilibrium point of system (3.1) is the trivial point, known as the yellow fever-free equilibrium point $(\mathscr{E}_1)$. This equilibrium point represents a situation in which all infected populations in humans and mosquitoes, as well as a recovered compartment in humans, do not exist. In this study, $\mathscr{E}_1$ is given as follows:

From the expressions of $S^{\dagger}$ and $V^{\dagger}$, $S^{\dagger}+V^{\dagger} = \frac{\Lambda_h}{\mu_h}$, which means that the existence of vaccination does not change the size of the total population at $\mathscr{E}_1$. Instead, it depends only on the number of newborns and the natural death rate. On the other hand, note that $\frac{S^{\dagger}}{V^{\dagger}} = \frac{\mu_h + \omega_h}{u_1}.$ Our aim was to avoid the rapid spread of yellow fever by expecting many vaccinated individuals. In this study, this aim can be read as minimizing this ratio. To reduce the ratio of $S^{\dagger}$ to $V^{\dagger}$, we can either increase the number of vaccinated individuals (increasing the value of $u_1$) or improve vaccine quality by extending vaccine duration (reducing the value of $\omega_h$). Additionally, note that the total number of mosquitoes at $\mathscr{E}_1$ is inversely proportional to the fumigation rate $u_3$. Hence, it is important to provide a large (but controlled) rate of fumigation to reduce the number of aedes mosquitoes. However, in many reports, uncontrolled fumigation may trigger a mosquito resistance to some fumigants ^[30].

Having the yellow fever-free equilibrium point in hand, as shown in Eq (3.2), we can calculate the basic reproduction number $(\mathscr{R}_0)$ of system (3.1). In the context of our problem, the basic reproduction number presents an expected number of secondary yellow fever cases caused by one primary case during one infection period in a completely susceptible population. Several studies have demonstrated that basic reproduction numbers play an important role in determining the existence or extinction of diseases in many epidemiological models ^[28,29,31]. These studies demonstrated that the disease will likely disappear if $\mathscr{R}_0 < 1$, and it must appear if $\mathscr{R}_0 > 1$. However, it is not always the case that $\mathscr{R}_0 < 1$ will always end in a situation where the disease disappears in the field. One of the causes is the occurrence of a backward bifurcation phenomenon in the model. When backward bifurcation occurs, $\mathscr{R}_0$ is no longer becomes a sufficient condition for the disappearance of a disease because the condition $\mathscr{R}_0 < 1$ causes a bistability phenomenon to occur. More examples of backward bifurcation phenomena arising from epidemiological models are available in ^[32].

Using the next-generation matrix approach ^[33], we have the basic reproduction number of the system (3.1) as follows:

with

The derivation of $\mathscr{R}_0$ is available in Appendix C. Using Theorem 2 in ^[34], the local stability criterion of $\mathscr{E}_1$ is given by the following theorem:

Theorem 3. The yellow fev-free equilibrium point of system (3.1) is locally asymptotically stable if $\mathscr{R}_0 < 1$, and unstable if $\mathscr{R}_0 > 1$.

Proof. Please see Appendix D for the proof.

3.3.   The yellow fever endemic equilibrium

In this section, we analyze the existence of the non-trivial equilibrium point of system (3.1). The other equilibrium point of system (3.1) besides the yellow fever-free equilibrium point is the yellow fever endemic equilibrium point, which is given by the following expression:

where

where $W^{\ddagger}$ is taken from the positive roots of the following quadratic polynomial of $W$:

with

Based on the expression of $F(W)$, we have the following theorem regarding the existence criteria of $\mathscr{E}_2$ in $\mathbb{R}_+^6$:

Theorem 4. The yellow fever endemic equilibrium $(\mathscr{E}_2)$ of system (3.1) always exists if $\mathscr{R}_0 > 1$, and no yellow fever endemic equilibrium otherwise.

Proof. Please see Appendix E for the proof.

3.4.   Forward bifurcation

In this section, we analyze the type of bifurcation of the yellow fever model in (3.1). In many classical models, forward bifurcation in epidemiological models indicates that $\mathscr{R}_0$ is a necessary condition to ensure disease extinction (^{[5,31,35,36,37]}). However, in several epidemiological models ^{[28,38,39,40]}, another type of bifurcation, namely, backward bifurcation, may appear. Under these circumstances, the condition of the basic reproduction number being less than unity is no longer became sufficient to guarantee the extinction of the disease. This is because when backward bifurcation phenomena arise, bistability phenomena may appear where stable endemic equilibrium exists together with the stable disease-free equilibrium.

To analyze the type of bifurcation of our model in (3.1), we use the result of Theorem 11 in ^[41], which is commonly known as the Castillo-Song bifurcation theorem. To apply this theorem, let us re-symbolize our variable from $S, V, E_1, E_2, I, T_1, T_2,$ and $R$ into $x_i$ for $i = 1, 2 \dots 10$, respectively. Hence, we now have the system (3.1) that reads as follows:

Next, we choose $\beta_h$ as the bifurcation parameter. Therefore, solving $\mathscr{R}_0^2 = 1$ with respect to $\beta_h$ (an argument similar to the proof of Theorem 4), we obtain the following:

where $Z_0 = u_{{1}} \left(\gamma_{{1}}+\mu_{{h}} \right) \left(\mu_{{h}}+ \left(\delta_{{1}}-\delta_{{2}}-u_{{2}}+\alpha_{{1}}-\alpha_{{2}} \right) p+\alpha_{{2}}+\delta_{{2}}+u_{{2}} \right)$$\gamma_{{2}}\beta_{{v}}\Lambda_{{v}}\Lambda_{{h}}\xi.$ Linearization of system (3.6) at $\mathscr{E}_2$ with $\beta_h = \beta_h^*$ gives

where

and $A = ((-{\mu_{{h}}}^{2}+(\xi\, u_{{1}}-\alpha_{{2}}-\delta_{{2}}-\omega_{{h}}-u_{{2}})\mu_{{h}}+$$((\delta_{{2}}+u_{{2}}-\alpha_{{1}}+\alpha_{{2}}-\delta_{{1}}) p-\alpha_{{2}}-\delta_{{2}}-u_{{2}})$$u_{{1}}\xi-\omega_{{h}} (\alpha_{{2}}+\delta_{{2}}+u_{{2}}))$$\gamma_{{1}}+ (-\mu_{{h}}+(\delta_{{2}}+u_{{2}}-\alpha_{{1}}+\alpha_{{2}}-\delta_{{1}})$$p-\delta_{{2}}-u_{{2}}-\alpha_{{2}}) u_{{1}}\mu_{{h}}\xi).$ From the direct calculation, it can be shown that $\mathcal{J}$ has a simple zero eigenvalue, while the other nine eigenvalues are negative. Hence, we can use the center manifold theory to analyze the bifurcation of our yellow fever model.

Following the application of the Castillo-Song bifurcation theorem, we must calculate the right and left eigenvectors related to the zero eigenvalue. The right eigenvector of $\mathcal{J}$ at eigenvalue zero is given as follows:

where

with

However, the left eigenvector of $\mathcal{J}$ for eigenvalue zero is given as follows:

where

Using the formula as follows:

we have the following:

where $M_1 = \left(\xi\, \gamma_{{1}}\gamma_{{2}}u_{{1}}+\xi\, \gamma_{{2}}\mu_{{h}}u_{{1}}+\gamma_{{1}}\gamma_{{2}}\mu_{{h}}+\gamma_{{1}}\gamma_{{2}}\omega_{{h}}+\gamma_{{1}}{\mu_{{h}}}^{2}+\gamma_{{1}}\mu_{{h}}\omega_{{h}}\right)$$\left(\alpha_{{3}}+\delta_{{3}}+\mu_{{h}}\right)\Lambda_{{v}}\Lambda_{{h}}\beta_{{v}}$ while $M_0 > 0$ cannot be shown in this study because of its long expression. Because we always have $\mathcal{A} < 0$ and $\mathcal{B} > 0$, we have the following results.

Theorem 5. The yellow fever model in system (3.1) always exhibits a forward bifurcation at $\mathscr{R}_0 = 1$.

Based on Theorem 5, we understand that $\mathscr{E}_2$ is locally asymptotically stable when $\mathscr{R}_0 > 1$ but is close to one. Additionally, based on our proposed model, we can conclude that a basic reproduction number less than unity will always ensure the extinction of yellow fever from the population.

3.5.   Sensitivity analysis, bifurcation and autonomous simulation

In this section, we perform some numerical experiments based on our previous results on Theorems 3, 4, and 5. First, we analyze the elasticity of each parameter in $\mathscr{R}_0$ to determine the most influential parameter for changing the size of $\mathscr{R}_0$. Based on this result, we continue our simulation by showing the level set of $\mathscr{R}_0$ with respect to each control variable to understand the sensitivity of each parameter. The forward bifurcation diagram shows the impact of $\mathscr{R}_0$ on the endemic size of the infected mosquito population. Finally, numerical experiments were conducted to demonstrate the dynamic behavior of the system (3.1) under certain scenarios.

To conduct the elasticity analysis of $\mathscr{R}_0$ with respect to any parameter $\rho$, we use the following equation ^[42]:

To perform these experiments, we use the following parameter values:

which gives $\mathscr{R}_0 = 1.504$. Using the Eq (3.9), the elasticity of $\mathscr{R}_0$ with respect to each parameter in the system (3.1) with the value of the above parameter is shown in Table 2. Note that because $\alpha_3$ and $\delta_3$ do not appear in the expression of $\mathscr{R}_0$, $\epsilon_{\mathscr{R}_0}^{\alpha_3}$ and $\epsilon_{\mathscr{R}_0}^{\delta_3}$ are both zero. The positive sign of $\epsilon_{\mathscr{R}_0}^{\rho}$ means that $\mathscr{R}_0$ increases when $\rho$ increases. However, when $\epsilon_{\mathscr{R}_0}^{\rho}$ is negative, $\mathscr{R}_0$ decreases as $\rho$ increases. Hence, all control variables can reduce the value of $\mathscr{R}_0$ whenever these control parameters increase. The value of $\epsilon_{\mathscr{R}_0}^{\rho}$ represents the percentage of increase or decrease of $\mathscr{R}_0$ for a $1\%$ increase or decrease in $\rho$. For example, because $\epsilon_{\mathscr{R}_0}^{u_3} = -0.2308$, a $1\%$ increase in $u_3$ will reduce $\mathscr{R}_0$ by $0.2308\%$. Additionally, note that $|\epsilon_{\mathscr{R}_0}^{u_3}| = \max\left\{|\epsilon_{\mathscr{R}_0}^{u_1}|, |\epsilon_{\mathscr{R}_0}^{u_2}|, |\epsilon_{\mathscr{R}_0}^{u_3}|\right\}$. Hence, we conclude that fumigation is more essential in reducing $\mathscr{R}_0$ compared to vaccination and hospitalization. These results are in accordance with the recommendations of the WHO regarding fumigation as one of the most recommended efforts to overcome yellow fever ^[3].

In addition to understanding the impact of all parameters in Eq (2.1) to $\mathscr{R}_0$, we must understand the impact of these parameters on the dynamics of all variables in Eq (2.1). Hence, we conducted a sensitivity analysis using three methods: non-normalization, half normalization, and full normalization using the SimBiology toolbox for MATLAB ^[43,44].

Suppose that the yellow fever model in Eq (2.1) has $10$ compartments, $x_{i}$ in $i = 1, 2, ..., 10$ represents $S, V, E_1, E_2, I, T_1, T_2, R, U$, and $W$, respectively. Furthermore, let $k_{j}$ for $j = 1, 2, ..., 20$ represent all the parameters in system (2.1). We assume that the model-balanced equations are expressed as a system of differential equations as follows:

where $x \in \mathbb{R}^{n}$, $k \in \mathbb{R}^{m}$, and $f_i(x, k)$ is the right-hand side of the system (2.1). The non-normalization sensitivities are given as follows:

where $\mathcal{S}^{x_{i}}_{k_{j}}$ is the sensitivity coefficient of each $x_{i}$ with respect to each parameter $k_{j}$. The half-normalization sensitivities are also given by the following equation:

Furthermore, the full-normalization sensitivities are defined as follows:

The results of the sensitivity analysis without normalization are shown in Figure 3. The results shown in Figure 3 show that the infection parameters $\beta_h$ and $\beta_v$ are very dominant in influencing changes in the populations of $S$ and $E_1$ in humans, and the populations of $U$ and $W$ on mosquitoes. This agrees with the results of the analysis in Table 2 regarding the effect of $\beta_h$ and $\beta_v$, which are also dominant on $\mathscr{R}_0$. Furthermore, it was also observed that the fumigation intervention ($u_3$) was more dominant in suppressing the number of infected mosquitoes, whereas vaccination ($u_2$) was dominant in determining the dynamics of healthy and vaccinated human populations. The simulation results in Figure 3 also confirm that hospitalization is the least effective among all control parameters.

We perform a half-normalization of our model in Figure 4. Note that $\beta_h$ and $\beta_v$ are still dominant in determining the dynamics of the (3.1) model compared to other parameters. Interestingly, fumigation is not only dominant in regulating population changes in mosquitoes but also in infected human populations. Clearly, the presence of fumigation has a significant impact on the population of $E_1$. This is consistent with the fact that the expected effect of fumigation is a reduction in the number of new infections in the human population. Figure 5 shows the sensitivity analysis using the full normalization approach. The results indicate the same result as using the non-and half-normalization approaches, where fumigation is more dominant compared to other control parameters.

Next, we show the bifurcation diagram of system (3.1) using the same parameter values as in previous numerical experiments. We use our results in Theorem 5, which indicates that our proposed model in system (3.1) always exhibits a forward bifurcation at $\mathscr{R}_0 = 1$. The bifurcation diagram is shown in Figure 6. Note that the yellow fever-free equilibrium is always locally asymptotically stable when $\mathscr{R}_0 < 1$. At $\mathscr{R}_0 = 1$, the yellow fever-free equilibrium becomes unstable, while the yellow fever endemic equilibrium starts to arise and is locally asymptotically stable. Note that a larger value of $\mathscr{R}_0$ will increase the size of infected mosquitoes in an endemic equilibrium.

To illustrate the dynamics of infected compartments with respect to the values of $\mathcal{R}_0$, we apply the Runge-Kutta method to solve system (3.1). In the first simulation we set the values of all parameter similar to the previous simulation, except $u_1 = 0.03, u_2 = 0.01, u_3 = 0.03$. This setting gives us $\mathcal{R}_0 < 1$. In Figure 7, we can see that all trajectories tend to the yellow fever free equilibrium point asymptotically, which in this case is given by

To simulate $\mathcal{R}_0 > 1$, we choose the same parameter value as before, except $\beta_h = 2\times10^{-5}, \beta_v = 2\times10^{-5}, u_1 = 0, u_2 = 0, u_3 = 0$. In Figure 8, we can see that all solutions tend to the endemic equilibrium point, which is given by

4.   Characterization and simulation of the optimal control problem

4.1.   Characterization of the optimal control problem

Optimal control has been used in several epidemiological models to analyze how the best form of intervention should be chosen to reduce the rate of spread of the disease, where the intervention is treated as a time-dependent variable ^{[45,46,47,48]}. In these studies, the optimal condition for the control variables concludes that the optimal dynamics of the controls depend on the population involved in the model. These various studies concluded that time-dependent interventions provide more optimal results in suppressing the spread of disease but at a more optimal cost.

Before we characterize our optimal control problem, we recall all the control variables introduced in Section 2 as follows:

● $u_1(t)$: Vaccination rate

● $u_2(t)$: Hospitalization rate

● $u_3(t)$: Fumigation rate.

The objective function is expressed using Eq (2.2). Our aim is to minimize the number of infected individuals with optimal interventions, which are read as follows:

where $\mathcal{U}$ is the admissible control set.

We use the Pontryagin maximum principle (PMP) ^[49] to determine the necessary condition for the existence of a triplet of our optimal controls. We define the Hamiltonian of the problem as follows:

where $\lambda_i$ for $i = 1, 2, \dots 10$ are the adjoint variables for the state variables $S, V, E_1, E_2, I, T_1, T_2, R, U$, and $W$, respectively. Hence, the necessary condition of our optimal control problem is given by the following theorem:

Theorem 6. For the optimal control variables $\hat{u}_1, \hat{u}_2,$, and $\hat{u}_3$, and the optimal solutions $\hat{S}, \hat{V}, \hat{E_1}, \hat{E_2}, \hat{I}, \hat{T_1},$ $\hat{T_2}, \hat{R}, \hat{U},$ and $\hat{W}$ of system (2.1), the adjoint variables $\lambda_i$ for $i = 1, 2, \dots 10$ that satisfy the following:

with transversality conditions $\lambda_i(T) = 0$ for $i = 1, 2, \dots 10$. Additionally, the controls $\hat{u_i}$ for $i = 1, 2, 3$ are given by the following:

Proof. The adjoint system in Eq (4.3) is taken from the following formula

with the transversality condition $\lambda_i(T) = 0$ for $i = 1, 2\dots 10$. To obtain the control characterization, we first take the derivative of $\mathcal{H}$ with respect to each control variable. Hence, we obtain the following:

Solve $\frac{\partial \mathcal{H}}{\partial u_i} = 0$ for each control variables, we have the following:

Hence, based on the standard argument that $\hat{u_i}$ should be bounded between $u_i^{\min}$ and $u_i^{\max}$ (in this study, we choose $u_i^{\min} = 0, u_i^{\max} = 1$ for each control variable), we have Eq (4.4).

4.2.   Numerical experiments

In this section, we perform numerical simulations for the optimal control problem. To summarize, we aim to minimize the cost function in Eq (2.2) subject to the yellow fever model in Eq (2.1) as the state system, and the adjoint system in Eq (4.3) with the terminal condition given, and the optimality condition in Eq (4.4). We performed numerical simulations using the forward-backward sweep method ^[50]. To perform these simulation, we chose the following initial conditions:

Additionally, to balance each component in the cost function Eq (2.2), we choose $\omega_3 = \omega_4 = \omega_5 = \omega_6 = \omega_7 = 1,$ $\omega_{10} = 0.1, \varphi_1 = 2500, \varphi_2 = 5000,$ and $\varphi_3 = 1500$. We use the same parameter values as in the previous section, which gives $\mathscr{R}_0 = 1.504$. Hence, without further improvement in the intervention, yellow fever will be endemic to the population. The autonomous simulation of system (3.1) without any intervention $(u_1(t) = u_2(t) = u_3(t) = 0)$ is shown as a red curve in Figures 9–13. We called this simulation as Strategy 0.

For optimal control experiments, we divide our experiments based on the combination of each form of intervention as follows:

● Strategy 1: Vaccination only $(u_1(t) \geq 0, u_2(t) = u_3(t) = 0)$.

● Strategy 2: Fumigation only $(u_1(t) = u_2(t) = 0, u_3\geq 0)$.

● Strategy 3: Vaccination and hospitalization only $(u_1(t)\geq 0, u_2(t)\geq 0, u_3(t) = 0)$.

● Strategy 4: Vaccination and fumigation only $(u_1(t)\geq 0, u_2(t) = 0, u_3(t)\geq 0)$.

● Strategy 5: All interventions implemented $(u_1(t)\geq 0, u_2(t)\geq 0, u_3(t)\geq 0)$.

Strategy 1

In the first simulation, the optimal control simulation used vaccination as the only means of intervention. The results are shown in Figure 9.

The implementation of vaccination with the greatest effort should be done at the beginning of simulation. It will reduce the number of infected people, reducing both direct contact between vulnerable mosquitoes and infected people and the number of infected mosquitoes. Reducing the level of intervention will cause the number of infected people and mosquitoes to increase near the end of the simulation. Yellow fever vaccines can cause side effects in rare cases; therefore, the total number of infected people rose immediately after vaccination before slowly decreasing toward the end.

Based on maximum vaccination at the start, this strategy has a cost of $632,377$ and prevents $444,868$ infections, leading to a total of $2,981,365$ recovered individuals.

Strategy 2

The optimal control strategy for this simulation only uses fumigation to control mosquito populations with the aim of breaking the yellow fever infection cycle. The results are shown in Figure 10.

This strategy involves maximum fumigation for nearly three months, before drastically decreasing the rate of fumigation toward the end owing to the rapidly increasing number of infected mosquitoes since the start. This significantly reduced the number of infected people and mosquitoes. Fumigation helps control the number of infected mosquitoes that can come into direct contact with susceptible individuals.

Although this strategy has a cost of $333,345$, which is almost half that of the first strategy, it prevents $784,547$ infections, far more than the first strategy. With this strategy, $1,819,666$ people will recover by the end of the period.

Strategy 3

The third optimal control simulation combines vaccination and treatment. People in $T_1$ receive treatment to aid in the recovery process. Treatment is given because yellow fever vaccines may cause side effects in rare cases. The results are shown in Figure 11.

The implementation of vaccination with the greatest effort should be done at the beginning of simulation. It would be accompanied by the side effects of the yellow fever vaccine. Treatment must counteract this effect. Treatment continued until day 14. This reduces the number of infected individuals, which without any intervention will increase until the 30$^{th}$ day. Vaccination is reduced because of the declining number of infected people.

This strategy has a higher cost than the previous two strategies ($664,808$) but is only slightly better at preventing infections ($458,867$) and yields a similar number of recovered people ($2,949,196$).

Strategy 4

The fourth optimal control simulation uses both vaccination and fumigation to combat yellow fever. The results are shown in Figure 12.

Maximum fumigation and vaccination are given simultaneously for the first 60 days to slow the rise in the number of infected people and mosquitoes.

Maximum intervention throughout the period is impossible because of cost constraints. Therefore, vaccination is reduced after the 65$^{th}$ day, and fumigation is reduced after the 81$^{st}$ day. This strategy causes the number of infected mosquitoes and infected humans to approach zero on the 14$^{th}$ and the 38$^{th}$ days, respectively.

This strategy has the highest cost of all four strategies ($1,933,426$) but yields the same number of prevented infections as the second strategy ($815,668$) and fewer recovered individuals than the second strategy ($1,551,319$).

Strategy 5

In the last simulation, the optimal control simulation simultaneously involved all interventions (vaccination, hospitalization, and fumigation) to control the spread of yellow fever. The simulation results are shown in Figure 13. Note that all interventions must be given to the fullest extent possible from the start of the simulation. Therefore, the number of infected humans and mosquitoes has decreased significantly since the first day. When the infected population in both human and mosquito populations stops increasing, the vaccination, hospitalization, and fumigation interventions will gradually decrease until they reach zero on the final day of the simulation.

With maximum intervention from the beginning of the simulation, the cost function is worth $2694.146$, which is much larger than the cost function value in Strategy 4. Nonetheless, when compared to Strategy 4, the number of infected humans successfully avoided was not significantly different, that is, 816 individuals. Similarly, the number of successfully recovered people is not significantly different from Strategy 4, which is 1551 individuals.

4.3.   Cost-effectiveness analysis

It is critical to determine the most cost-effective strategy among the five proposed strategies based on the results of our optimal control simulations. To accomplish this, we use three methods: the average cost-effectiveness ratio (ACER), infected averted ratio (IAR), and incremental cost-effectiveness ratio (ICER) ^[51,52,53].

Average cost-effectiveness ratio (ACER)

The ACER method calculates the average cost that should be spent for each number of avoided infected individuals. The ACER formula is as follows:

Therefore, a smaller ACER indicates a better strategy. The results of the ACER values for each strategy are shown in Figure 14. Note that the best strategy based on the ACER index implemented fumigation as a single form of intervention (Strategy 2), followed by Strategies 1, 3, 4, and 5. This means that using fumigation as a single form of intervention is the cheapest way, on average, to avoid many newly infected individuals. An interesting result shows that in our numerical experiment, vaccination is not cost-effective in reducing the expected number of newly infected individuals. This is because in our simulation, the number of infected individuals is already high (13% in the human population and 15% in the mosquito population). Furthermore, the simultaneous implementation of all strategies is not cost-effective because we must increase the already high hospitalization cost further.

Infected averted ratio (IAR)

The IAR method calculates the ratio between the number of infected averted individuals and the number of recovered individuals. The equation for IAR is as follows:

Based on the above equation, a larger IAR indicates a better strategy. Figure 15 shows that a combination of all means of interventions (Strategy 5) is the best strategy in terms of IAR, but it is only slightly different from Strategy 4. This is because when vaccination, hospitalization, and fumigation are implemented together, all means of intervention aim to reduce the endemic situation and avoid newly infected individuals. However, as we already discussed in the previous section, Strategy 5 comes at a high cost. The second-best strategy in terms of IAR is Strategy 4, followed by Strategies 2, 3, and 1. Note that implementing vaccination as a single single form of intervention is the least cost-effective strategy in the terms of IAR. This is because this type of intervention is more concerned with prevention rather than endemic reduction.

Incremental cost-effectiveness ratio (ICER)

The ICER was the final cost-effectiveness analysis method. The ICER compares the cost differences between two strategies and the number of infected averts. Hence, the equation for ICER is as follows:

Based on the numerical simulation, we rank all strategies in an increasing order based on the total number of infections averted in Table 3. Strategy 1 (intervention of vaccination only) has the lowest number of infections averted, while Strategy 5 (all interventions implemented) has the largest number of infections averted.

To compare strategies 1 and 3, the ICER is calculated as follows:

Note that ICER-3 $>$ ICER-1. Hence, we can exclude ICER-3 from the next calculation because it is more costly than ICER-3. Next, we continue our calculation by comparing the ICER between Strategies 1 and 2. The result of calculating ICER using the same method as before is shown in Table 4.

Table 4 shows that ICER-1$>$ICER-2, which means that Strategy 1 is more costly than Strategy 2. Hence, we exclude Strategy 1 from the next calculation. Up to this step, we have three last candidates for the best strategies, i.e., Strategies 2, 4, and 5. Next, we compare Strategies 2 and 4, the results of which are shown in Table 5.

Table 5 shows that Strategy 4 is more costly compared to Strategy 2 because ICER-4 $>$ ICER-2. Hence, we exclude Strategy 4 from the next comparison. The results are presented in Table 6.

Table 6 indicates that ICER-5 $>$ ICER-2, which means that Strategy 5 (all interventions implemented) is more costly compared to Strategy 2 (only fumigation strategy implemented) in reducing the spread of yellow fever. Additionally, Strategy 5 is dominated by Strategy 2, which means that Strategy 5 is less effective compared to Strategy 2. Hence, we can conclude that fumigation, as a single intervention to reduce the spread of yellow fever, is the most cost-effective strategy compared to other possible strategies.

5.   Conclusions

Yellow fever is an acute viral hemorrhagic fever caused by the yellow fever virus of the Flavivirus genus. It is spread by the bite of several types of female mosquitoes (Aedes, Haemagogus) ^[54]. It is estimated that more than 100,000 cases occur each year, with Africa accounting for more than 60% of all cases. There are three types of yellow fever based on species that interact during the transmission process. The first type is the jungle transmission, where the transmission process involves monkeys and mosquitoes, while the second type, namely, the urban type, involves mosquitoes and humans. The last type, namely the intermediate type, which commonly occurs in Africa, involves mosquitoes, monkeys, and humans in the transmission process ^[55]. Yellow fever has an incubation period of three to six days, and the symptoms vary, from headache, backache, vomiting, bleeding, and yellow eyes to death ^[56]. Several means of interventions have been implemented to control the spread of yellow fever, such as vaccination and fumigation. Vaccination for yellow fever is strongly recommended for individuals planning to travel to the yellow fever endemic area. However, there are a few cases where yellow fever vaccination has a side effect on some vulnerable individuals, such as the elderly, pregnant women, etc.^[25]. Infected individuals who already have severe symptoms should be hospitalized.

This work introduces a novel optimal control of the yellow fever model based on a ten-dimensional system of ordinary differential equations with a compartmental model approach. Our model includes three possible form of interventions: vaccination, hospitalization, and fumigation. We show that our model is mathematically and biologically well-defined. Our mathematical analysis shows that our model has two possible equilibriums: the yellow fever-free equilibrium point and the yellow fever endemic equilibrium point. We show that the basic reproduction number of our model plays an important role in determining the existence and stability of our equilibrium points. We discovered that the yellow fever-free equilibrium is always locally asymptotically stable when the basic reproduction number is less than one and unstable when it is greater than one. However, the yellow fever endemic equilibrium never exists when the basic reproduction number is less than one and exists (and is stable) when the basic reproduction number is greater than one. Additionally, we found that the endemic size in human and mosquito populations increased rapidly when the basic reproduction number increased. Using the Castillo-Song bifurcation theorem ^[41], we demonstrated that our model exhibits a forward bifurcation phenomenon at a basic reproduction number of one. Hence, we can conclude that our basic reproduction number becomes a necessary threshold to determine whether yellow fever will be endemic or disappear from the population.

An optimal control problem is then developed and characterized using the Pontryagin maximum principle ^[49]. We treated our form of interventions (vaccination, hospitalization, and fumigation) as time-dependent variables. We ran our simulations using a forward-backward method ^[50] for five types of strategies, depending on the possible combination of strategies. The simulation results demonstrated that all means of intervention can eradicate yellow fever. Our sensitivity analysis found that fumigation was the most elastic parameter for determining the basic reproduction number. A cost-effectiveness analysis was conducted to determine the best strategies to prevent the spread of yellow fever. Three types of cost-effectiveness indicators were used: ACER, IAR, and ICER. We found that vaccination intervention showed the best ACER results, implying that the vaccination strategy as a single intervention provides the lowest average cost required for an averted infection individual. On the other hand, the strategy of combining vaccination, hospitalization, and fumigation showed the best IAR results. Additionally, we found that using fumigation as a single form of intervention shows the best ICER result, which agrees with the ACER result. Therefore, the results of this study can be used by authorities to determine the best strategy to implement in the field to control the spread of yellow fever.

This study considers only three form of interventions: fumigation, vaccination, and hospitalizations. However, some other means of intervention can be chosen as an option, such as using bed nets or mosquito repellents to reduce the mosquito bite. Furthermore, massive vector control such as fumigation, is a high risk because it can trigger mosquitoes' resistance to fumigants. Hence, a mathematical model concerning more means of intervention and the possibility of mosquito resistance could be considered in future work. Last but not least, it is essential to fit our model to the real-life situation by estimating the best fit parameter. Hence, incidence data of yellow fever is essential for conducting parameter estimation into our model.

Acknowledgements

The author would like to thank the anonymous reviewers for his/her valuable comments. This research was funded by the Indonesian RistekBRIN with PDUPT research grant scheme (ID Number: NKB-159/UN2.RST/HKP.05.00/2021), 2021.

Conflict of interest

The authors declare that they have no conflict of interest.

Appendix

A.   Proof of Theorem 1

From yellow fever model in Eq (3.1), with a non-negative initial conditions on each variables, we have the following result:

The above calculation at the boundary of each variable indicates that the rates of each variable are non-negative on the boundary planes of $\mathbb{R}_+^{10}$. Hence, we can conclude that all vector fields in the boundary are always pointed inward. Therefore, whenever the initial condition is non-negative in $\mathbb{R}_+^{10}$, the solution will always remain non-negative for all times $t > 0$. Hence, the proof is complete.

B.   Proof of Theorem 2

Because $\frac{dN}{dt} = \Lambda_h - I \alpha_1 - T_1 \alpha_2 - T_2 \alpha_3 - \mu_h N$, we have

Solve the above equation with respect to $N(t)$, we have the following:

Therefore, for $t\to \infty$, $N(t)$ tends to $\dfrac{\Lambda_h}{\mu_h}$. A similar approach can be used to demonstrate that the total mosquito population $(M(t))$ is eventually bounded by $\dfrac{\Lambda_v}{\mu_v}$. Hence, the proof is complete.

C.   Derivation of $\mathscr{R}_0$

Taking only the infected compartment in system (3.1), the transmission $({\textbf{T}})$ and the transition $(\Sigma)$ matrix, which are evaluated in $\mathscr{E}_1$, are given by the following:

and

where $a_{22} = -(1-p) \gamma_2 - p \gamma_2 - \mu_h.$ Using the formula in ^[33], the next generation matrix of system (3.1) is given by ${\textbf{K}} = -E^t {\textbf{T}} \Sigma^{-1} E$, where

Hence, the basic reproduction number of system (3.1) as the spectral radius of K is given by the following:

with

D.   The proof of Theorem 3

Let us simplify the notion of $E_1, E_2, I, T_1, T_2, W, R, S, V, U$ as $x_i$ for $i = 1, 2, \dots 10$, respectively. We redefine system (3.1) as follows:

Now, we define our yellow-fever free equilibrium as

To apply Theorem 2 in ^[34], we have to show that the system (D.1) satisfies five axioms in ^[34]. Therefore, we need to separate the infection $(\mathcal{F})$ and the non-infection $(\mathcal{V})$ term in the symstem (D.1) as follows:

Note that $\mathcal{V}_i$ can be decomposed as $\mathcal{V}_i = \mathcal{V}_i^{-} - \mathcal{V}_i^+$, where $\mathcal{V}_i^{-}$ and $\mathcal{V}_i^+$ represent the out and in flow in compartment-$i$, respectively. Now, we are ready to proof the five axioms in ^[34].

1. If all compartment ($x_i$, $i = 1, 2, \dots 10$) are non-negative, then $\mathcal{F}_i, \mathcal{V}_i^{-},$ and $\mathcal{V}_i^+$ are always non-negative. By substituting $x_i\geq 0$ into $\mathcal{F}_i, \mathcal{V}_i^-,$ and $\mathcal{V}_i^+$, we can verify that $\mathcal{F}_i, \mathcal{V}_i^-,$ and $\mathcal{V}_i^+$ are always non-negative.

2. If all compartments are empty, then there is no out-flow rate from each compartment. It is easy to verify that, $\mathcal{V}_i^- = 0$ for $i = 1, 2, \dots 10$ when we set $x_i = 0$ for $i = 1, 2, \dots 10$.

3. For all non infected compartment ($x_i$, $i = 7, 8, 9, 10$), we have that $\mathcal{F}_i = 0$. This axiom satisfies directly by the expression of $\mathcal{F}$ at (D.3).(?)

4. If $x_i \in {{\bf{X}}_{\bf{s}}}$, then $\mathcal{F}_i = 0 = \mathcal{V}_i^+ = 0$ for $i = 1, 2, \dots 6$. By substituting ${{\bf{X}}_{\bf{s}}}$ into $\mathcal{F}_i$, $i = 1, 2, \dots 6$, we have that $\mathcal{F}_i = 0$ and $\mathcal{V}_i^+ = 0$, $i = 1, 2, \dots 6$.

5. If $\mathcal{F}(x) = 0,$ then all the eigenvalues of $Df({{\bf{X}}_{\bf{s}}})$ have a negative real part, where $Df({{\bf{X}}_{\bf{s}}})$ is the Jacobian matrix of (D.1) evaluated at ${{\bf{X}}_{\bf{s}}}$. By substituting ${{\bf{X}}_{\bf{s}}}$ into system (D.1), we obtain :

Then, we can evaluate $Df({{\bf{X}}_{\bf{s}}})$ of system (D.4). By standard calculation, we obtain the eigenvalues of $Df({{\bf{X}}_{\bf{s}}})$ as follows:

Since all parameter values are non-negative, we can see that all eigenvalues of $Df({{\bf{X}}_{\bf{s}}})$ are negative.

Hence, the proof is completed.

E.   The proof of Theorem 4

From the expression of $F(W)$ in Eq (3.5), we can see that $a_2$ is always positive, whereas $a_0 < 0 \iff \mathscr{R}_0^2 > 1$. Using the formula for the multiplication of roots in a second-degree polynomial, we have $W_1^{\ddagger} \times W_2^{\ddagger} = \frac{a_0}{a_2}$, where $W_1^{\ddagger}$ and $W_2^{\ddagger}$ is the root of $F(W)$. Hence, if $\mathscr{R}_0^2 > 1$, we have $W_1^{\ddagger} \times W_2^{\ddagger} < 0$, which indicates that there exists one positive root, and the other one is negative. Therefore, we have a unique yellow fever endemic equilibrium point when $\mathscr{R}_0^2 > 1 \iff \mathscr{R}_0 > 1$.

Now, we analyze the possibility of the existence of $\mathscr{E}_2$ when $\mathscr{R}_0^2 < 1$. To accomplish this, we use the properties of the second-degree polynomial to have two positive roots, that is, $W_1^{\ddagger} \times W_2^{\ddagger} = \frac{a_0}{a_2} > 0$, $W_1^{\ddagger} + W_2^{\ddagger} = -\frac{a_1}{a_2} > 0$, and $a_1^2 - 4a_2 a_0 \geq 0$. However, because the expressions of $a_2, a_1,$ and $a_0$ is not simple, the expression of $a_1^2 - 4a_2 a_0$ is not easy to analyze. Hence, we use another approach to analyze the existence of $\mathscr{E}_2$ when $\mathscr{R}_0 < 1$. To accomplish this, we use a gradient analysis of $W$ in $\mathscr{R}_0 = 1$ and $W = 0$. If we can find a condition such that $\frac{\partial W}{\partial \mathscr{R}_0} < 0$, we can have at least one positive root of $F(W)$ when $\mathscr{R}_0 < 1$. Because the expression of $\mathscr{R}_0$ contains a square root, and the condition of $\mathscr{R}_0 = 1$ is represented by $\mathscr{R}_0^2 = 1$, we will analyze the gradient of $W$ with respect to $\mathscr{R}_0^2$ at $\mathscr{R}_0^2 = 1$ instead of $\mathscr{R}_0 = 1$.

Solving the expression of $\mathscr{R}_0^2 = 1$ in Eq (3.3) with respect to $\beta_h$, and substituting it with $F(W)$, we have $a_2, a_1$, and $a_0$, which is an expression of $\mathscr{R}_0^2$. Taking the implicit derivative of $W$ with respect to $\mathscr{R}_0^2$ from $F(W)$, we have the following:

substitute $W = 0$, and solve it with respect to $\dfrac{\partial W^\ddagger}{\partial\mathscr{R}_0^2}$, we have the following:

where $\dfrac{\partial a_0}{\partial\mathscr{R}_0^2} = -2\, \mu_{{h}}(\mu_{{v}}+u_{{3}}) ^{2} (\gamma_{{2}}+\mu_{{h}})$$(\gamma_{{1}}+\mu_{{h}})(\alpha_{{2}}+\delta_{{2}}+\mu_{{h}}+u_{{2}})(\mu_{{h}}+\omega_{{h}}+u_{{1}})(\alpha_{{1}}+\delta_{{1}}+\mu_{{h}}) < 0$ and $a_1 > 0$. Hence, the gradient of $W$ at $\mathscr{R}_0^2 = 1$ is always positive without any condition. Therefore, we can conclude that there is no yellow fever endemic equilibrium when $\mathscr{R}_0 < 1$. Hence, the proof is complete.