Effective Control Strategies for Sex-Structured Transmission Dynamics of Visceral Leishmaniasis

1. Introduction

A chronic infectious disease caused by Leishmania infantum, visceral leishmaniasis (VL), also referred to as kala-azar, is the second largest cause of death of humans in the globe next to malaria [1]. An estimated 12 million individuals are affected with the disease worldwide in 98 endemic countries, and an additional 350 million are at risk of contracting it [2]. The vector-borne disease, VL, is transmitted through the bite of a very small flying insect female adult sand-fly. Warm and humid weather is suitable for the sand flies to flourish in[3]. Different interventions have been implemented in various endemic areas of the world to control the transmission of VL. These include personal protection through the use of clothes that cover as much of the body as possible, avoiding being outside between dusk and dawn, utilizing insecticide-treated bed nets, use of indoor and outdoor spraying of insecticides, using window and door screens to prevent sand flies bites [4,5] together with WHO-recommended therapies such as liposomal amphotericin B, Miltefosine, Pentavalent antimonials, and Paromomycin [6]. However, due to the costs of medicines and the administration route, the VL control strategy continues to be challenging. Mathematical and epidemiological models that consider different intervention mechanisms have been developed to minimize such challenges.

These models are indispensable tools for understanding, predicting, and controlling the spread of infectious diseases. They provide critical insights that inform policy decisions, optimize resource allocation, support emergency response efforts, and drive innovation in disease control strategies. Hence, numerous modeling studies of VL have been conducted to understand the impact of different intervention strategies and to describe and predict how the disease spreads (see [1,7,11,12,13,14,15,16] and the reference therein). While the mathematical models of VL are valuable for understanding the disease dynamics, they do not explicitly optimize control strategies. Optimal control models provide more rigorous and efficient tools for designing and implementing effective disease control strategies. However, some studies have been conducted to investigate the optimal control strategies for controlling the spread of VL.

For example, Pantha et al. [20] developed a mathematical model for VL transmission dynamics involving human, canine, and vector populations. The goal was to identify effective control measures for reducing VL burden and guiding public health interventions in endemic areas. They evaluated three control measures using optimal control theory: personal protection, insecticide spraying, and culling infected canine reservoirs. Simulation results showed that combining these interventions could eliminate VL from the region. Zhao et al. [19] developed a mathematical model for zoonotic visceral leishmaniasis (ZVL) in Brazil using a modified SEIR-type approach, focusing on dogs as reservoirs. They classified the dog population into susceptible, exposed, infected, and recovered categories and the human and sand fly populations into similar compartments. Mathematical analysis indicated the possibility of backward bifurcation. They introduced three time-dependent controls using optimal control theory: dog vaccination, sand fly control with insecticides, and personal human protection. Numerical simulations showed that vector control via insecticide spraying was the most effective strategy, while culling dogs was ineffective. However, the study lacked validation with real data and did not include cost-effectiveness analysis.

Ghosh et al. [21] developed a compartmental VL model considering human, animal, sand fly, and non-adult sand fly populations. The human population was divided into six classes, including asymptomatic, symptomatic, and PKDL-infected individuals. Global sensitivity analysis identified key parameters affecting the cost function, and the model was validated with 2013 weekly VL incidence data from South Sudan. The model used three time-dependent controls: treatment of infected individuals, non-adult stage vector control, and adult stage vector control. Cost-effectiveness analysis showed that combining treatment and non-adult stage vector control was the most cost-effective strategy. Augusto et al. [9] developed a deterministic model for malaria-VL co-infection to understand its dynamics and find cost-effective control strategies. Using the basic model by ELmojtaba [11], they introduced four time-dependent controls: personal human protection, insecticide-treated bed-nets, culling of infected reservoir animals, and insecticides. Their results show that combining personal protection, indoor residual spraying, and reservoir culling is the most cost-effective approach for controlling malaria-VL co-infection. In several endemic regions of the world, the results of numerous research have demonstrated that males are more severely impacted by VL than females. The observed gender disparity in visceral leishmaniasis is likely influenced by a combination of behavioral, biological, socioeconomic, and environmental factors. Conducting a sex-structured disease dynamics model for visceral leishmaniasis allows for a more comprehensive understanding of transmission dynamics. To address this issue, a sex-structured VL disease transmission model was conducted by Awoke et al. [10]. However, the research work did not include optimal control and cost-effectiveness analysis. The optimal control plays a crucial role in disease dynamics models. It offers insightful information about the possible effects of various intervention strategies (e.g., insecticide-treated bed nets, vector control, and medical treatment), and supports long-term planning goals by taking the sustainability and feasibility of intervention strategies over extended time horizons into account. Furthermore, it helps public health decision-makers prioritize disease control strategies during resource allocation and maximizes the use of limited resources (for example vaccines, drugs, and medical staff). Since VL treatments are very expensive, it is very important to study sex-structured VL models that incorporate disease control mechanisms and propose cost-effective intervention strategies.

Managing VL poses a financial burden on local communities in countries where the disease is prevalent, perpetuating a cycle of illness and poverty. Strategic allocation of public funds towards control programs could alleviate this burden by reducing the prevalence of VL [22]. To this end, employing dynamic models that incorporate an optimal control theory is essential for assessing intervention strategies and devising an economically efficient policy to contain the disease. The motivation for this research arises from the considerable public health challenge posed by VL, a disease that predominantly affects males. While many VL models, including those incorporating animal reservoirs and optimal control strategies, have been studied, there remains a gap in analyzing the cost-effectiveness of interventions in sex-structured populations. Building on the VL model developed by Awoke et al. [10], this article introduces four time-dependent control interventions: the use of insecticide-treated bed nets, vector control through residual spraying and environmental management, medical treatment, and infected animal culling. To the best of the knowledge of the authors, no study has yet explored the cost-effectiveness of such a comprehensive, sex-structured approach to VL control. To evaluate whether sex differences influence optimal control strategies, we incorporated parameters for the use of insecticide-treated bed nets and medical treatment. By addressing this gap, this research will provide crucial insights that can directly benefit both local communities and policymakers. A thorough cost-effectiveness analysis will identify the most cost-effective strategies for reducing VL transmission, helping to optimize limited resources in low-income regions where the disease is prevalent. This research will also equip policymakers with evidence-based recommendations, allowing them to prioritize interventions that offer the greatest health benefits for the lowest cost, ultimately contributing to more sustainable and effective public health strategies. Hence, our study aims to integrate effective preventive and control measures for visceral leishmaniasis (VL) into the existing sex-structured model. We intend to assess the cost-effectiveness of these interventions and propose an optimal control strategy to enhance decision-making processes.

The rest of the manuscript is organized as follows: Section 2, presents the model formulation, control strategies and analysis of key properties, including positivity, boundedness and the basic reproduction number with controls. This sections also introduces the optimal control problem and its characterization. Section 3 presents results and discussion. Conclusions are given in Section 4.

2. Methods

Controlling VL remains challenging in different endemic areas of the world for different reasons, such as the complexity of the transmission, lack of effective vaccine, cost of treatment, and socioeconomic factors. In addition, the disease mainly affects the male human population. Hence, there is a need to develop comprehensive cost-effective control strategies that will address these challenges. By optimizing control measures across different gender groups and considering trade-offs between intervention options, optimal control helps maximize the impact of interventions and improve disease control efforts. The most powerful VL disease transmission control mechanisms are the use of insecticide-treated bed-nets and medical treatment of infected humans together with vector control strategies. In this work, we want to consider these interventions as time-dependent control function and incorporate them into the sex-structured VL model proposed in [10]. The total human population is initially divided into male and female sub-populations, with each further classified into five compartments: susceptible males and females ($S_m$, $S_f$), exposed males and females ($E_m$, $E_f$), asymptomatic males and females ($A_m$, $A_f$), actively infected males and females ($I_m$, $I_f$), and recovered males and females ($R_m$, $R_f$). The sand fly population is split into susceptible ($S_v$) and infected ($I_v$) groups, while the reservoir population is categorized into susceptible ($S_r$), infected ($I_r$), and recovered ($R_r$) groups. The description of the model parameters are indicated in Table 1 and for the detail please refer to [10].

As described in [23,33], the possible interventions for visceral leishmaniasis can be categorized as: applying preventive mechanisms through the use of Insecticide Treated Bed Nets (ITNs), vector control with the use of insecticides, treatment of infected individuals and animal culling. The practical applications of these mechanisms, and assumptions used in this work are described below.

