Optimal control and cost-effectiveness analysis for leptospirosis epidemic

Abstract

This paper aims to apply an optimal control theory for the autonomous model of the leptospirosis epidemic to examine the effect of four time-dependent control measures on the model dynamics with cost-effectiveness. Pontryagin's Maximum Principle was used to derive the optimality system associated with the optimal control problem. Numerical simulations of the optimality system were performed for different control strategies and the results were presented graphically with and without controls. The optimality system was simulated using the Forward–Backward Sweep method in the Matlab programme. The numerical results revealed that the combination of all optimal control measures is the most effective strategy for minimizing the spread and impact of disease in the community. Furthermore, a cost-effectiveness analysis was performed to determine the most cost-effective strategy using the incremental cost-effectiveness ratio approach and we observed that the rodenticide control-only strategy is most effective to combat the spread of disease when available resources are limited.

Keywords:

• Leptospirosis
• leptospira
• optimal control
• numerical simulation
• cost-effectiveness

1. Introduction

Leptospirosis is a major zoonosis disease caused by bacterial pathogens, spirochete leptospira interrogans [11, 33, 36]. The disease is typically found in tropical and subtropical regions, especially in developing countries including South-East Asia countries and Sub-Sahara Africa but is increasingly being reported from many parts of the world [7, 9, 11]. Recently, the incidence of leptospirosis was estimated to be 1.03 million cases and 58,900 deaths annually worldwide [4, 10, 21]. Most of the outbreaks of the disease occur during the rainy season after happening of heavy rainfall and typhoons in the countries where it is endemic [34].

Rats (rodents) are the primary sources of pathogens worldwide, particularly in urban slum areas where the disease is endemic [10, 29, 33]. The high incidence of infections in the rat population is mostly to blame for the high concentration of leptospira bacteria shed in the slum surroundings [24]. The pathogen usually infects people and other animals by direct contact with the disease's source or indirectly contact with a contaminated environment, like the urine of an affected animal [11, 24, 33, 35]. Rarely does the disease spread from person to person [8]. Leptospirosis illness in humans can be seen in two stages. The acute (primary) stage is characterized by mild illness with nonspecific signs such as higher fever, headache, and conjunctival suffusion [11, 12]. While syndromes in the second (severe illness or immune) stage of leptospirosis include jaundice, kidney failure, haemorrhage (especially pulmonary), meningitis, cardiac arrhythmias, respiratory insufficiency, and hemodynamic collapse [8, 11].

Human vaccine against leptospirosis is not widely practiced and accessible only in a few developed countries to protect only against the serovar in the vaccine [11, 23]. However, leptospirosis infections can be minimized by providing personal protection and control intervention efforts. Personal protection measures include personal protective equipment such as wearing rubber boots, dressings to cover wounds or skin, goggles, and rubber gloves, especially for those whose, working environment exposes them to risk of infections and maintaining good personal hygiene. The growth of the concentration of spirochete leptospiries can be reduced by making environmental modifications like draining moist areas or improving the slum areas [8, 12]. Patients with early leptospirosis infection (mild illness) can be treated through antibiotic doxycycline, ampicillin or amoxicillin, azithromycin, or clarithromycin. While patients with severe illness can be treated through IV penicillin and ceftriaxone drugs [8, 23].

The study of the infectious diseases epidemiology requires the use of mathematical modelling tools because these tools can provide some insight into the transmission dynamics of disease and help in the identification of appropriate control measures [5]. Recently the burden of infectious diseases has led to looking for an optimal control technique to reduce the spread of infectious diseases in an infected population [16, 22]. Optimal control theory helps in evaluating the effectiveness of different preventive and control programmes as well as the cost-effectiveness of implementing the interventions.

A few researchers have studied a non-autonomous mathematical model of the leptospirosis epidemic by incorporating different preventive and control measures in their models. Authors in [25] examined a two strain compartmental model by incorporating time-dependent preventive and control measures to minimize the risk and spread of the disease in infected human and vector populations. They achieved the effectiveness as well as cost-effectiveness of the intervention strategies. While [24] proposed and achieved the effects of time-dependent control measures on rodent infection dynamics and concentration of leptospira in the environment using optimal control theory. They considered the two controls namely, rodenticide and resource reduction on rodent populations and the two permanent environment controls, namely the natural death rate of leptospira and carrying capacity controls. However, humans can acquire leptospirosis infection through multiple pathway in the community. Moreover, reducing rodents by using rodenticide will cause a high concentration of leptospira bacteria in the environment due to the high rate of shedding of leptospira from dead bodies of rodents. Consequently, preventive measures to reduce the concentration of the bacteria in the environment should be implemented. Likewise, the researchers in [28] also, used the optimal control theory on transmission dynamics of leptospirosis disease on infected human and vector populations to reduce the impact of the epidemic in the human population. Moreover, a few other researchers have been studying the deterministic autonomous model for transmission dynamics of leptospirosis epidemic using a compartmental approach in different forms [17, 19]. In recent times, authors in [11] developed and examined a deterministic mathematical model for transmission dynamics of leptospirosis described by nonlinear ordinary differential equations with eight state variables in human, rodent and bacterial populations. They demonstrated a detailed qualitative and quantitative analysis to show more insight into the spread of disease. They performed a sensitivity analysis to identify the model parameters that had an impact on the dynamics behaviour of their model. To corroborate their theoretical findings and the outcomes of their sensitivity analysis, they also carried out numerical simulations. Motivated by the above research, this paper aims to evaluate the effectiveness of multiple control measures and the cost-effectiveness of different optimal control interventions for the study in [11]. To achieve this, we incorporated four time-dependent preventive (or control) measures. These are: two measures on human population – appropriate personal protective equipment (PPE) to minimize the acquisition of the infection in humans and treatment control to reduce the number of infected individuals in the human population; one control measure on rodent population – rodenticide to reduce the total rodent population; and one preventive measure on pathogen – improvements in slum areas (or draining wet areas) to reduce the growth of the pathogenic spirochete leptospira bacteria in the environment.

The remaining sections of the paper are organized as follows. The general description of the autonomous model from the study [11] is presented in Section 2. Section 3 is devoted to the formulation of the optimal control problem of the leptospirosis epidemic model in [11], the proof of the standard results for the existence and characterization of the optimal control problem. Section 4 presents numerical simulations and cost-effectiveness analysis which includes numerical simulations of various strategic control strategies and their cost-effectiveness analysis. Concluding remarks are resented in Section 5.

2. The autonomous model