• Preventive mechanisms through the use of ITNs
  The WHO recommended VL preventive strategies include the use of insecticide treated sand fly nets [24]. In order to curb the spread of VL in a study by [25], the breakthrough in the reduction of the spread of VL was achieved after the introduction of ITNs. Due to its practical applicability, we are considering the use of ITNs as one of the control strategies. Using ITNs reduces the probability of humans being bitten. Furthermore, since the nets are treated with insecticides, they kill vectors. Therefore, we used the following two assumptions formulated in [9,28].

  a)

  The mathematical expression for the use of ITNs in order to reduce human-vector contact rate is described by the relation

  $$\begin{array}{cc}\hfill & C\left(u_1\right)=C_{\max }-u_1(C_{\max }-C_{\min }),\hfill \end{array}$$

  (1)

  where $C_{\max }$ and $C_{\min }$ are the maximum and minimum possible number of bites a host population receives per week, respectively. Since most sandfly biting activity occurs during the early evening before most people go to sleep and other limitations of the use of INTs, the maximum control effort is assumed to be $65\%$. Hence, $0\le u_1\left(t\right)\le u_{1\max }\le 1$. Applying different interventions, such as treated bed nets and medication, to different sexes (male and female human populations) in the control of VL might raise ethical considerations. However, due to biological and behavioral differences as well as the need to optimize resource allocation, considering differential interventions for both sexes might trigger an ethical question. Males and females often have different exposure risks to the sand fly. In VL-endemic areas, men may work outdoors in high-risk settings where they are more exposed to sand fly bites, while women may spend more time indoors. Tailoring interventions to these behavioral patterns, such as providing stronger protection for men who are more exposed outdoors, may be a targeted and effective strategy for disease control. In addition to these, the availability of limited resources in public health programs often necessitates prioritizing interventions for populations at greater risk. Since data reveals that men are at higher risk of contracting the disease than women, focusing certain interventions like availing treated bed nets or medications on that group can maximize the impact of disease control efforts. Due to these reasons we assumed that the intervention rates for female is reduced by a proportion ${\varphi }_1$ for the use of bed nets and by ${\varphi }_2$ for the use of medical treatment as compared to that of the male counterparts, where $0<{\varphi }_1,{\varphi }_2<1$. In this work, the expression for the use of ITNs is based on the concept formulated in Asfaw et al.[28]. Therefore, the forces of infection for the male human population (${\lambda }_m=c_{1m}{\beta }_{vh}{\displaystyle \frac{I_v}{N_m}}$), female human population (${\lambda }_f=c_{1f}{\beta }_{vh}{\displaystyle \frac{I_v}{N_f}}$), and vector population (${\lambda }_v=c_{3m}{\displaystyle \frac{{\beta }_{hv}}{N_m}}({\kappa }_1\zeta A_m+{\kappa }_2{I}_m)+c_{3f}{\displaystyle \frac{{\beta }_{hv}}{N_f}}(\zeta A_f+I_f)+c_4{\beta }_{rv}{\displaystyle \frac{I_r}{N_r}}$) are re-written as follows:

  $$\begin{array}{cc}\hfill & {\lambda }_m=(C_{\max }-u_{1}C){\beta }_{vh}{\displaystyle \frac{I_v}{N_m}},{\lambda }_f=(C_{\max }-{\varphi }_1{u}_{1}C){\beta }_{vh}{\displaystyle \frac{I_v}{N_f}},\hfill \\ \hfill & {\lambda }_v=(C_{\max }-u_{1}C){\displaystyle \frac{{\beta }_{hv}}{N_m}}({\kappa }_1\zeta A_m+{\kappa }_2{I}_m)+c_4{\beta }_{rv}{\displaystyle \frac{I_r}{N_r}}\hfill \end{array}$$

  (2)

  $$\begin{array}{c}\hfill +(C_{\max }-{\varphi }_1{u}_{1}C){\displaystyle \frac{{\beta }_{hv}}{N_f}}(\zeta A_f+I_f),{\lambda }_r=c_2{\beta }_{vr}{\displaystyle \frac{I_v}{N_r}},\end{array}$$

  (3)

  where $C=C_{\max }-C_{\min }$. Our model is based on the assumption that vector bites are distributed in fixed proportions among males, females, and reservoir hosts. However, in real-world scenarios, the use of bed nets by a subset of the human population could alter this distribution. Specifically, when some individuals are protected by bed nets, vectors may shift their biting behavior toward unprotected groups, potentially leading to increased exposure among unprotected groups such as reservoir hosts. We acknowledge that this adaptive biting behavior is not captured in our current model.

  b)

  Since bed nets are treated with insecticides, they have the capacity of killing vectors that landed on the net. Hence, if the natural death rate of vectors is ${\mu }_v$, ITNs usage increases the sand fly death rate by ${\alpha }_1{u}_1$. This is mathematically represented by

  $$\begin{array}{c}\hfill {\mu }_v\left(u_1\right)={\mu }_v+{\alpha }_1{u}_1,\end{array}$$

  (4)

  where ${\alpha }_1$ represents the capability of insecticide-treated bed-nets to kill sand flies and ${\alpha }_1{u}_1$ represents the death rate of adult sand flies due to insecticide-treated bed nets.

• Vector control: The incidence rate of VL is increasing in some endemic regions of the world, including Kenya, Ethiopia, and Sudan [32], as a result of the serious side effects of the current VL medications and the lack of an effective vaccination. Programs to eradicate leishmaniasis, therefore, mainly rely on vector control using insecticides and environmental management. The goal of environmental management is to eliminate or make places undesirable for sand flies to nest and rest, as well as to kill the pre-adult stage of vectors. Tree stump clearance, soil and crack filling, routine cleaning of the per-domicile area and animal shelters, and trash removal are common vector control methods through environmental management [31]. Additionally, the most effective method of controlling sand flies is the use of chemical interventions such as indoor residual spraying (IRS), which is associated with a decrease in the lifespan of sand flies, leading to a subsequent decline in their population [14,26,27]. By reducing the lifespan of sand flies, the vector is less likely to survive long enough to acquire infection and to infect a host again. As the death rate of vectors increases, they spend less time in an infected state. Taking this into account, we incorporate a time-dependent control $u_2$ that represent an additional vector control effort through the use of environmental management and indoor residual spraying into the VL model. If the current natural mortality rate of vectors is ${\mu }_v$, then the number of additional deaths of non-adult sand flies after the use of vector control intervention will be ${\alpha }_2{u}_2$, where ${\alpha }_2$ represents the effectiveness of insecticide spraying in killing sand flies. This implies that ${\mu }_v$ of Eq.(4) will be updated by an additional term ${\alpha }_2{u}_2$. Such that

  $$\begin{array}{c}\hfill {\mu }_v(u_1,u_2)={\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2.\end{array}$$

  (5)
  Furthermore, the recruitment rate of vectors (${\Lambda }_v$) will decrease by $1-u_2$ as the application of insecticides is expected to destroy the breeding cites and kills the non-adult stage sand flies. Hence, ${\Lambda }_v$ is updated by $(1-u_2){\Lambda }_v$, where $0\le u_2\left(t\right)\le u_2\max \le 1$. Note here that $u_1$ and $u_2$ must satisfy the condition that $0\le {\mu }_v(u_1,u_2)\le 1.$
• Treating infected individuals: Treatment is an important method to recover from the disease. Treatment has an effective role in preventing and controlling the spreading of the disease. Proper and timely treatment can reduce the total prevalence of the disease. However, treatment alone will not be enough to truly halt the transmission of VL [29,30]. The efficacy of VL drugs varies from one endemic area to another. We assume that the efficacy of VL drug is $ϵ$. If ${\tau }_1$ and ${\tau }_2$ are the recovery rates of male and female human host, then the mathematical expression for the formulation of controls with an additional treatment effort ($u_3$) is given by: ${\tau }_1={\tau }_{10}+ϵu_3$ and ${\tau }_2={\tau }_{20}+{\varphi }_2ϵu_3$. Where ${\tau }_0$ represents the natural recovery rate of the male human host, ${\tau }_0$ represents the natural recovery of female human hosts, and ${\varphi }_2$ is a modification parameter, as stated before, with  $0\le u_3\left(t\right)\le 1$ for every time $t>0$.
• Culling infected animals: Dogs are the most important reservoirs of VL infection in domestic settings, while other animals can serve as a source of VL infection [33]. As a result, it is believed that having infected dogs increases the risk of VL in people in endemic areas. Culling leishmania-seropositive reservoirs such as dogs has been recommended as a control measure in many VL endemic countries [34]. It is indicated in [20] and [9] that, culling of infected animals is a cost effective VL control strategy together with other interventions. Hence, we consider culling of infected dogs as one of the control strategy. Given ${\mu }_r$ and ${\delta }_r$, which represent the natural and disease-induced mortality rates of infected dogs, ${\mu }_r+{\delta }_r+u_4$ provides the mathematical expression for the rate at which infected dogs diminish as a result of further culling effort $u_4$, where $0\le u_4\left(t\right)\le u_{4\max }\le 1$ for every time $t>0$.

In general, by considering all these, we incorporate the above four VL intervention strategies in the VL model formulated and analyzed in [10]. Since most governments cannot implement intervention programs indefinitely, we considered the optimal control problem with a fixed terminal time $[0,T]$, where T is the final time. The corresponding system of differential equations can be written as in Eq.(6) once these control parameters have been applied to the model system.

Figure 1. Schematic diagram for VL compartmental model in which the human population is classified based on their sex.

Figure 1. Schematic diagram for VL compartmental model in which the human population is classified based on their sex.

$$F(\mathbf{x},\mathbf{u},t)=\left\{\begin{array}{cc}{\dot{S}}_m=\hfill & {\Lambda }_m-({\lambda }_m+{\mu }_h)S_m+{\psi }_1{R}_m\hfill \\ {\dot{E}}_m=\hfill & {\lambda }_m{S}_m-({\mu }_h+\eta )E_m\hfill \\ {\dot{A}}_m=\hfill & p\eta E_m-({\theta }_1+{\omega }_1+{\mu }_h)A_m\hfill \\ {\dot{I}}_m=\hfill & (1-p)\eta E_m+{\theta }_1{A}_m-\left({\tau }_{10}+\epsilon u_3+{\delta }_1+{\mu }_h\right)I_m\hfill \\ {\dot{R}}_m=\hfill & {\omega }_1{A}_m+\left({\tau }_{10}+\epsilon u_3\right)I_m-({\mu }_h+{\psi }_1)R_m\hfill \\ {\dot{S}}_f=\hfill & {\Lambda }_f-({\lambda }_f+{\mu }_h)S_f+{\psi }_2{R}_f\hfill \\ {\dot{E}}_f=\hfill & {\lambda }_f{S}_f-({\mu }_h+\eta )E_f\hfill \\ {\dot{A}}_f=\hfill & q\eta E_f-({\theta }_2+{\omega }_2+{\mu }_h)A_f\hfill \\ {\dot{I}}_f=\hfill & (1-q)\eta E_f+{\theta }_2{A}_f-\left({\tau }_{20}+\epsilon {\varphi }_2{u}_3+{\delta }_2+{\mu }_h\right)I_f\hfill \\ {\dot{R}}_f=\hfill & {\omega }_2{A}_f+\left({\tau }_{20}+\epsilon {\varphi }_2{u}_3\right)I_f-({\mu }_h+{\psi }_2)R_f\hfill \\ {\dot{S}}_v=\hfill & (1-u_2){\Lambda }_v-\left({\lambda }_v+{\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2\right)S_v\hfill \\ {\dot{I}}_v=\hfill & {\lambda }_v{S}_v-\left({\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2\right)I_v\hfill \\ {\dot{S}}_r=\hfill & {\Lambda }_r-({\lambda }_r+{\mu }_r)S_r+{\psi }_r{R}_r\hfill \\ {\dot{I}}_r=\hfill & {\lambda }_r{S}_r-({\mu }_r+{\delta }_r+\gamma +u_4)I_r\hfill \\ {\dot{R}}_r=\hfill & \gamma I_r-({\mu }_r+{\psi }_r)R_r\hfill \end{array}\right.$$

(6)

where $\mathbf{x}=(x_1,x_2,\cdots ,x_{15})=(S_m,E_m,A_m,I_m,R_m,S_f,E_f,A_f,I_f,R_f,S_v,I_v,S_r,I_r,R_r)$ and $\mathbf{u}=(u_1,u_2,u_3,u_4)$. with ${\lambda }_m$, ${\lambda }_f$, ${\lambda }_v$ and ${\lambda }_r$ are as indicated in equations (2) and (3). $N_m=S_m+E_m+A_m+I_m+R_m$, $N_f=S_f+E_f+A_f+I_f+R_f$, $N_r=S_r+I_r+R_r$, and $N_v=S_v+I_v$. The model system (6) is appended with initial conditions:

$$\begin{array}{cc}\hfill & S_{m0}>0,E_{m0}\ge 0,A_{m0}\ge 0,I_{m0}\ge 0,R_{m0}\ge 0,S_{f0}>0,E_{f0}\ge 0,A_{f0}\ge 0,\hfill \\ \hfill & I_{f0}\ge 0,R_{f0}\ge 0,S_{v0}>0,I_{v0}\ge 0,S_{r0}>0,I_{r0}\ge 0,\&R_{r0}\ge 0.\hfill \end{array}$$

(7)

2.1. Model Analysis

2.1.1. Boundedness and Positivity of Solutions with Constant Control Parameters

For the VL model system (6) to be epidemiologically meaningful, it is important to analyze that all its state variables are non-negative for all time and the region $\Omega$ defined in Theorem 1 is positively invariant.

Theorem 1.

The VL model system (6) is a dynamical system on the biologically feasible region Ω. Where,

$$\begin{array}{cc}\hfill \Omega =& \{(S_m,E_m,A_m,I_m,R_m,S_f,E_f,A_f,I_f,R_f,S_v,I_v,S_r,I_r,R_r)\in {\mathbb{R}}_{+}^{15}:S_m+E_m+A_m+I_m\hfill \\ \hfill & +R_m\le {\displaystyle \frac{{\Lambda }_m}{{\mu }_h}},S_f+E_f+A_f+I_f+R_f\le {\displaystyle \frac{{\Lambda }_f}{{\mu }_h}},S_v+I_v\le {\displaystyle \frac{{\Lambda }_v}{{\mu }_v}},S_r+I_r+R_r\le {\displaystyle \frac{{\Lambda }_r}{{\mu }_r}}\}.\hfill \end{array}$$

Proof. 

The proof of Theorem 1 proceeds in a similar manner to the proof of Theorem 1 in Awoke et al. [10]. For brevity, we omit the details here, as they follow similar reasonings and steps presented in [10].    □

2.1.2. Disease Free Equilibrium Point (DFE)

The model system (6) has a disease-free equilibrium (DFE) given by

$$\begin{array}{c}\hfill {\mathcal{E}}_0=\left({\displaystyle \frac{{\Lambda }_m}{{\mu }_h}},0,0,0,0,{\displaystyle \frac{{\Lambda }_f}{{\mu }_h}},0,0,0,0,{\displaystyle \frac{{\Lambda }_v}{{\mu }_v}},0,{\displaystyle \frac{{\Lambda }_r}{{\mu }_r}},0,0\right)\end{array}$$

(8)

2.1.3. Basic Reproduction Number with Constant Controls

As indicated in [35], the basic reproduction number (${\mathcal{R}}_0$) is defined as the average number of secondary infections that occur when one infective is introduced into a completely susceptible host population.

We calculate the basic reproduction number, ${\mathcal{R}}_0$ of the model system (6), using the next generation matrix approach proposed by [35]. To apply this method, we define $\mathcal{F}$ as the column-vector of rates of the appearance of new infections in each compartment, ${\mathcal{V}}^{+}$ is the column-vector of rates of transfer of individuals into the particular compartment; and ${\mathcal{V}}^{-}$ is the column-vector of rates of transfer of individuals out of the particular compartment. Then, $\dot{x}=\mathcal{F}-\mathcal{V}$ for each compartment, where, $\mathcal{V}={\mathcal{V}}^{-}-{\mathcal{V}}^{+}$. Considering the system (6) has eight infected classes, namely $E_m, E_f, I_v, I_r, A_m, A_f, I_m$ and $I_f$ we obtain

$$\mathcal{F}\left(x\right)=\left(\begin{array}{c}{\beta }_{vh}{\displaystyle \frac{I_v}{N_m}}{S}_m(C_{\max }-u_{1}C)\\ {\beta }_{vh}{\displaystyle \frac{I_v}{N_f}}{S}_f(C_{\max }-{\varphi }_1{u}_{1}C)\\ S_v(C_{\max }-u_{1}C){\displaystyle \frac{{\beta }_{hv}}{N_m}}\left({\kappa }_1\zeta A_m+{\kappa }_2{I}_m\right)+c_4{\beta }_{rv}{S}_v{\displaystyle \frac{I_r}{N_r}}\\ +S_v(C_{\max }-{\varphi }_1{u}_{1}C){\displaystyle \frac{{\beta }_{hv}}{N_f}}\left(\zeta A_f+I_f\right)\\ c_2{\beta }_{vr}{\displaystyle \frac{I_v}{N_r}}{S}_r\\ 0\\ 0\\ 0\\ 0\end{array}\right)$$

(9)

$$\mathrm{and}\mathcal{V}\left(x\right)=\left(\begin{array}{c}(\eta +{\mu }_h)E_m\\ (\eta +{\mu }_h)E_f\\ ({\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2)I_v\\ ({\mu }_r+\gamma +{\delta }_r+u_4)I_r\\ ({\theta }_1+{\omega }_1+{\mu }_h)A_m-p\eta E_m\\ ({\theta }_2+{\omega }_2+{\mu }_h)A_f-q\eta E_f\\ ({\tau }_{01}+{\delta }_1+{\mu }_h+ϵu_3)I_m-(1-p)\eta E_m-{\theta }_1{A}_m\\ ({\tau }_{02}+{\delta }_2+{\mu }_h+{\varphi }_2ϵu_3)I_f-{\theta }_2{A}_f-(1-q)\eta E_f\end{array}\right).$$

(10)

Evaluating the Jacobian of (9) at ${\mathcal{E}}_0$ gives

$$F=\left(\begin{array}{cccccccc}0& 0& f_{13}& 0& 0& 0& 0& 0\\ 0& 0& f_{23}& 0& 0& 0& 0& 0\\ 0& 0& 0& f_{34}& f_{35}\zeta {\kappa }_1& \zeta f_{36}& {\kappa }_2{f}_{35}& f_{36}\\ 0& 0& c_2{\beta }_{vr}& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right),$$

where

$$\begin{array}{cc}\hfill & f_{13}={\beta }_{vh}(C_{\max }-u_{1}C),f_{23}={\beta }_{vh}(C_{\max }-{\varphi }_1{u}_{1}C),f_{34}={\displaystyle \frac{c_4{\beta }_{rv}{\mu }_r{\Lambda }_v}{{\Lambda }_r({\mu }_v+{\alpha }_1{u}_1+\alpha u_2)}},\hfill \\ & f_{35}={\displaystyle \frac{{\beta }_{hv}{\mu }_h{\Lambda }_v(C_{\max }-u_{1}C)}{{\Lambda }_m({\mu }_v+{\alpha }_1{u}_1+\alpha u_2)}},f_{36}={\displaystyle \frac{{\beta }_{hv}{\mu }_h{\Lambda }_v(C_{\max }-{\varphi }_1{u}_{1}C)}{{\Lambda }_f({\mu }_v+{\alpha }_1{u}_1+\alpha u_2)}}.\hfill \end{array}$$

Similarly, the Jacobian of (10) with respect to $E_m, E_f, I_v, I_r, A_m, A_f, I_m$ and $I_f$ at ${\mathcal{E}}_0$ gives

$$V=\left(\begin{array}{cccccccc}{v}_{11}& 0& 0& 0& 0& 0& 0& 0\\ 0& v_{22}& 0& 0& 0& 0& 0& 0\\ 0& 0& v_{33}& 0& 0& 0& 0& 0\\ 0& 0& 0& v_{44}& 0& 0& 0& 0\\ -p\eta & 0& 0& 0& v_{55}& 0& 0& 0\\ 0& -q\eta & 0& 0& 0& v_{66}& 0& 0\\ -(1-p)\eta & 0& 0& 0& -{\theta }_1& 0& v_{77}& 0\\ 0& -(1-q)\eta & 0& 0& 0& -{\theta }_2& 0& v_{88}\end{array}\right),$$

where $v_{11}=v_{22}=\eta +{\mu }_h$,$v_{33}={\mu }_v+{\alpha }_1{u}_1+\alpha u_2$, $v_{44}={\mu }_r+{\delta }_r+\gamma +u_4$, $v_{55}={\theta }_1+{\omega }_1+{\mu }_h$, $v_{66}={\theta }_2+{\omega }_2+{\mu }_h$, $v_{77}={\tau }_{01}+{\delta }_1+{\mu }_h+\epsilon u_3$, $v_{88}={\tau }_{02}+{\delta }_2+{\mu }_h+\epsilon {\varphi }_2{u}_3$.

The basic reproduction number (${\mathcal{R}}_0$) is then given by the spectral radius of the next generation matrix, $F{V}^{-1}$ and it is found to be

$${\mathcal{R}}_0=\sqrt{{\mathcal{R}}_{0r}+{\mathcal{R}}_{0m}+{\mathcal{R}}_{0f}},$$

(11)

where

$$\begin{array}{cc}\hfill {\mathcal{R}}_{0r}=& {\displaystyle \frac{c_2{c}_4{\beta }_{vr}{\beta }_{rv}{\mu }_r{\Lambda }_v}{{\Lambda }_r{({\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2)}^2({\mu }_r+{\delta }_r+\gamma +u_4)}},\hfill \\ \hfill {\mathcal{R}}_{0m}=& {\displaystyle \frac{\eta {\mu }_h{\beta }_{vh}{\beta }_{hv}{\Lambda }_v{(C_{\max }-u_{1}C)}^2}{{\Lambda }_m(\eta +{\mu }_h){({\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2)}^2}}\left({\displaystyle \frac{p\zeta {\kappa }_1}{{\theta }_1+{\omega }_1+{\mu }_h}}+{\displaystyle \frac{p{\theta }_1{\kappa }_2+{\kappa }_2(1-p)({\theta }_1+{\omega }_1+{\mu }_h)}{({\theta }_1+{\omega }_1+{\mu }_h)({\tau }_{10}+{\delta }_1+{\mu }_h+\epsilon u_3)}}\right),\hfill \\ \hfill {\mathcal{R}}_{0f}=& {\displaystyle \frac{\eta {\mu }_h{\beta }_{vh}{\beta }_{hv}{\Lambda }_v{(C_{\max }-{\varphi }_1{u}_{1}C)}^2}{{\Lambda }_f(\eta +{\mu }_h){({\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2)}^2}}\left({\displaystyle \frac{q\zeta }{{\theta }_2+{\omega }_2+{\mu }_h}}+{\displaystyle \frac{q{\theta }_2+(1-q)({\theta }_2+{\omega }_2+{\mu }_h)}{({\theta }_2+{\omega }_2+{\mu }_h)({\tau }_{20}+{\delta }_2+{\mu }_h+\epsilon {\varphi }_2{u}_3)}}\right).\hfill \end{array}$$

Remark 1.

The formula for basic reproduction number with control $\left({\mathcal{R}}_0\right)$ has three parts and is epidemiologically interpreted as follows.

1. 

The quantity ${\mathcal{R}}_{0r}$ is associated with the contribution of infected reservoirs,

2. 

The quantity ${\mathcal{R}}_{0m}$ is associated with the contribution of infected male humans,

3. 

The quantity ${\mathcal{R}}_{0f}$ is associated with the contribution of infected female humans in the spread of the disease at its initial stage.

Theorem 2.

The Disease free equilibrium (DFE) ${\mathcal{E}}_0$ of the model system (6) with control is locally asymptotically stable whenever ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

Proof. 

This is an immediate consequence of Theorem 2 in [35].    □

The epidemiological implication of Theorem 2 is that VL can be effectively controlled in the host and reservoir populations (humans and animal reservoirs) if the initial number of infected humans, reservoirs and vectors are small enough (i.e., in the basin of attraction of the non-trivial disease-free equilibrium, ${\mathcal{E}}_0$).

To understand how modifications in controls influence the basic reproduction number (${\mathcal{R}}_0$), we perform partial differentiation of ${\mathcal{R}}_0$ with respect to each control ($u_1$, $u_2$, $u_3$ and $u_4$), by keeping the other parameters as a constant:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial u_1}}& =-{\displaystyle \frac{{\alpha }_1}{{\mathcal{R}}_0{({\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2)}^3}}[{\mathcal{R}}_{0r}+Q_1(C_{\max }-{\varphi }_1{u}_{1}C)(C({\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2)\hfill \\ \hfill & +(C_{\max }-{\varphi }_1{u}_{1}C))+Q_2(C_{\max }-u_{1}C)\left((C_{\max }-u_{1}C)+{\displaystyle \frac{C}{{\alpha }_1}}({\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2)\right)]<0,\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial u_2}}& =-{\displaystyle \frac{{\alpha }_2}{{\mathcal{R}}_0({\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2)}}\left({\mathcal{R}}_{0r}({\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2)+{\mathcal{R}}_{0m}+{\mathcal{R}}_{0f}\right)<0,\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial u_3}}& =-{\displaystyle \frac{ϵ\eta {\mu }_h{\beta }_{vh}{\beta }_{hv}{(C_{\max }-u_{1}C)}^2{\Lambda }_v\left(p{\theta }_1{\kappa }_2+{\kappa }_2(1-p)({\theta }_1+{\omega }_1+{\mu }_h)\right)}{2{\mathcal{R}}_0{\Lambda }_m(\eta +{\mu }_h){({\mu }_h+{\alpha }_1{u}_1+{\alpha }_2{u}_2)}^2({\theta }_1+{\omega }_1+{\mu }_h){({\tau }_{10}+{\delta }_1+{\mu }_h+\epsilon u_3)}^2}}<0,\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{R}}_0}{\partial u_4}}& =-{\displaystyle \frac{c_2{c}_4{\beta }_{vr}{\beta }_{rv}{\mu }_r{\Lambda }_v}{2{\mathcal{R}}_0{\Lambda }_r{({\mu }_v+{\alpha }_1{u}_1+{\alpha }_2{u}_2)}^2{({\mu }_r+{\delta }_r+\gamma +u_4)}^2}}<0,\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill Q_1& ={\displaystyle \frac{\eta {\mu }_h{\beta }_{vh}{\beta }_{hv}{\Lambda }_v}{{\Lambda }_m(\eta +{\mu }_h)}}\left({\displaystyle \frac{p\zeta {\kappa }_1}{{\theta }_1+{\omega }_1+{\mu }_h}}+{\displaystyle \frac{p{\theta }_1{\kappa }_2+{\kappa }_2(1-p)({\theta }_1+{\omega }_1+{\mu }_h)}{({\theta }_1+{\omega }_1+{\mu }_h)({\tau }_{10}+{\delta }_1+{\mu }_h+\epsilon u_3)}}\right),\hfill \\ \hfill Q_2& ={\displaystyle \frac{\eta {\mu }_h{\beta }_{vh}{\beta }_{hv}{\Lambda }_v}{{\Lambda }_f(\eta +{\mu }_h)}}\left({\displaystyle \frac{q\zeta }{{\theta }_2+{\omega }_2+{\mu }_h}}+{\displaystyle \frac{q{\theta }_2+(1-q)({\theta }_2+{\omega }_2+{\mu }_h)}{({\theta }_2+{\omega }_2+{\mu }_h)({\tau }_{20}+{\delta }_2+{\mu }_h+\epsilon {\varphi }_2{u}_3)}}\right).\hfill \end{array}$$

The rate of change of ${\mathcal{R}}_0$ with respect to each controls ($u_1$,$u_2$, $u_3$ and $u_4$) are negative. This shows that increasing the values of the controls decreases the value of the basic reproduction number (${\mathcal{R}}_0$). Therefore, employing ITNs for male and female human populations, medical treatment for infected male and infected female human populations, vector control, and removing infected animals will decrease the average number of new infections and, hence, the prevalence of the disease. We can observe from Figure 2 that the use of additional effective medical treatment decreases the disease’s prevalence.

2.2. Optimal Control

2.2.1. Definition of Cost Function

We want to find the controls that minimize the total number of infected human populations and the cost of controls. This implies that finding optimal values of $u_1^{*}$, $u_2^{*}$, $u_3^{*}$ and $u_4^{*}$ that minimize the cost function $J(u_1,u_2,u_3,u_4)$ where,

$$\begin{array}{cc}\hfill J(u_1,u_2,u_3,u_4)=& {\int }_0^T[B_1{I}_m\left(t\right)+B_2{I}_f\left(t\right)+B_3{A}_m\left(t\right)+B_4{A}_f\left(t\right)+A_1{u}_3{I}_m\hfill \\ \hfill & +A_2{\varphi }_2{u}_3{I}_f+{\displaystyle \frac{1}{2}}(D_1{u}_1^2+D_2{u}_2^2+D_3{u}_3^2+D_4{u}_4^2)]\mathrm{d}t.\hfill \end{array}$$

(12)

subject to the differential equations (6), where T is the final time, the constants $B_1$, $B_2$, $A_1$, $A_2$, $D_1$, $D_2$, $D_3$ and $D_4$ are values that will balance the units of measurement and also may indicate the importance of one type of intervention over the other. $B_1{I}_m$, $B_2{I}_f$ represent the cost on the population of actively infected male and female individuals as well as $B_3{A}_m$ and $B_4{A}_f$ represent the cost on the population of asymptomatically infected male and female individuals. Moreover, the term $A_1{I}_m{u}_3$ and $A_2{\varphi }_2{I}_f{u}_3$ represent the individual treatment costs of male and female human host that receive treatment respectively, ${\displaystyle \frac{D_1}{2}}{u}_1^2$ represent the cost of implementing the control on ITNs, ${\displaystyle \frac{D_2}{2}}{u}_2^2$ represents the cost of vector control, ${\displaystyle \frac{D_3}{2}}{u}_3^2$ represent community level treatment cost related to the preparation and facility set-up for treatment and ${\displaystyle \frac{D_4}{2}}{u}_4^2$ represent the cost of removing infected reservoir animals. While the quadratic form of the cost function is mathematically advantageous, it also has practical biological and economic significance. The quadratic dependence on control efforts reflects the fact that as the intensity of interventions increases, costs tend to rise disproportionately. Several factors contribute to this phenomenon. For example, initial interventions, such as moderate dog culling or bed net distribution, might be relatively affordable. However, achieving broader coverage often requires significantly more resources, leading to higher logistical and operational expenses. Additionally, while culling a small number of infected reservoir animals may be manageable, scaling up this effort could encounter increased community resistance, further elevating costs. The quadratic formulation of control costs is widely used in objective functions of optimal control problems, offering a familiar and well-understood approach within the field. We seek to find optimal controls $u_1^{*}$, $u_2^{*}$, $u_3^{*}$ and $u_4^{*}$ such that

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

(13)

where the admissible set

$$\begin{array}{cc}\hfill \mathcal{U}=& \{(u_1,u_2,u_3,u_4)\in {\mathbb{R}}^4:u_1\left(t\right),u_2\left(t\right),u_3\left(t\right),u_4\left(t\right)\mathrm{are} \mathrm{Lebesgue} \mathrm{integrable},\hfill \\ \hfill & u_1\left(t\right)\in [0,u_{1\max }],u_2\left(t\right)\in [0,u_{2\max }],u_3\left(t\right)\in [0,1],u_4\left(t\right)\in [0,u_{4\max }].\}\hfill \end{array}$$

Considering the objective function in (12), we have the optimal control problem:

$$\begin{array}{cc}\hfill \underset{u_1,u_2,u_3,u_4}{\min }J=& {\int }_0^T[B_1{I}_m\left(t\right)+B_2{I}_f\left(t\right)+B_3{A}_m\left(t\right)+B_4{A}_f\left(t\right)+A_1{u}_3{I}_m\hfill \\ \hfill & +A_2{\varphi }_2{u}_3{I}_f++{\displaystyle \frac{1}{2}}(D_1{u}_1^2+D_2{u}_2^2+D_3{u}_3^2+D_4{u}_4^2)]\mathrm{d}t\hfill \end{array}$$

(14)

subject to the system in Eq. (6) and bounded controls.

2.2.2. Existence and Characterization of Optimal Control Solution

Examining the conditions that can assure the existence of a solution to our optimum control issue will be the first step. The following theorem shows the existence of the optimal control solution to the above optimal control problem [36,41].

Theorem 3.

There exists an optimal control ${\mathbf{u}}^{*}$, and corresponding solution vector ${\mathbf{x}}^{*}$ to the state initial value problem (6) that minimizes the objective function $J\left(\mathbf{u}\right)$ of (12) over the set of admissible controls $\mathcal{U}$.

Proof. 

The non-trivial requirements on the set of admissible controls $\mathcal{U}$ and on the set of end conditions are verified from Fleming & Rishel’s theorem [36].

• The set of all solutions to the system (6) with corresponding control functions in $\mathcal{U}$ is non-empty.
• The state system can be written as a linear function of the control variables with coefficients dependent on time and the state variables.
• The integrand L in (12) from objective functional with $L(\mathbf{x},\mathbf{u},t)=B_1{I}_m\left(t\right)+B_2{I}_f\left(t\right)+B_3{A}_m\left(t\right)+B_4{A}_f\left(t\right)+A_1{u}_3{I}_m\left(t\right)+A_2{\varphi }_2{u}_3{I}_f\left(t\right)+{\displaystyle \frac{1}{2}}\left(D_1{u}_1^2+D_2{u}_2^2+D_3{u}_3^2+D_4{u}_4^2\right)$ is convex on $\mathcal{U}$ and additionally satisfies $L(\mathbf{x},\mathbf{u},t)\ge {\delta }_1{\left|(u_1,u_2,u_3,u_4)\right|}^{\beta }-{\delta }_2$ where ${\delta }_1>0$ and $\beta >1$.
  In order to establish condition A, we refer to Picard-Lindelöf’s theorem from [37]. If the solutions to the state equations are bounded, and the state equations are continuous and Lipschitz with respect to the state variables, then for every admissible control $\mathcal{U}$, there exists a unique solution to the system. The total male human population is bounded between $N_{0m}$ and ${\displaystyle \frac{{\Lambda }_m}{{\mu }_h}}$, while the female human, reservoir and sand fly populations are bounded above ${\displaystyle \frac{{\Lambda }_f}{{\mu }_h}}$, ${\displaystyle \frac{{\Lambda }_r}{{\mu }_r}}$ and ${\displaystyle \frac{{\Lambda }_v}{{\mu }_v}}$ respectively, with lower bounds $N_{0f}$, $N_{0r}$ and $N_{0v}$. With these bounds, the state system is continuous, bounded, and Lipschitz with respect to the state variables, proving condition A holds. Condition B is verified by observing the linear dependence of the state equations on controls $u_1,u_2$, $u_3$ and $u_4$. Finally, to verify condition C by definition from [39,40] any constant, linear and quadratic functions are convex. So $u_1^2$, $u_2^2$,$u_3^2$ and $u_4^2$ are convex on $\mathcal{U}$. Since linear combination of convex functions are also convex, the integrand $L(\mathbf{x},\mathbf{u},t)$ is convex on $\mathcal{U}$.
  To prove the bound on the L, we adopted the method in [38]. We note that by the definition of $\mathcal{U}$, we have

  $$\begin{array}{cc}\hfill & D_4{u}_4^2\le D_4\mathrm{since}{u}_4\in [0,1]\hfill \\ \hfill & {\displaystyle \frac{D_4}{2}}{u}_4^2\le {\displaystyle \frac{D_4}{2}},\mathrm{implies}\mathrm{that}{\displaystyle \frac{D_4}{2}}{u}_4^2-{\displaystyle \frac{D_4}{2}}\le 0\hfill \\ \hfill & L(\mathbf{x},\mathbf{u},t)=B_1{I}_m\left(t\right)+B_2{I}_f\left(t\right)+B_3{A}_m\left(t\right)+B_4{A}_f\left(t\right)+A_1{u}_3{I}_m\hfill \\ \hfill & +A_2{\varphi }_2{u}_3{I}_f+{\displaystyle \frac{1}{2}}\left(D_1{u}_1^2+D_2{u}_2^2+D_3{u}_3^2+D_4{u}_4^2\right)\hfill \\ \hfill & \ge {\displaystyle \frac{1}{2}}\left(D_1{u}_1^2+D_2{u}_2^2+D_3{u}_3^2+D_4{u}_4^2\right)-{\displaystyle \frac{D_4}{2}}\hfill \\ \hfill & \ge \min {\displaystyle \frac{1}{2}}\left(D_1,D_2,D_3,D_4\right)\left(u_1^2+u_2^2+u_3^2+u_4^2\right)-{\displaystyle \frac{D_4}{2}}\hfill \\ \hfill & \ge \min {\displaystyle \frac{1}{2}}\left(D_1,D_2,D_3,D_4\right)|\left(u_1,u_2,u_3,u_4\right){|}^2-{\displaystyle \frac{D_4}{2}}\hfill \\ \hfill \therefore & L(\mathbf{x},\mathbf{u},t)\ge {\delta }_1|\left(u_1,u_2,u_3,u_4\right){|}^{\beta }-{\delta }_2\hfill \\ & \mathrm{where}{\delta }_1=\min {\displaystyle \frac{1}{2}}\left(D_1,D_2,D_3,D_4\right),{\delta }_2={\displaystyle \frac{D_4}{2}}\&\beta =2.\hfill \end{array}$$

   □

2.2.3. Characterization of Optimal Control Solution

To formulate the necessary conditions for optimality, we need to define the Hamiltonian function of the optimal control problem (14). The Hamiltonian equation with marginal cost function, state variables and adjoint variables is given by

$$\begin{array}{cc}\hfill & H(\mathbf{x},\mathbf{u},\mathbf{h},t)=B_1{I}_m\left(t\right)+B_2{I}_f\left(t\right)+B_3{A}_m\left(t\right)+B_4{A}_f\left(t\right)+A_1{u}_3{I}_m+A_2{\varphi }_2{u}_3{I}_f\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}(D_1{u}_1^2+D_2{u}_2^2+D_3{u}_3^2+D_4{u}_4^2)+\sum _{i=1}^{15}{h}_{i}F(\mathbf{x},\mathbf{u},\mathbf{h},t)\hfill \end{array}$$

By using Pontryagin’s Maximum Principle [36], we derive necessary conditions for the optimal controls and corresponding states. That is, it is now possible to determine the optimal control variables, $u_1^{*},u_2^{*},u_3^{*}$ and $u_4^{*}$, from the necessary conditions. Now, the necessary conditions are given as follows.

• Optimality Conditions: The minimization of the Hamiltonian H with respect to the control variables $u_1,u_2,u_3,u_4$, is the first condition that we will examine from the Pontryagin’s Maximum principle. Since the cost function is convex, if the optimal control occurs in the interior region we must have ${\displaystyle \frac{\partial H}{\partial u_i}}=0$ for $i=1,2,3,4$. Therefore

  (i)

  for the control $u_1$ we must have,

  $$\begin{array}{cc}\hfill {\displaystyle \frac{\partial H}{\partial u_1}}& =0,\Rightarrow {\overline{u}}_1=C{\beta }_{hv}{\displaystyle \frac{S_v}{D_1}}({\displaystyle \frac{{\varphi }_1}{N_f}}(h_{12}-h_{11})(\zeta A_f+I_f)+{\displaystyle \frac{1}{N_m}}(h_{12}-h_{11})({\kappa }_1\zeta A_m\hfill \\ \hfill & +{\kappa }_2{I}_m\left)\right)+C{\beta }_{vh}{\displaystyle \frac{I_v}{D_1}}\left({\displaystyle \frac{S_m}{N_m}}(h_2-h_1)+{\varphi }_1{\displaystyle \frac{S_f}{N_f}}(h_7-h_6)\right)+{\displaystyle \frac{{\alpha }_1}{D_1}}(S_v{h}_{11}+I_v{h}_{12}).\hfill \end{array}$$

  (ii)

  for the control $u_2$ we must have,

  $$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial H}{\partial u_2}}=0,\Rightarrow {\overline{u}}_2={\displaystyle \frac{1}{D_2}}[{\Lambda }_v{h}_{11}+{\alpha }_2(S_v{h}_{11}+I_v{h}_{12})].\hfill \end{array}$$

  (iii)

  for the control $u_3$ we must have,

  $$\begin{array}{cc}\hfill {\displaystyle \frac{\partial H}{\partial u_3}}=0,\Rightarrow {\overline{u}}_3& ={\displaystyle \frac{1}{D_3}}[-A_1{I}_m-A_2{\varphi }_2{I}_f+\epsilon (h_4-h_5)I_m+\epsilon {\varphi }_2(h_9-h_{10})I_f]\hfill \end{array}$$

  (iv)

  the control $u_4$ we must have,

  $$\begin{array}{c}\hfill {\displaystyle \frac{\partial H}{\partial u_2}}=D_4{u}_4-I_r{h}_{14}\Rightarrow {\overline{u}}_4={\displaystyle \frac{I_r}{D_4}}{h}_{14}.\end{array}$$

  And therefore, the optimal controls on the given bounded intervals are given by

  $$\begin{array}{cc}\hfill u_1^{*}& =\min \{u_{1\max },\max \{0,{\overline{u}}_1\}\},u_2^{*}=\min \{{\displaystyle \frac{1}{{\alpha }_2}}(1-{\mu }_v-\alpha 1u_{1\max }),\max \{0,{\overline{u}}_2\}\},\hfill \\ \hfill u_3^{*}& =\min \{1,\max \{0,{\overline{u}}_3\}\},u_4^{*}=\min \{1-{\mu }_r-{\delta }_r,\max \{0,{\overline{u}}_4\}\}.\hfill \end{array}$$
• The adjoint (co-state) equations: According to the second condition of Pontryagin’s Maximum Principle, we need to have the optimal controls ${\displaystyle \frac{\partial H}{\partial x_i}}=-{\displaystyle \frac{d{h}_i}{dt}}$ for each $i=1,2,\cdots ,15$. As a result, we must compute and resolve the system,

  $$\begin{array}{cc}\hfill & h_1^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial S_m}}{h}_2^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial E_m}}{h}_3^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial A_m}}{h}_4^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial I_m}}{h}_5^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial R_m}}\hfill \\ \hfill & h_6^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial S_f}}{h}_7^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial E_f}}{h}_8^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial A_f}}{h}_9^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial I_f}}{h}_{10}^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial R_f}}\hfill \\ \hfill & h_{11}^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial S_v}}{h}_{12}^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial I_v}}{h}_{13}^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial S_r}}{h}_{14}^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial I_r}}{h}_{15}^{\prime }\left(t\right)=-{\displaystyle \frac{\partial H}{\partial R_r}}.\hfill \end{array}$$

  (15)