The developed model of the leptospirosis epidemic in [11] includes two host populations (namely, human and rodent), and a bacterial population. The human population at time t, represented by $N_h(t)$, is subdivided into four distinct compartments, given as: susceptible humans $S_h(t)$, latently infected humans $E_h(t)$, infectious humans $I_h(t)$, and humans recovered from the disease $R_h(t)$. Thus, the total human population is given by (1) $\begin{array}{r}{N}_h(t)=S_h(t)+E_h(t)+I_h(t)+R_h(t).\end{array}$(1) Whereas, the total rodent population at time t, denoted by $N_v(t)$ consists of the following distinct classes: susceptible rodents $S_v(t)$, infected rodents $I_v(t)$, and rodents recovered from the disease $R_v(t)$. Thus, $N_v(t)$ is given by (2) $\begin{array}{r}{N}_v(t)=S_v(t)+I_v(t)+R_v(t).\end{array}$(2) Assumed the concentration of the pathogen, leptospira interrogans in the environment at time t, is represented by $B_l(t)$. Therefore, the model includes eight distinct compartments. The susceptible humans and rodents are recruited, respectively, at rates Λ and Π. The leptospira interrogans in the environment $B_l(t)$, and rodents in the infected class $I_v(t)$ are assumed of transmitting the rodent-born leptospirosis to susceptible individuals at the transmission rate ${\beta }_2$ and ${\beta }_1$, respectively. It is also, assumed that the susceptible rodents may be infected through contact with individuals in the infectious class $I_h(t)$, at the transmission rate ${\beta }_3$. The progression rate of latently infected individuals to infectious class is denoted by θ. The rate of recovery of infectious individuals and infected rodents, respectively, are denoted as δ and σ. The parameters γ and ρ, respectively, denote the disease-waning immunity of humans and rodents. Here, the disease-induced death rate of humans is represented by α. The human, rodent, and bacteria mortality rates are denoted by $\mu ,{\mu }_v$, and ${\mu }_b$, respectively. The constant pathogen concentration in the environment is denoted by κ. Furthermore, the population $B_l$ increases its size in a contaminated environment from the release of bacteria by infected humans and rodents with the rate of ${\tau }_1$ and ${\tau }_2$, respectively. Therefore, the transmission dynamics autonomous model for the leptospirosis epidemic from [11] is given by the following set of nonlinear differential equations: (3) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-({\lambda }_h+\mu )S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& ={\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =\delta I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-{\mu }_b{B}_l,\end{array}\end{array}$(3) with the initial conditions: $S_h(0)=S_{h0}\ge 0,E_h(0)=E_{h0}\ge 0,I_h(0)=I_{h0}\ge 0,R_h(0)=R_{h0}\ge 0,S_v(0)=S_{v0}\ge 0,I_v(0)=I_{v0}\ge 0,R_v(0)=R_{v0}\ge 0,B_l(0)=B_{l0}\ge 0,\mathrm{where},{\lambda }_h=\frac{{\beta }_2{B}_l}{\kappa +B_l}+{\beta }_1{I}_v$.

In [11], the researchers presented and rigorously analysed both qualitative and quantitative behaviours of the autonomous model (3(3) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-({\lambda }_h+\mu )S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& ={\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =\delta I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-{\mu }_b{B}_l,\end{array}\end{array}$(3) ) to describe the spread of the epidemic. The results for non-negativity and boundedness of the model solutions were established, and stability analysis of steady state was carried out. Based on the next-generation matrix approach, the basic reproduction number $(R_0)$ of the ‘formulated model’ was computed. $R_0$ from [11] is found as (4) $\begin{array}{r}{R}_0=\rho (G)=\frac{1}{2}(\frac{{\beta }_2\mathrm{\Lambda }\theta {\tau }_1}{\kappa \mu {\mu }_b{ϵ}_1{ϵ}_2}+\sqrt{{\left(\frac{{\beta }_2\mathrm{\Lambda }\theta {\tau }_1}{\kappa \mu {\mu }_b{ϵ}_1{ϵ}_2}\right)}^2+4\frac{{\beta }_3\mathrm{\Pi }\mathrm{\Lambda }\theta ({\beta }_2{\tau }_2+{\beta }_1\kappa {\mu }_b)}{\kappa \mu {\mu }_v{\mu }_b{ϵ}_1{ϵ}_2{ϵ}_4}}).\end{array}$(4) where, ${ϵ}_1=\theta +\mu ,{ϵ}_2=\alpha +\delta +\mu ,{ϵ}_4=\sigma +{\mu }_v.$

The autonomous system undergoes forward bifurcation at $R_0=1$ as confirmed in [11], this shows that the leptospirosis can be eliminated in the population as long as when $R_0<1$ using possible intervention strategies. For a model that experiences a forward bifurcation, the condition $R_0<1$ is a necessary and sufficient condition for the eradication of the disease [15]. The most influencing model parameters on $R_0$ were also determined using sensitivity analysis. Further, numerical simulations were carried out to show the stability behaviour of a unique endemic equilibrium and to assess the chaining effect of the most sensitive parameters on the model dynamics.

3. Optimal control analysis of the leptospirosis model

In the work by [11], the most sensitive parameters in their basic reproduction determined as: contact rate among infected rodents and susceptible individuals ${\beta }_1$, and transmission rate of rodent infection ${\beta }_3$, human recovery rate δ, and mortality rate of rodents ${\mu }_v$ among others. The effect of these parameters on the disease transmission dynamics were simulated numerically. Moreover, the basic reproduction number of the model was found as, $R_0\approx 2.8892>1$. Based on the findings from the sensitivity analysis in [11], this study extends the autonomous mathematical model (3(3) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-({\lambda }_h+\mu )S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& ={\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =\delta I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-{\mu }_b{B}_l,\end{array}\end{array}$(3) ) into a non-autonomous one to demonstrate the impact of four time-dependent control functions on the dynamics of the epidemic. These control functions are introduced at a specified time t with $t\in [0,T]$, as follows, where T is the final time and its value is fixed.

1. $u_1(t)$: An appropriate preventive measure on susceptible humans; appropriate personal protective equipment (PPE) (wearing rubber boots, waterproof overalls to cover wounds or skin, goggles, rubber gloves) and treat unsafe drinking water. This is to minimize the acquisition of the infection in a community.

2. $u_2(t)$: The treatment control on infectious humans (doxycycline 100mg q12h PO, Ampicillin 500-750 mg q6h PO, Amoxicillin 500 mg q6h PO, or penicillin G 1.5 MU IV q6h).

3. $u_3(t)$: Rodenticide control measure (RCM). This is to reduce the number of rodents (rats) by poisoning (or rodenticide).

4. $u_4(t)$: The control sanitation rate of the environment; improvements in slum areas (elimination of trash, improve drainage) and environmental modifications (maintaining areas around human habitations clean). This is to reduce the growth of the pathogen Leptospira. We assumed that the proportional number of rodents in each compartment of the rodent population are removed with a constant control rate of rodenticide.

The non-autonomous system of the autonomous system (3(3) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-({\lambda }_h+\mu )S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& ={\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =\delta I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-{\mu }_b{B}_l,\end{array}\end{array}$(3) ) through incorporation of the above four control measures is formulated as: (5) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& =(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +{\eta }_1{u}_2+\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =(\delta +{\eta }_1{u}_2)I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\eta }_2{u}_3+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\eta }_2{u}_3+{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\eta }_2{u}_3+{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l,\end{array}\end{array}$(5) with the initial conditions: $S_h(0)\ge 0,E_h(0)\ge 0,I_h(0)\ge 0,R_h(0)\ge 0,S_v(0)\ge 0,I_v(0)\ge 0,R_v(0)\ge 0,\mathrm{and}{B}_l(0)\ge 0$, where, ${\lambda }_h=\frac{{\beta }_2{B}_l}{\kappa +B_l}+{\beta }_1{I}_v$.

The parameters ${\eta }_1,{\eta }_2,$ and ${\eta }_3$ are the constant control rates of treatment, rodenticide and the environmental sanitation, respectively. The incidence rate of infection in humans and the recruitment rate of rodent population are reduced by the factors $(1-u_1(t))$ and $(1-u_3(t))$, respectively. The human infectious class is reduced by recovery at the rate of ${\eta }_1{u}_2(t)$. The number of rodents in each compartment of the rodent populations $S_v,I_v$ and $R_v$ are decreased at the rate ${\eta }_2{u}_3(t)$. Further, the concentration of bacterial population decreased due to control strategy ${\eta }_3{u}_4(t)$. The aim of this study is to minimize the number of infectious individuals, total rodent population, bacterial population and to minimize the costs of implementing the control interventions $u_k(t),1\le k\le 4$. For this, we use the following objective functional given by: (6) $\begin{array}{r}J(u_k)={\int }_0^T(w_1{I}_h+w_2{B}_L+w_3(S_v+I_v+R_v)+\frac{1}{2}\sum _{k=1}^4{z}_k{u}_k^2)\mathrm{d}t\end{array}$(6) subject to the non-autonomous system (5(5) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& =(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +{\eta }_1{u}_2+\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =(\delta +{\eta }_1{u}_2)I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\eta }_2{u}_3+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\eta }_2{u}_3+{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\eta }_2{u}_3+{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l,\end{array}\end{array}$(5) ), where, the coefficients $w_1,w_2$ and $w_3$ are positive weight constants of infectious individual, pathogenic population and the total number of rodent population $N_v$, respectively while the constants $z_1,z_2,z_3,z_4$ are relative cost positive weights of enforcing control interventions for each individual/vector or unit area. The terms $w_1{I}_h,w_2{B}_L$ and $w_3(S_v+I_v+R_v)$, respectively, denote the costs associated with $I_h,B_l$ and $N_v$, respectively. Further, the expressions $\frac{z_1}{2}{u}_1^2,\frac{z_2}{2}{u}_2^2,\frac{z_3}{2}{u}_3^2$ and $\frac{z_4}{2}{u}_4^2$ respectively, represent the costs of implementing the PPE, treatment control, reducing the size of $N_v$ and sanitation rate of the environment in a community. Assumed the quadratic functions for costs on the controls, as in other literatures [1, 6, 14, 26–28]. Our goal is to seek an optimal control quadruple, $u^{\ast }(t)=(u_1^{\ast }(t),u_2^{\ast }(t),u_3^{\ast }(t),u_4^{\ast }(t))$, such that quadruple (7) $\begin{array}{r}J(u_1^{\ast },u_2^{\ast },u_3^{\ast },u_4^{\ast })=\underset{(u_1,u_2,u_3,u_4)\in U}{inf}J(u_1,u_2,u_3,u_4),\end{array}$(7) where, $U=\{(u_1(t),u_2(t),u_3(t),u_4(t)):0\le u_k(t)\le 1,t\in [0,T]\}$ is the non-empty control set and $\mathrm{each}{u}_k(t)$ is Lebesgue measurable function, k=1,2,3,4.

3.1. Existence of an optimal control

Theorem 3.1

Given the objective functional J (6(6) $\begin{array}{r}J(u_k)={\int }_0^T(w_1{I}_h+w_2{B}_L+w_3(S_v+I_v+R_v)+\frac{1}{2}\sum _{k=1}^4{z}_k{u}_k^2)\mathrm{d}t\end{array}$(6) ), defined on the control set U, and subject to the model equation (5(5) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& =(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +{\eta }_1{u}_2+\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =(\delta +{\eta }_1{u}_2)I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\eta }_2{u}_3+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\eta }_2{u}_3+{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\eta }_2{u}_3+{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l,\end{array}\end{array}$(5) ), then there exists an optimal control quadruple $u^{\ast }=(u_1^{\ast },u_2^{\ast },u_3^{\ast },u_4^{\ast })$ corresponding state $(S_h^{\ast },E_h^{\ast },I_h^{\ast },R_h^{\ast },S_v^{\ast },I_v^{\ast },R_v^{\ast },B_m^{\ast })$ that holds (7(7) $\begin{array}{r}J(u_1^{\ast },u_2^{\ast },u_3^{\ast },u_4^{\ast })=\underset{(u_1,u_2,u_3,u_4)\in U}{inf}J(u_1,u_2,u_3,u_4),\end{array}$(7) ), provided that the following hypotheses are satisfied [5, 13, 26, 27, 32].

1. The admissible control set is closed and convex,

2. The right hand-side of the non-autonomous system (5(5) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& =(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +{\eta }_1{u}_2+\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =(\delta +{\eta }_1{u}_2)I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\eta }_2{u}_3+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\eta }_2{u}_3+{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\eta }_2{u}_3+{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l,\end{array}\end{array}$(5) ) is bounded by a linear function in the state and control variables,

3. The integrand of the objective functional in (6(6) $\begin{array}{r}J(u_k)={\int }_0^T(w_1{I}_h+w_2{B}_L+w_3(S_v+I_v+R_v)+\frac{1}{2}\sum _{k=1}^4{z}_k{u}_k^2)\mathrm{d}t\end{array}$(6) ) is convex with respect to the control variables,

4. There exists $J_1,J_2>0$ and $J_3>1$ such that the integrand of the objective functional is bounded below by $\begin{array}{r}{J}_1{(\sum _{k=1}^4\mid u_k\mid )}^{\frac{J_3}{2}}-J_2.\end{array}$

Proof.

Let $x=(S_h,E_h,I_h,R_h,S_v,I_v,R_v,B_l)$ be the state variables of the system (5(5) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& =(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +{\eta }_1{u}_2+\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =(\delta +{\eta }_1{u}_2)I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\eta }_2{u}_3+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\eta }_2{u}_3+{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\eta }_2{u}_3+{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l,\end{array}\end{array}$(5) ). Then, we need to verify the four hypotheses stated by Theorem 3.1.

(a) Let $U_0=\{(u_1,u_2,u_3,u_4):0\le u_k\le 1,k=1,2,3,4\}$ be the control set. Then, $U_0$ is closed and convex by definition.

(b) Let $g(x,u)$ be the right-hand side of the system (5(5) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& =(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +{\eta }_1{u}_2+\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =(\delta +{\eta }_1{u}_2)I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\eta }_2{u}_3+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\eta }_2{u}_3+{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\eta }_2{u}_3+{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l,\end{array}\end{array}$(5) ), where $u=(u_1,u_2,u_3,u_4{)}^T$. Thus, the expression $g(x,u)$ can be written as $g(x,u)=g_1(x)+g_2(x)u$, where $\begin{array}{rl}{g}_1(x)& =\left[\begin{array}{c}\mathrm{\Lambda }+\gamma R_h-{\lambda }_h{S}_h-\mu S_h\\ {\lambda }_h{S}_h-(\theta +\mu )E_h\\ \theta E_h-(\alpha +\delta +\mu )I_h\\ \delta I_h-(\gamma +\mu )R_h\\ \mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\mu }_v)S_v\\ {\beta }_3{I}_h{S}_v-(\sigma +{\mu }_v)I_v\\ \sigma I_v-(\rho +{\mu }_v)R_v\\ {\tau }_1{I}_h+{\tau }_2{I}_v-{\mu }_b{B}_l\end{array}\right],\\ g_2(x)& =\left[\begin{array}{cccc}({\displaystyle \frac{{\beta }_2{B}_l}{\kappa +B_l}}+{\beta }_1{I}_v)S_h& 0& 0& 0\\ -({\displaystyle \frac{{\beta }_2{B}_l}{\kappa +B_l}}+{\beta }_1{I}_v)S_h& 0& 0& 0\\ 0& -{\eta }_1{I}_h& 0& 0\\ 0& {\eta }_1{I}_h& 0& 0\\ 0& 0& -(\mathrm{\Pi }+{\eta }_2{S}_v)& 0\\ 0& 0& -{\eta }_2{I}_v& 0\\ 0& 0& -{\eta }_2{R}_v& 0\\ 0& 0& 0& -{\eta }_3{B}_l\end{array}\right].\end{array}$From definition of Euclidean norm [5, 32], we obtain $\begin{array}{r}\parallel g_1(x)\parallel \le \parallel g_1(x)\parallel +\parallel g_2(x)\parallel \parallel u\parallel \le m+n\parallel u\parallel ,\end{array}$where, m and n are positive constants given by, (8) $\begin{array}{r}\begin{array}{rl}m& =\sqrt{\max \{m_1,m_2,m_3,m_4,m_5\}({\mathrm{\Lambda }}^2{\mathrm{\Pi }}^2+{\mathrm{\Lambda }}^2+{\mathrm{\Pi }}^2+{\mathrm{\Lambda }}^2\mathrm{\Pi }+\mathrm{\Lambda }\mathrm{\Pi })},\mathrm{and}\\ n& =\sqrt{\max \{n_1,n_2,n_3,n_4,n_5\}({\mathrm{\Lambda }}^2{\mathrm{\Pi }}^2+{\mathrm{\Lambda }}^2+{\mathrm{\Pi }}^2+{\mathrm{\Lambda }}^2\mathrm{\Pi }+\mathrm{\Lambda }\mathrm{\Pi })},\end{array}\end{array}$(8) with (9) $\begin{array}{r}\begin{array}{rl}{m}_1& =\frac{1}{{\mu }^2}(2{\beta }_1^2+2{\eta }_1^2+{\eta }_3^2{{\tau }^{\ast }}^2),m_2=\frac{1}{{\mu }_v^2}(({\mu }_v+{\eta }_2{)}^2+2{\eta }_2^2+{\eta }_1^2{{\tau }^{\ast }}^2),\\ m_3& =\frac{{\beta }_1^2}{{\mu }^2{\mu }^2},m_4=\frac{4{\beta }_1{\beta }_2}{{\mu }_v^2\mu },m_5=\frac{2{\eta }_3^2{{\tau }^{\ast }}^2}{{\mu }_v\mu },\\ n_1& =\frac{1}{{\mu }^2}({\beta }_2^2+(\mu +\gamma {)}^2+{\theta }^2+{\delta }^2+{\tau }_1^2),n_2=\frac{1}{{\mu }_v^2}(({\mu }_v+\rho {)}^2+{\sigma }^2+{\tau }_2^2),\\ n_3& =\frac{{\beta }_1^2+{\beta }_3^2}{{\mu }^2{\mu }^2},n_4=\frac{2{\beta }_1{\beta }_2}{{\mu }_v^2\mu },n_5=\frac{2{\tau }_1{\tau }_2}{{\mu }_v\mu }.\end{array}\end{array}$(9) (c) Let $y=(I_h,B_l,S_v,R_v,I_v)$. Then, integrand of the objective functional (6(6) $\begin{array}{r}J(u_k)={\int }_0^T(w_1{I}_h+w_2{B}_L+w_3(S_v+I_v+R_v)+\frac{1}{2}\sum _{k=1}^4{z}_k{u}_k^2)\mathrm{d}t\end{array}$(6) ) is the Lagrangian of the form denoted by $L(y,u)$ defined as (10) $\begin{array}{r}L(y,u)=L_1(y)+L_2(u),\end{array}$(10) where, $L_1(y)=w_1{I}_h+w_2{B}_L+w_3(S_v+I_v+R_v),\mathrm{and}{L}_2(u)=\frac{1}{2}\sum _{k=1}^4{z}_k{u}_k^2$. Note that $L_2(u)$ is a finite linear combination of the non-negative control functions ${\psi }_k=\frac{1}{2}{u}_k^2,k=1,2,3,4$. Thus, $L_2(u)$ is convex on the control variable u. (d) Since the Lagrangian in (10(10) $\begin{array}{r}L(y,u)=L_1(y)+L_2(u),\end{array}$(10) ) is the sum of non-negative terms $L_1(y)$ and $L_2(u)$, the last hypothesis is shown as follow: (11) $\begin{array}{r}L(y,u)=L_1(y)+\frac{1}{2}\sum _{k=1}^4{z}_k{u}_k^2\ge \frac{1}{2}\sum _{k=1}^4{z}_k{u}_k^2\ge \frac{J_1}{2}{(\sum _{k=1}^4\mid u_k{\mid }^2)}^{\frac{J_3}{2}}-J_2,\end{array}$(11) where, $J_1=\min \{z_1,z_2,z_3,z_4\},J_2>0\mathrm{and}{J}_3=2.$

3.2. Characterization of the optimal control

To solve numerical solutions of our standard optimal problem, we need to generate necessary conditions that the optimal control quadruple and state must satisfy. Such conditions are derived from an associated problem of minimizing point-wise with respect to controls $u_k,k=1,2,3,4$, a Hamiltonian H, of the optimal problem. The Hamiltonian H, is formulated from the system (5(5) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& =(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +{\eta }_1{u}_2+\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =(\delta +{\eta }_1{u}_2)I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\eta }_2{u}_3+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\eta }_2{u}_3+{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\eta }_2{u}_3+{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l,\end{array}\end{array}$(5) ) and from the objective functional (6(6) $\begin{array}{r}J(u_k)={\int }_0^T(w_1{I}_h+w_2{B}_L+w_3(S_v+I_v+R_v)+\frac{1}{2}\sum _{k=1}^4{z}_k{u}_k^2)\mathrm{d}t\end{array}$(6) ) using the method presented in Pontryagin's Maximum Principle (PMP) [22], which is based on PMP [30], as follow: (12) $\begin{array}{r}\begin{array}{rl}H& =w_1{I}_h+w_2{B}_L+w_3(S_v+I_v+R_v)\\ & +\frac{1}{2}\sum _{k=1}^4{z}_k{u}_k^2+{\lambda }_1[\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h]\\ & +{\lambda }_2[(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h]\\ & +{\lambda }_3[\theta E_h-(\alpha +\delta +\eta u_2+\mu )I_h]\\ & +{\lambda }_4[(\delta +\eta u_2)I_h-(\gamma +\mu )R_h]\\ & +{\lambda }_5[(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\mu }_v+{\eta }_2{u}_3)S_v]\\ & +{\lambda }_6[{\beta }_3{I}_h{S}_v-(\sigma +{\mu }_v+{\eta }_2{u}_3)I_v]\\ & +{\lambda }_7[\sigma I_v-(\rho +{\mu }_v+{\eta }_2{u}_3)R_v]\\ & +{\lambda }_8[{\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l],\end{array}\end{array}$(12) where, ${\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5{\lambda }_6,{\lambda }_7$ and ${\lambda }_8$ are the adjoint (co-state) variables which determine the adjoint system and are associated for the states $S_h,E_h,I_h,R_h,S_v,I_v,R_v$ and $B_l$ respectively. Then, we apply PMP to the Hamiltonian (12(12) $\begin{array}{r}\begin{array}{rl}H& =w_1{I}_h+w_2{B}_L+w_3(S_v+I_v+R_v)\\ & +\frac{1}{2}\sum _{k=1}^4{z}_k{u}_k^2+{\lambda }_1[\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h]\\ & +{\lambda }_2[(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h]\\ & +{\lambda }_3[\theta E_h-(\alpha +\delta +\eta u_2+\mu )I_h]\\ & +{\lambda }_4[(\delta +\eta u_2)I_h-(\gamma +\mu )R_h]\\ & +{\lambda }_5[(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\mu }_v+{\eta }_2{u}_3)S_v]\\ & +{\lambda }_6[{\beta }_3{I}_h{S}_v-(\sigma +{\mu }_v+{\eta }_2{u}_3)I_v]\\ & +{\lambda }_7[\sigma I_v-(\rho +{\mu }_v+{\eta }_2{u}_3)R_v]\\ & +{\lambda }_8[{\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l],\end{array}\end{array}$(12) ) to generate the necessary conditions for our control problem. Further, the PMP [22] and the existence of the optimal control from (see Theorem 4.1, [13]) can be used to obtain the following theorem.

Theorem 3.2

If $X^{\ast }=(S_h^{\ast },E_h^{\ast },I_h^{\ast },R_h^{\ast },S_v^{\ast },I_v^{\ast },R_v^{\ast },B_m^{\ast })$ is an optimal state of the corresponding non-autonomous system (5(5) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& =(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +{\eta }_1{u}_2+\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =(\delta +{\eta }_1{u}_2)I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\eta }_2{u}_3+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\eta }_2{u}_3+{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\eta }_2{u}_3+{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l,\end{array}\end{array}$(5) ) and $U^{\ast }=(u_1^{\ast },u_2^{\ast },u_3^{\ast },u_4^{\ast })$ is an optimal control quadruple such that (7(7) $\begin{array}{r}J(u_1^{\ast },u_2^{\ast },u_3^{\ast },u_4^{\ast })=\underset{(u_1,u_2,u_3,u_4)\in U}{inf}J(u_1,u_2,u_3,u_4),\end{array}$(7) ) holds, then there exist adjoint variables ${\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5,{\lambda }_6,{\lambda }_7,{\lambda }_8$ satisfying (13) $\begin{array}{rl}\frac{\mathrm{d}{\lambda }_1}{\mathrm{d}t}& =({\lambda }_1-{\lambda }_2)(1-u_1)(\frac{{\beta }_2{B}_l}{\kappa +B_l}+{\beta }_1{I}_v)+{\lambda }_1\mu ,\\ \frac{\mathrm{d}{\lambda }_2}{\mathrm{d}t}& =({\lambda }_2-{\lambda }_3)\theta +{\lambda }_2\mu ,\\ \frac{\mathrm{d}{\lambda }_3}{\mathrm{d}t}& =-w_1+{\lambda }_3(\alpha +\delta +{\eta }_1{u}_2+\mu )-{\lambda }_4(\delta +{\eta }_1{u}_2)+({\lambda }_5-{\lambda }_6){\beta }_3{S}_v-{\lambda }_8{\tau }_1,\\ \frac{\mathrm{d}{\lambda }_4}{\mathrm{d}t}& =({\lambda }_4-{\lambda }_1)\gamma +{\lambda }_4\mu ,\\ \frac{\mathrm{d}{\lambda }_5}{\mathrm{d}t}& =-w_3+({\lambda }_5-{\lambda }_6){\beta }_3{I}_h+{\lambda }_5({\mu }_v+{\eta }_2{u}_3),\\ \frac{\mathrm{d}{\lambda }_6}{\mathrm{d}t}& =-w_3+({\lambda }_1-{\lambda }_2)(1-u_1){\beta }_1{S}_h+{\lambda }_6(\sigma +{\mu }_v+{\eta }_2{u}_3)-{\lambda }_7\sigma -{\tau }_2{\lambda }_8,\\ \frac{\mathrm{d}{\lambda }_7}{\mathrm{d}t}& =-w_3+({\lambda }_7-{\lambda }_5)\rho +{\lambda }_7({\mu }_v+{\eta }_2{u}_3),\\ \frac{\mathrm{d}{\lambda }_8}{\mathrm{d}t}& =-w_2+({\lambda }_1-{\lambda }_2)(1-u_1)\frac{{\beta }_2\kappa S_h}{(\kappa +B_l{)}^2}+{\lambda }_8({\eta }_3{u}_4+{\mu }_b),\end{array}$(13) and with transversality conditions (or final time conditions), (14) $\begin{array}{r}{\lambda }_k(T)=0,k=1,2,\dots ,8(\text{where T is the final time}).\end{array}$(14) Furthermore, the optimal controls $u_i^{\ast },i=1,2,3,4$ that minimizes the objective function $J(u_1,u_2,u_3,u_4)$ over U are characterized by (15) $\begin{array}{r}\begin{array}{rl}{u}_1^{\ast }& =min\{1,max\{\frac{({\lambda }_2-{\lambda }_1)S_h(\frac{{\beta }_2{B}_l}{\kappa +B_l}+{\beta }_1{I}_v)}{z_1},0\}\},\\ u_2^{\ast }& =min\{1,max\{\frac{({\lambda }_3-{\lambda }_4){\eta }_1{I}_h}{z_2},0\}\},\\ u_3^{\ast }& =min\{1,max\{\frac{{\lambda }_5\mathrm{\Pi }+{\eta }_2({\lambda }_5{S}_v+{\lambda }_6{I}_v+{\lambda }_7{R}_v)}{z_3},0\}\},\\ u_4^{\ast }& =min\{1,max\{\frac{{\lambda }_8{\eta }_3{B}_l}{z_4},0\}\}.\end{array}\end{array}$(15)

Proof.

We apply the standard results in [30] (PMP) to derive the adjoint relations, the transversality conditions and the optimal control quadruple.

Thus, by taking negatives of partial derivatives of the Hamiltonian function H given in (12(12) $\begin{array}{r}\begin{array}{rl}H& =w_1{I}_h+w_2{B}_L+w_3(S_v+I_v+R_v)\\ & +\frac{1}{2}\sum _{k=1}^4{z}_k{u}_k^2+{\lambda }_1[\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h]\\ & +{\lambda }_2[(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h]\\ & +{\lambda }_3[\theta E_h-(\alpha +\delta +\eta u_2+\mu )I_h]\\ & +{\lambda }_4[(\delta +\eta u_2)I_h-(\gamma +\mu )R_h]\\ & +{\lambda }_5[(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\mu }_v+{\eta }_2{u}_3)S_v]\\ & +{\lambda }_6[{\beta }_3{I}_h{S}_v-(\sigma +{\mu }_v+{\eta }_2{u}_3)I_v]\\ & +{\lambda }_7[\sigma I_v-(\rho +{\mu }_v+{\eta }_2{u}_3)R_v]\\ & +{\lambda }_8[{\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l],\end{array}\end{array}$(12) ) with respect to the associated state variables $S_h^{\ast },E_h^{\ast },I_h^{\ast },R_h^{\ast },S_v^{\ast },I_v^{\ast },R_v^{\ast },B_m^{\ast }$ respectively, yields the adjoint equations (13(13) $\begin{array}{rl}\frac{\mathrm{d}{\lambda }_1}{\mathrm{d}t}& =({\lambda }_1-{\lambda }_2)(1-u_1)(\frac{{\beta }_2{B}_l}{\kappa +B_l}+{\beta }_1{I}_v)+{\lambda }_1\mu ,\\ \frac{\mathrm{d}{\lambda }_2}{\mathrm{d}t}& =({\lambda }_2-{\lambda }_3)\theta +{\lambda }_2\mu ,\\ \frac{\mathrm{d}{\lambda }_3}{\mathrm{d}t}& =-w_1+{\lambda }_3(\alpha +\delta +{\eta }_1{u}_2+\mu )-{\lambda }_4(\delta +{\eta }_1{u}_2)+({\lambda }_5-{\lambda }_6){\beta }_3{S}_v-{\lambda }_8{\tau }_1,\\ \frac{\mathrm{d}{\lambda }_4}{\mathrm{d}t}& =({\lambda }_4-{\lambda }_1)\gamma +{\lambda }_4\mu ,\\ \frac{\mathrm{d}{\lambda }_5}{\mathrm{d}t}& =-w_3+({\lambda }_5-{\lambda }_6){\beta }_3{I}_h+{\lambda }_5({\mu }_v+{\eta }_2{u}_3),\\ \frac{\mathrm{d}{\lambda }_6}{\mathrm{d}t}& =-w_3+({\lambda }_1-{\lambda }_2)(1-u_1){\beta }_1{S}_h+{\lambda }_6(\sigma +{\mu }_v+{\eta }_2{u}_3)-{\lambda }_7\sigma -{\tau }_2{\lambda }_8,\\ \frac{\mathrm{d}{\lambda }_7}{\mathrm{d}t}& =-w_3+({\lambda }_7-{\lambda }_5)\rho +{\lambda }_7({\mu }_v+{\eta }_2{u}_3),\\ \frac{\mathrm{d}{\lambda }_8}{\mathrm{d}t}& =-w_2+({\lambda }_1-{\lambda }_2)(1-u_1)\frac{{\beta }_2\kappa S_h}{(\kappa +B_l{)}^2}+{\lambda }_8({\eta }_3{u}_4+{\mu }_b),\end{array}$(13) ): (16) $\begin{array}{r}\begin{array}{rl}\frac{\mathrm{d}{\lambda }_1}{\mathrm{d}t}& =-\frac{\mathrm{\partial }H}{\mathrm{\partial }{S}_h},\frac{\mathrm{d}{\lambda }_2}{\mathrm{d}t}=-\frac{\mathrm{\partial }H}{\mathrm{\partial }{E}_h},\frac{\mathrm{d}{\lambda }_3}{\mathrm{d}t}=-\frac{\mathrm{\partial }H}{\mathrm{\partial }{I}_h},\frac{\mathrm{d}{\lambda }_4}{\mathrm{d}t}=-\frac{\mathrm{\partial }H}{\mathrm{\partial }{R}_h},\\ \frac{\mathrm{d}{\lambda }_5}{\mathrm{d}t}& =-\frac{\mathrm{\partial }H}{\mathrm{\partial }{S}_v},\frac{\mathrm{d}{\lambda }_6}{\mathrm{d}t}=-\frac{\mathrm{\partial }H}{\mathrm{\partial }{I}_v},\frac{\mathrm{d}{\lambda }_7}{\mathrm{d}t}=-\frac{\mathrm{\partial }H}{\mathrm{\partial }{R}_v},\frac{\mathrm{d}{\lambda }_8}{\mathrm{d}t}=-\frac{\mathrm{\partial }H}{\mathrm{\partial }{B}_l},\end{array}\end{array}$(16) and with transversality conditions, ${\lambda }_k(T)=0,k=1,2,\dots ,8$ (see the proof the these conditions from theorem1 in [22]). Finally, the optimal characterization descried in (15(15) $\begin{array}{r}\begin{array}{rl}{u}_1^{\ast }& =min\{1,max\{\frac{({\lambda }_2-{\lambda }_1)S_h(\frac{{\beta }_2{B}_l}{\kappa +B_l}+{\beta }_1{I}_v)}{z_1},0\}\},\\ u_2^{\ast }& =min\{1,max\{\frac{({\lambda }_3-{\lambda }_4){\eta }_1{I}_h}{z_2},0\}\},\\ u_3^{\ast }& =min\{1,max\{\frac{{\lambda }_5\mathrm{\Pi }+{\eta }_2({\lambda }_5{S}_v+{\lambda }_6{I}_v+{\lambda }_7{R}_v)}{z_3},0\}\},\\ u_4^{\ast }& =min\{1,max\{\frac{{\lambda }_8{\eta }_3{B}_l}{z_4},0\}\}.\end{array}\end{array}$(15) ) is obtained by setting partial derivatives of the H, with respect to the control measures $u_1,u_2,u_3$ and $u_4$ to zero (optimal condition) [30] in the interior of the control set U. Thus, we have (17) $\begin{array}{r}\frac{\mathrm{\partial }H}{\mathrm{\partial }{u}_i}=0,\text{ for }{u}_i^{\ast }(\text{where }i=1,2,3,4)(\mathbf{\text{optimal condition}}).\end{array}$(17) Solving the Equation (17(17) $\begin{array}{r}\frac{\mathrm{\partial }H}{\mathrm{\partial }{u}_i}=0,\text{ for }{u}_i^{\ast }(\text{where }i=1,2,3,4)(\mathbf{\text{optimal condition}}).\end{array}$(17) ) for each optimal control $u_1^{\ast },u_2^{\ast },u_3^{\ast }$ and $u_4^{\ast }$ we obtain the characterization of optimal controls which is the same as descried in (15(15) $\begin{array}{r}\begin{array}{rl}{u}_1^{\ast }& =min\{1,max\{\frac{({\lambda }_2-{\lambda }_1)S_h(\frac{{\beta }_2{B}_l}{\kappa +B_l}+{\beta }_1{I}_v)}{z_1},0\}\},\\ u_2^{\ast }& =min\{1,max\{\frac{({\lambda }_3-{\lambda }_4){\eta }_1{I}_h}{z_2},0\}\},\\ u_3^{\ast }& =min\{1,max\{\frac{{\lambda }_5\mathrm{\Pi }+{\eta }_2({\lambda }_5{S}_v+{\lambda }_6{I}_v+{\lambda }_7{R}_v)}{z_3},0\}\},\\ u_4^{\ast }& =min\{1,max\{\frac{{\lambda }_8{\eta }_3{B}_l}{z_4},0\}\}.\end{array}\end{array}$(15) ). Let ${\varphi }_1=\frac{({\lambda }_2-{\lambda }_1)S_h(\frac{{\beta }_2{B}_l}{\kappa +B_l}+{\beta }_1{I}_v)}{z_1},{\varphi }_2=\frac{({\lambda }_3-{\lambda }_4){\eta }_1{I}_h}{z_2},{\varphi }_3=\frac{{\lambda }_5\mathrm{\Pi }+{\eta }_2({\lambda }_5{S}_v+{\lambda }_6{I}_v+{\lambda }_7{R}_v)}{z_3},{\varphi }_4=\frac{{\lambda }_8{\eta }_3{B}_l}{z_4}$.

Therefore, the optimal controls $u_i^{\ast },i=1,2,3,4$ are given by in compact form as (18) $\begin{array}{r}\begin{array}{rl}{u}_1^{\ast }& =\min \{1,\max \{{\varphi }_1,0\}\},u_2^{\ast }=\min \{1,\max \{{\varphi }_2,0\}\},\\ u_3^{\ast }& =\min \{1,\max \{{\varphi }_3,0\}\},u_4^{\ast }=\min \{1,\max \{{\varphi }_4,0\}\}.\end{array}\end{array}$(18)

Remark 3.1

1. $\frac{d{x}_k}{\mathrm{d}t}=\frac{\mathrm{\partial }H}{\mathrm{\partial }{\lambda }_k}\text{ from the equation~(12)}$ (where, $x_k$ are the state variables and ${\lambda }_k$ are the associated adjoint variables $k=1,2,3,\dots ,8$).

2. The optimality system of our model consists of the non-autonomous system (5(5) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& =(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +{\eta }_1{u}_2+\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =(\delta +{\eta }_1{u}_2)I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\eta }_2{u}_3+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\eta }_2{u}_3+{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\eta }_2{u}_3+{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l,\end{array}\end{array}$(5) ), adjoint equations (13(13) $\begin{array}{rl}\frac{\mathrm{d}{\lambda }_1}{\mathrm{d}t}& =({\lambda }_1-{\lambda }_2)(1-u_1)(\frac{{\beta }_2{B}_l}{\kappa +B_l}+{\beta }_1{I}_v)+{\lambda }_1\mu ,\\ \frac{\mathrm{d}{\lambda }_2}{\mathrm{d}t}& =({\lambda }_2-{\lambda }_3)\theta +{\lambda }_2\mu ,\\ \frac{\mathrm{d}{\lambda }_3}{\mathrm{d}t}& =-w_1+{\lambda }_3(\alpha +\delta +{\eta }_1{u}_2+\mu )-{\lambda }_4(\delta +{\eta }_1{u}_2)+({\lambda }_5-{\lambda }_6){\beta }_3{S}_v-{\lambda }_8{\tau }_1,\\ \frac{\mathrm{d}{\lambda }_4}{\mathrm{d}t}& =({\lambda }_4-{\lambda }_1)\gamma +{\lambda }_4\mu ,\\ \frac{\mathrm{d}{\lambda }_5}{\mathrm{d}t}& =-w_3+({\lambda }_5-{\lambda }_6){\beta }_3{I}_h+{\lambda }_5({\mu }_v+{\eta }_2{u}_3),\\ \frac{\mathrm{d}{\lambda }_6}{\mathrm{d}t}& =-w_3+({\lambda }_1-{\lambda }_2)(1-u_1){\beta }_1{S}_h+{\lambda }_6(\sigma +{\mu }_v+{\eta }_2{u}_3)-{\lambda }_7\sigma -{\tau }_2{\lambda }_8,\\ \frac{\mathrm{d}{\lambda }_7}{\mathrm{d}t}& =-w_3+({\lambda }_7-{\lambda }_5)\rho +{\lambda }_7({\mu }_v+{\eta }_2{u}_3),\\ \frac{\mathrm{d}{\lambda }_8}{\mathrm{d}t}& =-w_2+({\lambda }_1-{\lambda }_2)(1-u_1)\frac{{\beta }_2\kappa S_h}{(\kappa +B_l{)}^2}+{\lambda }_8({\eta }_3{u}_4+{\mu }_b),\end{array}$(13) ) and the terminal conditions (14(14) $\begin{array}{r}{\lambda }_k(T)=0,k=1,2,\dots ,8(\text{where T is the final time}).\end{array}$(14) ) together with the characterization of the optimal controls (18(18) $\begin{array}{r}\begin{array}{rl}{u}_1^{\ast }& =\min \{1,\max \{{\varphi }_1,0\}\},u_2^{\ast }=\min \{1,\max \{{\varphi }_2,0\}\},\\ u_3^{\ast }& =\min \{1,\max \{{\varphi }_3,0\}\},u_4^{\ast }=\min \{1,\max \{{\varphi }_4,0\}\}.\end{array}\end{array}$(18) ).

3.3. Numerical simulations

In this section, we performed numerical simulations of the optimality system of the leptospirosis epidemic (5(5) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& =(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +{\eta }_1{u}_2+\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =(\delta +{\eta }_1{u}_2)I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\eta }_2{u}_3+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\eta }_2{u}_3+{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\eta }_2{u}_3+{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l,\end{array}\end{array}$(5) ), (13(13) $\begin{array}{rl}\frac{\mathrm{d}{\lambda }_1}{\mathrm{d}t}& =({\lambda }_1-{\lambda }_2)(1-u_1)(\frac{{\beta }_2{B}_l}{\kappa +B_l}+{\beta }_1{I}_v)+{\lambda }_1\mu ,\\ \frac{\mathrm{d}{\lambda }_2}{\mathrm{d}t}& =({\lambda }_2-{\lambda }_3)\theta +{\lambda }_2\mu ,\\ \frac{\mathrm{d}{\lambda }_3}{\mathrm{d}t}& =-w_1+{\lambda }_3(\alpha +\delta +{\eta }_1{u}_2+\mu )-{\lambda }_4(\delta +{\eta }_1{u}_2)+({\lambda }_5-{\lambda }_6){\beta }_3{S}_v-{\lambda }_8{\tau }_1,\\ \frac{\mathrm{d}{\lambda }_4}{\mathrm{d}t}& =({\lambda }_4-{\lambda }_1)\gamma +{\lambda }_4\mu ,\\ \frac{\mathrm{d}{\lambda }_5}{\mathrm{d}t}& =-w_3+({\lambda }_5-{\lambda }_6){\beta }_3{I}_h+{\lambda }_5({\mu }_v+{\eta }_2{u}_3),\\ \frac{\mathrm{d}{\lambda }_6}{\mathrm{d}t}& =-w_3+({\lambda }_1-{\lambda }_2)(1-u_1){\beta }_1{S}_h+{\lambda }_6(\sigma +{\mu }_v+{\eta }_2{u}_3)-{\lambda }_7\sigma -{\tau }_2{\lambda }_8,\\ \frac{\mathrm{d}{\lambda }_7}{\mathrm{d}t}& =-w_3+({\lambda }_7-{\lambda }_5)\rho +{\lambda }_7({\mu }_v+{\eta }_2{u}_3),\\ \frac{\mathrm{d}{\lambda }_8}{\mathrm{d}t}& =-w_2+({\lambda }_1-{\lambda }_2)(1-u_1)\frac{{\beta }_2\kappa S_h}{(\kappa +B_l{)}^2}+{\lambda }_8({\eta }_3{u}_4+{\mu }_b),\end{array}$(13) ), (14(14) $\begin{array}{r}{\lambda }_k(T)=0,k=1,2,\dots ,8(\text{where T is the final time}).\end{array}$(14) ) & (18(18) $\begin{array}{r}\begin{array}{rl}{u}_1^{\ast }& =\min \{1,\max \{{\varphi }_1,0\}\},u_2^{\ast }=\min \{1,\max \{{\varphi }_2,0\}\},\\ u_3^{\ast }& =\min \{1,\max \{{\varphi }_3,0\}\},u_4^{\ast }=\min \{1,\max \{{\varphi }_4,0\}\}.\end{array}\end{array}$(18) ) using the Forward – Backward Sweep method (FBSM) implemented in Matlab programme see details in the book by Lenhart and Workman [22]. We first solve the system (5(5) $\begin{array}{r}\{\begin{array}{rl}\frac{d{S}_h}{\mathrm{d}t}& =\mathrm{\Lambda }+\gamma R_h-(1-u_1){\lambda }_h{S}_h-\mu S_h,\\ \frac{d{E}_h}{\mathrm{d}t}& =(1-u_1){\lambda }_h{S}_h-(\theta +\mu )E_h,\\ \frac{d{I}_h}{\mathrm{d}t}& =\theta E_h-(\alpha +\delta +{\eta }_1{u}_2+\mu )I_h,\\ \frac{d{R}_h}{\mathrm{d}t}& =(\delta +{\eta }_1{u}_2)I_h-(\gamma +\mu )R_h,\\ \frac{d{S}_v}{\mathrm{d}t}& =(1-u_3)\mathrm{\Pi }+\rho R_v-({\beta }_3{I}_h+{\eta }_2{u}_3+{\mu }_v)S_v,\\ \frac{d{I}_v}{\mathrm{d}t}& ={\beta }_3{I}_h{S}_v-(\sigma +{\eta }_2{u}_3+{\mu }_v)I_v,\\ \frac{d{R}_v}{\mathrm{d}t}& =\sigma I_v-(\rho +{\eta }_2{u}_3+{\mu }_v)R_v,\\ \frac{d{B}_l}{\mathrm{d}t}& ={\tau }_1{I}_h+{\tau }_2{I}_v-({\eta }_3{u}_4+{\mu }_b)B_l,\end{array}\end{array}$(5) ) with an initial guess for the control variables over simulated time [0, 400] using forward RKM-4. Then, due to the terminal conditions (14(14) $\begin{array}{r}{\lambda }_k(T)=0,k=1,2,\dots ,8(\text{where T is the final time}).\end{array}$(14) ) we use backward RKM-4 method in the time to solve the adjoint system (13(13) $\begin{array}{rl}\frac{\mathrm{d}{\lambda }_1}{\mathrm{d}t}& =({\lambda }_1-{\lambda }_2)(1-u_1)(\frac{{\beta }_2{B}_l}{\kappa +B_l}+{\beta }_1{I}_v)+{\lambda }_1\mu ,\\ \frac{\mathrm{d}{\lambda }_2}{\mathrm{d}t}& =({\lambda }_2-{\lambda }_3)\theta +{\lambda }_2\mu ,\\ \frac{\mathrm{d}{\lambda }_3}{\mathrm{d}t}& =-w_1+{\lambda }_3(\alpha +\delta +{\eta }_1{u}_2+\mu )-{\lambda }_4(\delta +{\eta }_1{u}_2)+({\lambda }_5-{\lambda }_6){\beta }_3{S}_v-{\lambda }_8{\tau }_1,\\ \frac{\mathrm{d}{\lambda }_4}{\mathrm{d}t}& =({\lambda }_4-{\lambda }_1)\gamma +{\lambda }_4\mu ,\\ \frac{\mathrm{d}{\lambda }_5}{\mathrm{d}t}& =-w_3+({\lambda }_5-{\lambda }_6){\beta }_3{I}_h+{\lambda }_5({\mu }_v+{\eta }_2{u}_3),\\ \frac{\mathrm{d}{\lambda }_6}{\mathrm{d}t}& =-w_3+({\lambda }_1-{\lambda }_2)(1-u_1){\beta }_1{S}_h+{\lambda }_6(\sigma +{\mu }_v+{\eta }_2{u}_3)-{\lambda }_7\sigma -{\tau }_2{\lambda }_8,\\ \frac{\mathrm{d}{\lambda }_7}{\mathrm{d}t}& =-w_3+({\lambda }_7-{\lambda }_5)\rho +{\lambda }_7({\mu }_v+{\eta }_2{u}_3),\\ \frac{\mathrm{d}{\lambda }_8}{\mathrm{d}t}& =-w_2+({\lambda }_1-{\lambda }_2)(1-u_1)\frac{{\beta }_2\kappa S_h}{(\kappa +B_l{)}^2}+{\lambda }_8({\eta }_3{u}_4+{\mu }_b),\end{array}$(13) ) by using the final time conditions for adjoints and the values from previous step for controls and states. Finally, we updated the control values by averaging the previous value and new value from the characterization of controls (18(18) $\begin{array}{r}\begin{array}{rl}{u}_1^{\ast }& =\min \{1,\max \{{\varphi }_1,0\}\},u_2^{\ast }=\min \{1,\max \{{\varphi }_2,0\}\},\\ u_3^{\ast }& =\min \{1,\max \{{\varphi }_3,0\}\},u_4^{\ast }=\min \{1,\max \{{\varphi }_4,0\}\}.\end{array}\end{array}$(18) ) and we use to solve the state and then the adjoint system. This process is repeated until the required convergence occur for the current state, adjoint, and control values. The parameters values used in the simulations are given in Table 1. Those parameters values have taken from previously published papers and others are estimated as shown in the Table 1. The parameters values in the Table 1 yields the basic reproduction number $R_0\approx 2.8892>1$. The initial conditions of the state variables are chosen as $S_h(0)=270,S_v(0)=510,E(0)=20,I_h(0)=10,I_v(0)=10,R_h(0)=R_v(0)=0,B_l(0)=100$. In addition, the control rate and weight constant values are chosen as; $w_1=7,w_2=1,w_3=7,z_2=7,z_2=7,z_3=10,z_4=7$. Next, we consider different control strategies to elucidate their effect on the transmission dynamics of leptospirosis epidemic. The following five most effective optimal strategies are implemented and presented graphically with and without the strategies.

• Strategy A: Rodenticide control measure only (RCM) $(u_3\ne 0,u_1=u_2=u_4=0)$,

• Strategy B: Combination of personal protective equipment (PPE) and RCM controls $(u_1\ne 0,u_3\ne 0,u_2=u_4=0)$,

• Strategy C: Combination PPE, treatment and RCM controls $(u_1\ne 0,u_2\ne 0,u_3\ne 0,u_4=0)$,

• Strategy D: Combination PPE, RCM and improvement of slum environment controls $(u_1\ne 0,u_3\ne 0,u_4\ne 0,u_2=0)$,

• Strategy E: Combination of all controls measures $(u_1\ne 0,u_2\ne 0,u_3\ne 0,u_4\ne 0)$.

Table 1. Description of parameter values used in the simulations.

Parameter	Description	Value	Unit	Source
Λ	The recruitment rate of susceptible humans	$\mu \times N0$	${\text{Humans Day}}^{-1}$	Assumed,
Π	The recruitment rate of susceptible rodents	2	${\text{Humans day}}^{-1}$	[11, 28]
${\beta }_1$	The transmission coefficient between $S_h$ and $I_v$	$0.00033$	$(\text{Rodents day}{)}^{-1}$	Assumed
${\beta }_2$	The transmission coefficient between $S_h$ and $B_b$	$0.0815$	${\mathrm{Day}}^{-1}$	Assumed
${\beta }_3$	The transmission coefficient between $S_v$ and $I_h$	$0.0007$	$(\text{Humans day}{)}^{-1}$	Assumed
μ	The natural death rate of susceptible humans	$0.000046$	${\mathrm{Day}}^{-1}$	[18, 25]
α	The disease related death rate	$0.04$	${\mathrm{Day}}^{-1}$	Assumed
θ	The rate at which latently infected individual
	transfer to infectious	$0.092$	${\mathrm{Day}}^{-1}$	[11, 20]
γ	Leptospirosis disease waning immunity from $R_h$	$0.089$	${\mathrm{Day}}^{-1}$	[11, 18]
δ	The rate at which infectious individual transfer
	to human recovery class	$0.072$	${\mathrm{Day}}^{-1}$	Assumed
σ	The rate of recovery from infectious rodents	$0.064$	${\mathrm{Day}}^{-1}$	Assumed
ρ	Leptospirosis disease waning immunity from $R_v$	$0.083$	${\mathrm{Day}}^{-1}$	Assumed
${\mu }_v$	The natural death rate of the vectors	$0.0029$	${\mathrm{Day}}^{-1}$	[3, 11]
${\mu }_b$	The natural death rate of the bacteria	$0.05$	${\mathrm{Day}}^{-1}$	[11, 24]
${\tau }_1$	The rate of increase of bacterial population due to $I_h$	$0.06$	${\displaystyle \frac{\text{No. of pathogen }}{\text{Human day}}}$	Assumed
${\tau }_2$	The rate of increase of bacterial population due to $I_v$	$0.2$	${\displaystyle \frac{\text{No. of pathogen}}{\text{Rodents day}}}$	Assumed
κ	The concentration of pathogenic population	10000	$\text{No. of pathogen}$	[11]
${\eta }_1$	The constant control rate of treatment	$0.09$	Dimensionless	Assumed
${\eta }_2$	The constant control rate of rodenticide	$0.2$	Dimensionless	[11, 24]
${\eta }_3$	The constant environmental sanitation control rate	$0.09$	Dimensionless	Assumed

Over the simulation period, the significant effects of all the five strategies on disease spread are illustrated in Figures 1–5. It can be seen for each of the different strategies implemented that in Figures 1(a), 2(a), 3(a), 4(a), and 5 (a ) the number of susceptible individuals in the population decreases more rapidly without optimal controls. On the other hand, Figures 1(b) – 1(e), 2(b) – 2(e), 3(b) – 3(e), 4(b) – 4(e), and 5(b) – 5(e) indicated that the number size of infected classes in human population ($E_h\mathrm{\&}{I}_h$), the size of bacterial population, and the rodent classes, respectively decrease more sharply in presence of control strategies than without control strategies. Further, the control profiles for all strategies are shown in Figures 6(a) – 6(e). We assume that all control measures would be attained the maximum effort (100%). In particular, the control profiles for strategy E in Figure 6 (e) revealed that the optimal use of $u_1$, $u_2$, $u_3$ and $u_4$ are at maximum effort $(100\mathrm{\%})$ for the first 15, 25, 399.2, and 23 days, respectively and then each control reduces slowly to the lower bound. The control profile in Figure 6 (a) revealed that the optimal use of the RCM is at the upper bound for 399.2 days and then gradually reduced to lower bound for the rest of simulation time. Similar, conclusions can be drawn for control profiles of the remaining optimal control strategies.

Figure 1. The Dynamics of Leptospirosis Disease using the Strategy A ($u_3\ne 0,u_1=u_2=u_4=0$). (a) Variation of susceptible individuals with and without control Strategy A. (b) Variation of infected individuals $(E_h\mathrm{\&}{I}_h)$ with and without control Strategy A. (c) Variation of bacterial population with and without control Strategy A. (d) Variation of susceptible rodents with and without control Strategy A. (e) Variation of infected and recovered rodents with and without control Strategy A.

Figure 2. The Dynamics of Leptospirosis Disease using the Strategy B ($u_1\ne 0,u_3\ne 0,u_2=u_4=0$). (a) Variation of susceptible individuals with and without control Strategy B. (b) Variation of infected individuals $(E_h\mathrm{\&}{I}_h)$ with and without control Strategy B. (c) Variation of bacterial population with and without control Strategy B. (d) Variation of susceptible rodents with and without control Strategy B. (e) Variation of infected and recovered rodents with and without control Strategy B.

Figure 3. The Dynamics of Leptospirosis Disease using the Strategy C ($u_1\ne 0,u_2\ne 0,u_3\ne 0,u_4=0$). (a) Variation of susceptible individuals with and without control Strategy C. (b) Variation of infected individuals $(E_h\mathrm{\&}{I}_h)$ with and without control Strategy C. (c) Variation of bacterial population with and without control Strategy C. (d) Variation of susceptible rodents with and without control Strategy C. (e) Variation of infected and recovered rodents with and without control Strategy C.

Figure 4. The Dynamics of Leptospirosis Disease using the Strategy D ($u_1\ne 0,u_3\ne 0,u_4\ne 0,u_2=0$). (a) Variation of susceptible individuals with and without control Strategy D. (b) Variation of infected individuals $(E_h\mathrm{\&}{I}_h)$ with and without control Strategy D. (c) Variation of bacterial population with and without control Strategy D. (d) Variation of susceptible rodents with and without control Strategy D. (e) Variation of infected and recovered rodents with and without control Strategy D.

Figure 5. The Dynamics of Leptospirosis Disease using the Strategy E ($u_1\ne 0,u_2\ne 0,u_3\ne 0,u_4\ne 0$). (a) Variation of susceptible individuals with and without control Strategy E. (b) Variation of infected individuals $(E_h\mathrm{\&}{I}_h)$ with and without control Strategy E. (c) Variation of bacterial population with and without control Strategy E. (d) Variation of susceptible rodents with and without control Strategy E. (e) Variation of infected and recovered rodents with and without control Strategy E.

Figure 6. Plots showing the Control profiles of the optimal control Strategies (A-E). (a) Control profile of the Strategy A. (b) Control profile of the Strategy B. (c) Control profile of the Strategy C. (d) Control profile of the Strategy D. (e) Control profile of the Strategy E.

4. Cost-effectiveness analysis

In this section, we implemented two approaches namely, the average cost-effectiveness ratio $(\mathrm{ACER})$ and the incremental cost-effectiveness ratio $(\mathrm{ICER})$ to identify the most cost-effective strategy among all optimal interventions implemented in the last section following the approach in [1, 2, 14].

4.1. Average cost-effectiveness ratio $(\mathrm{ACER})$

The average cost-effectiveness ratio $(\mathrm{ACER})$ evaluates a single optimal strategy against its baseline option. The $(\mathrm{ACER})$ of a particular optimal strategy is computed by the formula: (19) $\begin{array}{r}\mathrm{ACER}=\frac{\text{Total cost produced by a strategy}}{\text{Total number of infections reduced by the strategy}}\end{array}$(19) where, the total cost produced by the particular strategy is approximated from (20) $\begin{array}{r}{T}_C=\frac{1}{2}{\int }_0^T\left(\sum _{k=1}^4{z}_k{u}_k^2\right)\mathrm{d}t\end{array}$(20) and, the number of infections reduced due the particular strategy is estimated as the difference between the total number of infected individuals without control and the total number of infected individuals with control in simulation period. The strategy with the least $\mathrm{ACER}$ is the most cost effective. The $\mathrm{ACER}$ value for the control strategies is calculated in Table 2 using the parameter values in Table 1. From the Table 2, we conclude that the strategy A is the most cost effective of all for this particular study. Thus, according to $\mathrm{ACER}$ the control strategies from the most to the least are listed as strategy A, strategy B, strategy C, strategy D and strategy E.

Table 2. Total infection averted, total cost and $\mathrm{ACER}$.

Strategy	Infection averted	Total Cost ($)	ACER
Strategy A	67,556.7	1999.9	0.0296
Strategy B	67,883.8	2131.2	0.0314
Strategy C	68,183.9	2234.1	0.0328
Strategy D	69,609.4	2200.8	0.0316
Strategy E	69,799.6	2294.2	0.033

4.2. The incremental cost-effectiveness ratio (ICER)

Next, we consider the incremental cost-effectiveness ratio ($\mathrm{ICER}$) which is used to compare the difference between the costs and health outcomes of any two alternative intervention strategies incrementally that competing for the same limited resources [5, 26, 31]. Thus, $\mathrm{ICER}$ value for two strategies $S_1$ and $S_2$ (say) is calculated by the formula: (21) $\begin{array}{r}\mathrm{ICER}=\frac{\text{Difference in total intervation Costs in Strategies}{S}_1\mathrm{and}{S}_2}{\text{Difference in the total number of infection averted in Strategies}{S}_1\mathrm{and}{S}_2}.\end{array}$(21) The numerator in the formula 21(21) $\begin{array}{r}\mathrm{ICER}=\frac{\text{Difference in total intervation Costs in Strategies}{S}_1\mathrm{and}{S}_2}{\text{Difference in the total number of infection averted in Strategies}{S}_1\mathrm{and}{S}_2}.\end{array}$(21) includes the differences in the cost of interventions, costs of disease averted, costs of averted productivity loses or costs of prevented cases among others. While the denominator determine the differences in health outcomes which includes the total number of infections averted or the number of susceptibility cases prevented [1]. To implement the $\mathrm{ICER}$ the interventions strategies are then ranked according to their increasing order of total number of infection averted as shown in Table 2. The strategy to be discarded from the list of alternative interventions at each step is that corresponding to the highest $\mathrm{ICER}$ value [31]. We first compared the $\mathrm{ICER}$ value for strategies A and B. Their $\mathrm{ICER}$ values are computed as follows: $\begin{array}{r}\mathrm{ICER}(\mathrm{A})=\frac{1999.9-0}{67,556.7-0}=0.0296,\mathrm{ICER}(\mathrm{B})=\frac{2131.2-1999.9}{67,883.8-67,556.7}=\mathrm{0.401.}\end{array}$It follows that, strategy B is more costly and less effective than strategy A. Therefore, it is better to exclude strategy B from the set of alternative interventions and then the strategy A is compared with strategy C by calculating $\mathrm{ICER}$ as follow. $\begin{array}{r}\mathrm{ICER}(\mathrm{A})=\frac{1999.9-0}{67,556.7-0}=0.0296,\mathrm{ICER}(\mathrm{C})=\frac{2234.1-1999.9}{68,183.9-67,556.7}=\mathrm{0.373.}\end{array}$From the values of $\mathrm{ICER}(\mathrm{A})$ and $\mathrm{ICER}(\mathrm{C})$, it can be observed that strategy C is expensive than strategy A. In other word strategy A is less costly and more effective than strategy C. Therefore, strategy C should be removed from the list of alternative interventions and then the strategy A is compared with strategy D. $\begin{array}{r}\mathrm{ICER}(\mathrm{A})=\frac{1999.9-0}{67,556.7-0}=0.0296,\mathrm{ICER}(\mathrm{D})=\frac{2200.8-1999.9}{69,609.4-67,556.7}=\mathrm{0.0979.}\end{array}$Similarly, from the values of $\mathrm{ICER}(\mathrm{A})$ and $\mathrm{ICER}(\mathrm{D})$, it can be observed that strategy A is less costly and more effective than strategy D. Finally, we compared strategy A with strategy E to determine the least cos-effective strategy. $\begin{array}{r}\mathrm{ICER}(\mathrm{A})=\frac{1999.9-0}{67,556.7-0}=0.0296,\mathrm{ICER}(\mathrm{E})=\frac{2294.2-1999.9}{69,799.6-67,556.7}=\mathrm{0.1312.}\end{array}$The result of comparison between the two strategies indicates that strategy E is strongly dominated the strategy A. Thus, the strategy E is more costly and less effective than the strategy A. Consequently, the strategy A is most cost-effective of all the strategies for control of leptospirosis disease infection under consideration for this particular study. This result agrees with the result of the ACER method obtained earlier. Furthermore, this result is confirmed in Figures 7(a) and 7(b).

5. Conclusion

In this paper, we extended the model formulated in [11] into an optimal control problem. In formulation of an optimal control model, we incorporated four time-dependent control function to reduce the spread of the disease in a population. These include personal protection (personal protective equipment (PPE) and personal hygiene), treatment control on infected individuals class, rodenticide on each classes of the rodent, and improvements in slum areas (or draining wet areas). The existence of the optimal control quadruple was verified based on standard results and the Pontryagin's Maximum Principle was used to derive the necessary conditions of the optimality system. The numerical simulations of the optimality system were performed for different control strategies and presented graphically with and without controls. The optimal problem was simulated using the Forward–Backward Sweep method in the Matlab programme. The numerical result revealed that all the implemented control strategies have a significant effect on minimizing the leptospirosis infection in the population, although the combination of all preventive and control measures is most effective. However, our result of cost-effectiveness analysis shows that the rodenticide control only strategy is the most cost-effective strategy among all optimal interventions to reduce the number of infected individuals, rodents and a load of pathogenic in the environment. Thus, the rodenticide control only strategy is recommended to combat the infection of leptospirosis epidemic in the population where there are limited resources.

Figure 7. The total Cost produced and ACER for each Strategy. (a) Total cost Produced by control Strategies. (b) ACER of control Strategies.

Acknowledgements

The authors would like to thank Pan African University for providing the necessary support.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that has been used is supported by previously published papers. The data supporting this work are from previously published articles and they have been cited in this paper.

Additional information

Funding

The author(s) reported there is no funding associated with the work featured in this article.