• The transversality conditions

  $$h_1\left(T\right)=h_2\left(T\right)=h_3\left(T\right)=\cdots =h_{15}\left(T\right)=0.$$

  (16)

Since the model functions exhibit convexity with respect to the control variables, and considering the predetermined boundedness of the state and adjoint (or co-state) functions, the resulting optimal solution is unique for short duration of time T [42,43,44].

2.3. Numerical Methods

We numerically calculated optimal control strategies based on the iterative method used in [45] in the interval [0,156] weeks or for three years. Given initial estimates for the controls and initial conditions for the states, the state differential equations (6) are then solved forward in time using the $4^{\mathrm{th}}$ order Runge-Kutta method. According to the results of state values and the given final value in Eq. (16), we solved adjoint values from the adjoint equations backward in time, using the fourth-order Runge-Kutta method. Both updates of state values and adjoint values were utilized to calculate optimal control strategies. The same process is repeated until the current state, adjoint, and control values converge sufficiently [47]. All our numerical simulations are performed using the MATLAB software. The summary of the algorithm based on [46] is written as follows:

Algorithm 1: Forward-Backward-Sweep Algorithm

It is challenging to pinpoint the precise cost of pharmaceutical therapy for VL due to a variety of factors, including geographical variations, shifts in medication prices, the nature of the health-care system, and indirect expenses. The majority of the total cost per patient is frequently attributed to prescription prices. However, other elements, including hospitalization, diagnostics, and supportive care, also play a role. Depending on the patient’s location and the availability of WHO drug subsidies, the most effective VL medication, amphotericin B, could cost $\$200$ per treatment in many endemic counties [49]. Therefore, in our simulation, $A_1=\$200$ and $A_2=\$200$ are used. The cost of culling an infected dog and its disposal varies based on the dog’s weight, the type of service, and the location of the facility. Considering these, we estimate the coefficient for the cost of dog culling to be $\$150$. Hence $D_4=150$ is used. The cost of Insecticide-Treated Nets (ITNs) for controlling sand flies, particularly in VL endemic regions, is highly variable based on delivery systems, geography, and program scale. Since there is no pre-estimated cost for ITNs and environmental management as well as indoor residual spraying cost to control sand flies, we have assumed the total direct and indirect cost of those interventions for our numerical simulation to be $B_1=18$, $B_2=15$, $B_3=12$, $B_4=10$, $D_1=300$, $D_2=200$ and $D_3=400$.

For numerical simulations, the initial population values are taken from the fitted values at the time of last data recorded in [10]. These are $S_m\left(0\right)=920,091$, $E_m\left(0\right)=594$, $A_m\left(0\right)=136$, $I_m\left(0\right)=1,129$, $R_m\left(0\right)=7,020$, $S_f\left(0\right)=380,814$, $E_f\left(0\right)=69$, $A_f\left(0\right)=4$, $I_f\left(0\right)=175$, $R_f\left(0\right)=895$, $S_v\left(0\right)=1,254$, $I_v\left(0\right)=130$, $S_r\left(0\right)=115,290$, $I_r=10,463$ and $R_r=919$. The other model parameter values are as indicated in Table 2.

3. Results and Discussion

3.1. Various Interventions/Control Strategies

We have introduced different possible combinations of strategies in our model and numerically compared their effects on the total number of infected populations. The time profile of infected classes when no intervention is applied are indicated in Figure 3.

We can observe from Figure 3 that, in the absence of interventions, the disease continues to spread and remains persistent within the community. This persistence highlights the robustness of the transmission dynamics and the difficulty in achieving natural eradication. Without targeted measures, the disease remains endemic and may lead to long-term health and socio-economic impacts. Therefore, it becomes essential to implement a range of intervention strategies, such as insecticide-treated bed nets, vector control through indoor residual spraying, treatment programs, and screening and culling of infected reservoirs, to reduce disease prevalence significantly. In the following sections, we will see the impact of different interventions in reducing the number of infected populations.

Despite the ethical concerns, we simulated the total prevalence and corresponding costs to assess the impact of applying control interventions at varying rates for both male and female humans. When the proportion of treatment and ITNs used for the female human population as compared to their male counterpart is set at $50\%$, $75\%$, or $100\%$, there is no significant difference in the total prevalence or intervention costs (see Figure 4). However, when the proportion is reduced to 10%, a noticeable difference in prevalence is observed. Therefore, to optimize resource use and reduce prevalence, it is advisable to take this proportion into consideration.

We now analyse the effect of various possible combinations of intervention mechanisms, in reducing the burden of the disease. In each of these combinations, we assumed that only $50\%$ of the controls $u_1$ and $u_3$ are implemented to the female humans. This is based on the information from Figure 4.

3.1.1. Strategy A (All Controls i.e., the Combination of the Use of Treated Bed-Nets, Vector Control, Treatment of Infected Individuals and Culling Infected Reservoir Animals

In this approach, we applied all the four control variables, $u_1$, $u_2$, $u_3$ and $u_4$, to optimize the objective function J in Eq.(12). When comparing strategy A to the scenario without any control measures, the number of infected male humans ($I_m$), the infected sand fly population ($I_v$), and the number of infected reservoirs ($I_r$) show an exponential decline (refer to Figure 5 (a), (c) and (d) respectively). Whereas the number of infected female human population increases for the first about 16 weeks and goes down and becomes optimum at week 130 (Figure 5 (b)). When analyzing the control profile (Figure 5 (e)), it is clear that the first three controls ($u_1$, $u_2$ and $u_3$) remain at their maximum levels ($100\%$) for approximately 126 weeks. In contrast, the fourth control ($u_4$) remains at the maximum level for the first five weeks and drops from its peak by $60\%$ at the fifth week and gradually decreases to zero by the end of the intervention period.

3.1.2. Strategy B (the Combination of the Use of Medical Treatment for Infected Individuals and Animal Culling Without Using Insecticide-Treated Bed Nets and Vector Controls ($u_1=0$ and $u_2=0$)

In this strategy, we applied the combination of medical treatment($u_3$) and culling infected animals ($u_4$) without the use of insecticide-treated bed-net ($u_1$) and vector controls ($u_2$). When comparing Strategy B to an uncontrolled scenario, it can be observed that the population of infected male humans increases during the first 8 weeks, then drops by 62.5% over the following about 4 weeks, then decreases slowly until 133 weeks. Then the number of infected male increases exponentially for the rest of the intervention period (Figure 6 (a)). Similarly, the population of infected female humans initially rises for the first 24 weeks, then decreases until weeks 133 and goes up like the number of males (Figure 6 (b)). This pattern is due to the therapy starting at a low level for the first 8 weeks before reaching its maximum (see Figure 6 (e)). As we can see in Figure 6 (c) and (d), although this strategy decreases the number of reservoir and infected sand fly, it can’t decrease these variables to their lowest possible level like that of Strategy A. When the therapy intensifies, the number of infected individuals starts to decrease, but as the treatment level reduces, the infected populations of both male and female humans begin to rise again. Maintaining animal culling at its maximum level for the first 133 weeks ensures that the number of infected reservoirs remains at an optimal low level throughout this period.

3.1.3. Strategy C (the Combination of the Use of Insecticide-Treated Bed Nets and Vector Control Without the Use of Medical Treatment and Animal Culling ($u_3=0$ and $u_4=0$)

In strategy C, we employed a combination of insecticide-treated bed nets and vector control, while other interventions like medical treatment and culling of infected animals were set to zero. This strategy results in a significant reduction by $97.31\%$ of the number of infected males (Figure 7 (a)) from the baseline by the end of the intervention period. Although it also reduces the number of infected females significantly (by $88.85\%$) as compared to the scenario without any control measures, the decline is lower than that of males (Figure 7 (b)). This could be due to the $50\%$ use of the control $u_1$ for the female human population. With the use of this Strategy, one can observe that the number of infected vectors and infected reservoirs also decrease by $96.14\%$ and $93.68\%$ respectively compared to the baseline at the end of the intervention period 7 (c) and (d)). Figure 7 (e) indicates that the control profile for insecticide-treated bed nets and vector control should continue at their full intensity for the whole planning period.

3.1.4. Strategy D (the Combination of Use of Vector Control and Animal Culling Without Insecticide Treated Bed-Nets and Medical Treatment ($u_1=0$ and $u_3=0$)

Here, we employed the combination of the use of vector control and animal culling without insecticide-treated bed nets and medical treatment. When comparing the effect of Strategy D to the uncontrolled scenario, there is an exponential decline in the number of infected male and female humans, as well as sand fly populations (see Figure 8 (a), (b) and (c)). The number of infected reservoirs increases for the first three weeks and decreases until week 14, then rises (Figure 7 (d)). To achieve this, vector control is maintained at its maximum levels for the whole planning period, and animal culling is also at its maximum level for the first 25 weeks before going to the minimum, as shown in Figure 7 (e)).

3.1.5. Strategy E (the Combination of the Use of Insecticide-Treated Bed-Nets, Medical Treatment, and Animal Culling Without the Use of Vector Control ($u_2=0$)

In this strategy, the combination of the use of ITNs, vector control, and animal culling is employed without the use of medical treatment.

We can observe from Figures 9 (a), (b), (c), and (d)) that this strategy significantly reduces the number of all infected classes like that of Strategy A. However, it takes a long time for the number of infected sand flies to reach the minimum level as compared to that of strategy A due to the absence of a vector control strategy in this approach. The number of infected female human populations also does not reach the minimum possible level in the intervention period as to that of Strategy A. The profile of control graphs is also similar to Strategy A.

3.1.6. Strategy F (the Combination of the Use of Insecticide-Treated Bed Nets, Vector Control and Medical Treatment Without the Use of Animal Culling ($u_4=0$)

In this approach, we employed the combination of the three controls $u_1$, $u_2$, and $u_3$ without using infected animal culling. In this strategy, the number of infected reservoirs increases exponentially for the first 7.8 weeks and turns to decrease. Whereas the number of infected male, female, and sand fly populations decrease significantly (Figures 10 (a), (b), (c) and (d)).

The intervention Strategies have been implemented for 156 weeks or three years following the VL data recorded in the study by Awoke et al. [10].

When all control measures are implemented, the total disease prevalence reaches its optimal level. As shown in Figure 11, Strategy F and Strategy E also prove to be effective in reducing the total prevalence, ranking just below Strategy A in terms of the level of total prevalence and their corresponding total cost. On the other hand, Strategy D does not seem to help in controlling the disease, as both the disease prevalence and associated costs increase during the planning period unlike the others. The cost-effectiveness of these strategies will be analyzed in the subsection below.

3.1.7. Cost-Effectiveness Analysis of the Proposed Strategies

Cost-effectiveness analysis is a method used to compare the relative costs and outcomes of different interventions or strategies. In the context of optimal control strategies, cost-effectiveness analysis evaluates whether the benefits achieved by implementing a particular control strategy outweigh its costs. In order to understand the cost-effectiveness of each of the optimal combinations, the incremental cost-effectiveness ratio (ICER) is calculated for each strategy. The incremental cost-effectiveness ratio (ICER) is the additional cost per health outcome [8]. It is the ratio of the difference in total cost between one strategy and the next best alternative to the difference in total number of averted infections through each strategy [46,48]. The number of infections averted is the difference between the total number of infected individuals without control and the total count of infected individuals with control obtained during the control period. To implement the ICER, we simulate the model corresponding to each of the six intervention strategies and apply the formula in Eq.(17). Using these simulation results, we rank the control strategies in increasing order of effectiveness based on infection averted as indicated in Table 3. Accurately estimating the costs associated with interventions that prevent infections and connecting those costs to the Quality Life Adjusted Years (QALYs) gained is crucial for ICER calculation. However, associating the computation of ICER with QALYs is complicated due to the indirect costs, opportunity costs, and variations in health-care systems, as well as other factors. As a result, we did not express ICER in terms of QALYs in this study.

$$\begin{array}{c}\hfill \mathrm{I}\mathrm{C}\mathrm{E}\mathrm{R}={\displaystyle \frac{\mathrm{D}\mathrm{i}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{t}\mathrm{o}\mathrm{f}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{g}\mathrm{i}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{j}}{\mathrm{D}\mathrm{i}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{g}\mathrm{i}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{j}}}\end{array}$$

(17)

After arranging the Strategies in ascending order based on the total infections averted, the Incremental Cost-Effectiveness Ratio (ICER) is calculated using the method outlined by [9]. The ICER values are then compared pairwise to identify the most cost-effective strategy.

$$\begin{array}{cc}\hfill \mathrm{I}\mathrm{C}\mathrm{E}\mathrm{R}D& ={\displaystyle \frac{1.9090\times {10}^9}{33,743}}=56,574.69\hfill \\ \hfill \mathrm{I}\mathrm{C}\mathrm{E}\mathrm{R}B& ={\displaystyle \frac{1.0760\times {10}^9-1.9090\times {10}^9}{123,876.3-33,743}}=-924.18\hfill \\ \hfill \mathrm{I}\mathrm{C}\mathrm{E}\mathrm{R}C& ={\displaystyle \frac{1.4173\times {10}^8-1.0760\times {10}^9}{153,269.1-123,876.3}}=-68,393.28\hfill \\ \hfill \mathrm{I}\mathrm{C}\mathrm{E}\mathrm{R}E& ={\displaystyle \frac{5.0742\times {10}^7-1.4173\times {10}^8}{163,154.8-153,269.1}}=-9,204\hfill \\ \hfill \mathrm{I}\mathrm{C}\mathrm{E}\mathrm{R}F& ={\displaystyle \frac{3.3250\times {10}^7-5.0742\times {10}^7}{163,299.81-163,154.8}}=-120,626.16\hfill \\ \hfill \mathrm{I}\mathrm{C}\mathrm{E}\mathrm{R}E& ={\displaystyle \frac{3.2144\times {10}^7-3.3250\times {10}^7}{163,304.15-163,299.81}}=-679.65.\hfill \end{array}$$

The pairwise comparisons of the calculated ICER for the different scenarios and corresponding total cost are given in Table 4.

From ICER of D and ICER of B, it is evident that Strategy B’s ICER is lower than that of Strategy D. This demonstrates that Method D is less successful in saving lives and more costly. As a result, Strategy B saves more lives than Strategy D. Strategy D was omitted as a competing strategy, and we calculated the ICER for the remaining strategies, B, C, E, F and A, following the same procedure as in the above, and the result is as shown in Table 5.

$$\begin{array}{cc}\hfill \mathrm{ICER}B& ={\displaystyle \frac{1.0760\times {10}^9}{123,876.3}}=8,686.08\hfill \\ \hfill \mathrm{ICER}C& ={\displaystyle \frac{1.4173\times {10}^8-1.0760\times {10}^9}{153,269.1-123,876.3}}=-31,785.68\hfill \\ \hfill \mathrm{ICER}E& ={\displaystyle \frac{5.0742\times {10}^7-1.4173\times {10}^8}{163,154.8-153,269.1}}=-9,204\hfill \\ \hfill \mathrm{ICER}F& ={\displaystyle \frac{3.3250\times {10}^7-5.0742\times {10}^7}{163,299.81-163,154.8}}=-120,626.16\hfill \\ \hfill \mathrm{ICER}A& ={\displaystyle \frac{3.2144\times {10}^7-3.3250\times {10}^7}{163,340.15-163,299.81}}=-27,416.96.\hfill \end{array}$$

As we can see from Table 5, when we compare the ICERs of Strategy B and Strategy C, Strategy C has a lower ICER than Strategy B. This shows that Strategy B is strongly dominated by the lower ICER over Strategy C. In other words, compared to Strategy C, Strategy B is less successful and costs more. As a result, Strategy B is excluded in the list of options. Following the same procedures, we have computed the ICERs values as indicated in Table 6.

From Table 6, by comparing the ICERs of Strategy C and Strategy E, we can observe that Strategy E has a lower ICER and is hence cost effective than Strategy C.

By repeating this process, we can identify the next most cost-effective strategy. From this analysis, Strategy E is found to be the next in line, followed by Strategy F. A final comparison of Strategies F and A shows that Strategy F is strongly dominant, which means it is both more expensive and less effective than Strategy A. Therefore, Strategy A, which is the combination of all four controls, i.e., insecticide-treated bed nets, vector control, treatment of infected human populations, and infected animal culling, emerges as the most cost-effective disease control strategy, as it has the lowest ICER.

The cost-effectiveness analysis results of this study align with those of a study conducted by Agusto et al. [9] even if our control Strategy proposes only $50\%$ of the controls in insecticide-treated bed-nets and treatment to be provided to the female humans as compared to the one to be given to the male counterpart. Their findings indicate that the most cost-effective intervention strategy for controlling the transmission of malaria-VL co-infection involves using personal protection measures, insecticide-treated bed nets, culling infected animal reservoirs, and applying insecticides.

According to the calculated ICER, Strategy F, which involves all control strategies except animal culling, is the second cost-effective intervention. This finding is consistent with the results of Zhao et al. [19], who state that animal culling is an ineffective VL control strategy. The third most cost-effective strategy, according to our comparison result, includes using insecticide-treated be nets, medical treatment, and animal culling without using vector control or Strategy E.

4. Conclusions

In this paper, a mathematical model for the dynamics of sex-structured visceral leishmaniasis transmission model is presented and analyzed. We have added four time-dependent control variables to establish optimal disease control strategies. These control variables include treating infected male and female human populations, vector control, using bed-nets treated with insecticide and culling of infected animals. Positivity and boundedness of the solution, as well as the basic reproduction number with constant control, are computed. We derived and analyzed the necessary conditions for the optimal control problem by using Pontryagin’s maximum principle. The numerical simulations show that the overall prevalence of the disease will decrease over time with the inclusion of each of these control techniques into the model. Cost-effectiveness analysis is employed to evaluate the economic implications of the various combinations of the four control measures for VL. To investigate the cost-effectiveness of every possible strategies, we have computed the Incremental Cost-Effectiveness Ratio (ICER). We did a pairwise comparison of the ICER of each strategy and identified the strategy with a smaller ICER, which averts more infections in humans and utilizes the smallest possible cost of intervention.

The result shows a comprehensive evaluation of various intervention strategies aimed at reducing disease transmission and prevalence. The findings show that the most effective and cost-efficient approach is Strategy A, which combines all the four key controls: insecticide-treated bed nets, vector control, medical treatment of infected individuals, and culling of infected animal reservoirs. This strategy leads to a significant and sustained reduction in the number of infected male and female humans, vectors (e.g., sand flies), and animal reservoirs. The intervention’s success is attributed to the simultaneous targeting of multiple transmission pathways, ensuring that disease prevalence decreases rapidly and remains low throughout the intervention period.

Cost-effectiveness analysis, using the Incremental Cost-Effectiveness Ratio (ICER), further validates Strategy A as the optimal solution, providing the greatest reduction in disease prevalence for the lowest cost. Strategy E, which excludes vector control, and Strategy F, which omits animal culling, are also effective but rank lower in cost-effectiveness. These strategies reduce disease prevalence but do not achieve the same level of control as Strategy A. The study highlights that strategies like B (medical treatment and culling without bed nets or vector control) and D (vector control and culling without bed nets or medical treatment) are less efficient. Although they reduce some components of the infected population, they fail to maintain low infection levels over the intervention period and are more costly, making them suboptimal.

Our findings emphasize the importance of a multi-pronged approach for VL control. Strategies that focus on just one or two interventions, or that limit the application of controls to certain populations, are less successful and more expensive in the long term. The analysis was done by applying a 50% control strategy to the female humans, compared to their male counterparts. As can be seen from the evaluation, it can still effectively eliminate the burden of VL from the entire population within the planning period. Integrated intervention strategies, like Strategy A provides the most effective means of reducing both disease prevalence and long-term socio-economic impacts. Furthermore, the study underscores the need for a careful balance between cost and effectiveness, as shown by the ICER analysis, to optimize resource allocation in public health programs to control persistent infectious diseases.