Cost Effectiveness Analysis of Optimal Rabies Control Strategies

1. Introduction

Rabies is a zoonotic disease that affects Ethiopia’s livestock industry and public health, according to Abera et al. (2024). The burden and distribution of rabies in the country on both humans and animals were researched. By gathering data from the Ethiopian Public Health Institute (EPHI) and the Ministry of Agriculture over five years (2018–2022), the authors carried out an in-depth descriptive analysis of rabies. They showed a higher prevalence and broader spread of rabies throughout the country in humans as well as animals. The results account for the increasing number of suspected cases of human rabies exposure and mortality as animal outbreaks continue to spread. The animal species most frequently impacted were dogs, then camels, cattle, and horses.

Ruan Shigui (2017), Created a simple susceptible, exposed, infectious, and recovered (SEIR) type model to represent the transmission of the rabies virus from dogs to humans as well as between dogs and other dogs. The constructed model was used to reproduce data on human rabies cases in China from 1996 to 2010. Following that, the author included both domestic and stray dogs in the basic model and used it to mimic human rabies data from Guangdong Province, China. Furthermore, from January 2004 to December 2010, the Chinese Ministry of Health released monthly data on human rabies cases. He used this data and created an SEIR model with periodic transmission rates. With the system of ordinary differential equations, Musaili and Chepkwony (2020) updated the SIR model of infectious diseases to take into consideration the transmission of the rabies virus in dogs and included public health education as a control mechanism. They identified the equilibrium points for both endemic and disease-free conditions by computing the reproduction number. An SEIR and SEIV model was presented by Eze et al. (2020) to explain the dynamics of rabies virus transmission in humans and dogs. They computed the disease-free and endemic equilibrium points, as well as the basic reproduction number. They also obtained a control solution for the model, which predicts that increasing pre-exposure prophylaxis in humans and dogs as well as public education help eliminate rabies-related deaths. However, the results indicated that, if complete eradication of the disease is desired, pre- and post-exposure prophylaxis in humans along with vaccination in the dog population. Any combined approach that includes vaccinating the dog population increases the disease’s ability to be eradicated. A mathematical model that represented the dynamics of the rabies virus transmission between dogs and human beings was developed by Thongtha and Modnak (2021). They also used optimal control theory and equilibrium state analysis to try to minimize the cost of rabies vaccinations. A mathematical model was developed by Pantha et al. (2021) to explain the dynamics of rabies transmission in Nepal. Researchers can properly estimate parameters associated with rabies transmission in Nepal, which makes use of annual dog-bite data which was obtained from Nepal over ten years. They discovered that the main factors influencing the number of reproductions are dog-related elements. The findings suggest that jackals might have contributed to the ongoing epidemic of rabies in Nepal, although to a smaller degree than dogs. Additionally, based on the model, control measures significantly decrease the number of cases, though jackal vaccinations are unlikely to be as successful as dog-related preventative measures.

Motivation for the present work

The majority of studies in the literature mainly focused on developing SEIR models for both dog and human populations. In general, SEIR models of dogs assume that the cause of death for all infected dogs is rabies disease. However, medical literature reported that the behavior of rabies in the infected class is not constant but differed linearly. That is, the infected class is divided into three phases: (i) prodromal (ii) furious, and (iii) paralytic [21]. In the first phase i.e., prodromal phase dogs become shy and seek out isolated areas. In the second phase i.e., furious phase dogs become extremely irritable and aggressive. In the end phase or third phase i.e., the paralysis phase dogs lose their capacity to move normally, chew, and swallow.

Here, it is clear that the rrabies-infecteddog’s behavior is not the same always but differs linearly with time. The nature of the rabies disease we intend to include in the current study. In this study, the infected compartment is divided into (i) the prodromal phase and (ii) the furious phase compartments. However, the paralysis phase is omitted from consideration. This is because the dogs in this phase are highly inactive. That is, these paralytic dogs cannot move from one place to another place and therefore they can rarely transmit the disease.

However, it is assumed that rabies only causes mortality among dogs in the furious phase rather than the prodromal phase and the model divided the infection rate into two categories: dog-to-dog transmission rates during the prodromal phase (${\beta }_{dP}$) and the furious phase (${\beta }_{dF}$).

Furthermore, the human population is divided into SEIR compartmental structures. Unlike infected dogs, infected humans do not transmit the disease to others. Hence, in the case of humans splitting of infected class is not meaningful. Hence, the infected human class is left without splitting in this study.

This investigation is motivated by therefore mentioned considerations and is aimed to address these gaps.

In Section 2, the model assumptions are listed. State variables and model parameters are described and presented in tabular forms. A model flowchart is also included. A model consisting of nine compartments is developed: (i) five compartments structure$S_d{E}_d{I}_{dP}{I}_{dF}{R}_d$is for the dog population and (ii) four compartments structure $S_h{E}_h{I}_h{R}_h$is for the human population.

Section 3: defines the model’s boundedness and positivity, as well as equilibrium points of the model are determined. Both the basic reproductive number (R0) and the effective reproductive number (Re) are computed using the next-generation matrix approach. It shows that the disease-free equilibrium (DFE) is globally stable.

Section 4: The sensitivity analysis and Interpretation of Sensitivity Indices presented in the figures show the numerical simulation of the relationship between the sensitive parameters and the effective reproductive number.

Section 5: Pontryagin’s Maximum Principle is applied and the fundamentals of optimal control are examined.

Section 6: the simulation results are presented. The rate of per-exposed and post-exposed vaccination for both populations are considered as control mechanisms.

Section 7: A cost-effective analysis using the Incremental Cost Effectiveness Ratio (ICER) was conducted. Population dynamics are demonstrated through graphs.

The paper ends in Section 8 with recommendations and conclusions.

2. Description and Model Formulation

Here, a system of continuous nonlinear deterministic ordinary differential equations is used to create a mathematical model. The main objective of this model is to evaluate the way per-exposure and post-exposure vaccinations are effective rabies control methods for both dogs and humans. The dog and the human model are the two subpopulations that are involved in the spread of the rabies virus. The dog population is divided into five compartments: (i) susceptible dogs$S_d$ (ii) exposed dogs $E_d$ (iii) infected dogs in the prodromal phase $I_{dP}$ (iv) infected dogs in the furious phase $I_{dF}$ (v) recovered dogs $R_d$ due to pre-exposure vaccination. This $S_d{E}_d{I}_{dP}{I}_{dF}{R}_d$dog model is an extension of the existing SEIR model used to describe the dynamics of dog rabies.

Similarly, the human population is divided into four compartments consisting of (i) susceptible humans $S_h$ (ii) exposed humans$E_h$ (iii) infected humans $I_h$, and (iv) recovered humans $R_h$ due to post-exposure vaccination.

Both susceptible dogs and susceptible human individuals are at risk of catching the disease from rabid dogs. Both exposed dogs and humans have been exposed to the virus but they do not yet show any symptoms. That is, at this stage, they are not capable of transferring the disease.

Infected individuals, both dogs and humans, are unlikely to recover from rabies. If the exposed humans are immunized or given treatment before they become infectious then they will recover from the disease.

The model includes nine ordinary differential equations and is based on the following assumptions:

• Susceptible populations of dogs and humans are recruited at a rate level $A_d and A_h$.
• The rabies infection does not transfer from infected humans to susceptible dogs.
• Rabies transmission from humans to humans was ignored.
• Exposed dogs cannot transfer the rabies disease either to dogs or to humans.
• Giving post-exposed prophylaxis for dogs and humans.
• Giving pre-exposed prophylaxis for dogs and humans.
• For dogs, giving both pre-exposed and post-exposure treatment or prophylaxis.
• The primary way of rabies virus transmission is (i) from infectious dogs to susceptible dogs and (ii) from infectious dogs to susceptible humans.
• Individuals in each compartment have an equal natural death rate.
• All the model parameters are positive quantities.

Based on the aforementioned assumptions, the dynamics of the disease in both dog and human populations are demonstrated through a flowchart and are shown in Figure 1.

2.1. Model Equations

The ordinary differential equation below illustrates the paths of infection based on the transmission model shown in Figure 1.

$$\left.\begin{array}{c} \frac{d{S}_d}{dt}={\mathrm{A}}_d-{\left({{\beta }_{dP}I}_{dP}+{{\beta }_{dF}I}_{dF}\right)S}_d-\left({\mu }_d+v_d\right)S_d {+\alpha }_d{R}_d \\ \frac{d{E}_d}{dt}={({\beta }_{dP}{I}_{dp}+{{\beta }_{dF}I}_{dF})S}_d-\left({\mu }_d+{\delta }_d+{\theta }_d\right)E_d\\ \frac{d{I}_{dp}}{dt}={\delta }_d{E}_d-\left({\mu }_d+{\rho }_d+{\mathit{c}}_{\mathit{d}}\right)I_{dP}\\ \frac{d{I}_{dF}}{dt}={\rho }_d{I}_{dP}-\left({\mu }_d+{\epsilon }_d+{\mathit{c}}_{\mathit{d}}\right)I_{dF}\\ \frac{d{R}_d}{dt}=v_d{S}_d+{\theta }_d{E}_d-\left({{\mu }_d+\alpha }_d\right)R_d \\ \frac{d{S}_h}{dt}={\mathrm{A}}_h-{({{\beta }_{hP}I}_{dP}+{{\beta }_{hF}I}_{dF})S}_h-\left({\mu }_h+v_h\right)S_h{+\alpha }_h{R}_h\\ \frac{d{E}_h}{dt}={({{\beta }_{hP}I}_{dP}+{{\beta }_{hF}I}_{dF})S}_h-\left({{\mu }_h+\delta }_h+{\theta }_h\right)E_h\\ \frac{d{I}_h}{dt}={\delta }_h{E}_h-\left({\mu }_h+{\epsilon }_h\right)I_h \\ \frac{d{R}_h}{dt}=v_h{S}_h+{\theta }_h{E}_h-\left({{\mu }_h+\alpha }_h\right)R_h \end{array}\right\}$$

(1)

Table 1. Description of variables and parameters of the system (1).

Variables and Parameters

Description

$S_d$ $E_d$ $I_{dp}$ $I_{dF}$ $R_d$ $S_h$ $E_h$ $I_h$ $R_h$ $A_d$ ${\beta }_{dP}$ ${\beta }_{dF}$ ${\delta }_d$ ${\rho }_d$ $v_d$ ${\theta }_d$ ${\epsilon }_d$ ${\mu }_d$ ${\alpha }_d$ $A_h$ ${\beta }_{hP}$ ${\beta }_{hF}$ ${\delta }_h$ $v_h$ ${\theta }_h$ ${\epsilon }_h$ ${\mu }_h$ ${\alpha }_h$

Susceptible dog populations Exposed dog population Infectious dog population with prodromal phase Infectious dog population with the furious phase Recovered dog population Susceptible human populations Exposed human populations Infectious human populations Recovered human population Dogs annual crop of newborn puppies Prodromal phase dog -to-dog transmission rate Furious phase dog-to-dog transmission rate The incubation period of dog populations Rate of prodromal to furious stage Pre-exposure prophylaxis (vaccine) for dogs Post-exposure prophylaxis (treatment) for dogs The death rate of dogs due to rabies The natural death rate of dogs Loss of immunity in dogs Human annual birth Prodromal phase dog-to-human transmission rate Furious phase dog -to-Human transmission rate The incubation period of human populations Pre-exposure prophylaxis (vaccine) for humans Post-exposure prophylaxis (treatment) for humans The death rate of humans due to rabies Natural death rate of humans Rate of recovery to Susceptible human

Table 1. Description of variables and parameters of the system (1).

Variables and Parameters

Description

$S_d$ $E_d$ $I_{dp}$ $I_{dF}$ $R_d$ $S_h$ $E_h$ $I_h$ $R_h$ $A_d$ ${\beta }_{dP}$ ${\beta }_{dF}$ ${\delta }_d$ ${\rho }_d$ $v_d$ ${\theta }_d$ ${\epsilon }_d$ ${\mu }_d$ ${\alpha }_d$ $A_h$ ${\beta }_{hP}$ ${\beta }_{hF}$ ${\delta }_h$ $v_h$ ${\theta }_h$ ${\epsilon }_h$ ${\mu }_h$ ${\alpha }_h$

Susceptible dog populations Exposed dog population Infectious dog population with prodromal phase Infectious dog population with the furious phase Recovered dog population Susceptible human populations Exposed human populations Infectious human populations Recovered human population Dogs annual crop of newborn puppies Prodromal phase dog -to-dog transmission rate Furious phase dog-to-dog transmission rate The incubation period of dog populations Rate of prodromal to furious stage Pre-exposure prophylaxis (vaccine) for dogs Post-exposure prophylaxis (treatment) for dogs The death rate of dogs due to rabies The natural death rate of dogs Loss of immunity in dogs Human annual birth Prodromal phase dog-to-human transmission rate Furious phase dog -to-Human transmission rate The incubation period of human populations Pre-exposure prophylaxis (vaccine) for humans Post-exposure prophylaxis (treatment) for humans The death rate of humans due to rabies Natural death rate of humans Rate of recovery to Susceptible human

3. Model Analysis

3.1. Positivity

In this subsection, we show that the solutions of the system of model equations (1) are positive. This is stated in the form of a theorem accompanied by proofs

Positivity of dog population

Theorem 1:

Every solution of the system of model equations representing dog populations given in (1) together with the initial conditions exists in the interval$[\begin{array}{cc}0,& \mathrm{\infty })\end{array}$. Also, the solutions are positive.

Proof: 

Here we show that the solutions for the dog population equations given in model system (1) exist Also, we show that they are non-negative i.e., $S_d\left(t\right)>0, E_d\left(t\right)\ge 0, I_{dP}\left(t\right)\ge 0, I_{dF}\left(t\right)\ge 0, R_d\left(t\right)\ge 0$. That is, if the initial conditions $S_d\left(0\right), E_d\left(0\right), I_{dP}\left(0\right), I_{dF}\left(0\right), R_d\left(0\right)$ are non-negative then so are the variables $S_d\left(t\right), E_d\left(t\right), I_{dP}\left(t\right), I_{dF}\left(t\right), R_d\left(t\right)$ for all $t\ge 0$.

The solutions $\left\{S_d\left(t\right), E_d\left(t\right), I_{dP}\left(t\right), I_{dF}\left(t\right), R_d\left(t\right)\right\}$ of the system (1) together with the initial conditions exits and unique on $\left[\begin{array}{cc}0,& k\end{array}\right)$ where $0<k<+\infty$since the right-hand side of the system is completely continuous and locally Lipshitzian on $C$.

We now consider the dog population equations given in (1), one by one, and show that the solutions of the dog variables are non-negative i.e., $S_d\left(t\right)>0, E_d\left(t\right)\ge 0, I_{dP}\left(t\right)\ge 0, I_{dF}\left(t\right)\ge 0, R_d\left(t\right)\ge 0$ for all $t\ge 0$.

Positivity of susceptible dog population: Consider the model equation for susceptible dog population$\mathrm{d}{\mathrm{S}}_{\mathrm{d}}/\mathrm{d}\mathrm{t}={\mathrm{A}}_{\mathrm{d}}-{\left({\mathsf{\beta }}_{\mathrm{d}\mathrm{P}}{\mathrm{I}}_{\mathrm{d}\mathrm{p}}+{\mathsf{\beta }}_{\mathrm{d}\mathrm{F}}{\mathrm{I}}_{\mathrm{d}\mathrm{F}}\right)\mathrm{S}}_{\mathrm{d}}-\left({\mathsf{\mu }}_{\mathrm{d}}+{\mathrm{v}}_{\mathrm{d}}\right){\mathrm{S}}_{\mathrm{d}}{+\mathsf{\alpha }}_{\mathrm{d}}{\mathrm{R}}_{\mathrm{d}}$. It can be expressed without loss of generality, after elimination of the positive term $A_d$which appears on the right-hand side, as an inequality as $d{S}_d/dt\ge -{\left({\beta }_{dP}{I}_{dp}+{\beta }_{dF}{I}_{dF}\right)S}_d-\left({\mu }_d+v_d\right)S_d$.This inequality can also be expressed as$d{S}_d/S_d\ge {\pi }_{1}dt.$ Here, the function ${\pi }_1$denoting the expression ${\pi }_1=-\left({\beta }_{dP}{I}_{dp}+{\beta }_{dF}{I}_{dF}\right)-\left({\mu }_d+v_d\right)$can be negative zero or positive valued. Now, using the variables separation method and upon integrating, the solution of $d{S}_d/S_d\ge {\pi }_{1}dt$can be obtained as $S_d\left(t\right)\ge S_d\left(0\right) e^{\int {\pi }_{1}dt}.$Here, the integral constant $S_d\left(0\right)$is the initial susceptible dog population and is assumed to be a non-negative quantity $S_d\left(0\right)\ge 0$. Similarly, an exponential function is always a non-negative quantity irrespective of the value of the exponent i.e., $e^{\int {\pi }_{1}dt}\ge 0$. Hence, we conclude that$S_d\left(t\right)\ge 0$. That is, the susceptible dog population size is a positive quantity.

Positivity of exposed dog population: Consider the model equation for exposed dog population$d{E}_d/dt\ge -\left({\mu }_d+{\delta }_d+{\theta }_d\right)E_d$. Alternately, it can be expressed as $d{E}_d/E_d\ge {\pi }_{2}dt$. Here, the function ${\pi }_2$ denoting the expression ${\pi }_2=-\left({\mu }_d+{\delta }_d+{\theta }_d\right)$can be negative zero or positive valued. Now, using the variables separation method and applying integration, the solution of$d{E}_d/E_d\ge {\pi }_{2}dt$can be obtained as $E_d\left(t\right)=E_d\left(0\right) e^{\int {\pi }_{2}dt}.$ Here, the integral constant $E_d\left(0\right)$is the initial population of exposed dogs and is assumed to be a non-negative quantity $E_d\left(0\right)\ge 0$. Similarly, an exponential function is always a non-negative quantity irrespective of the value of the exponent i.e., $e^{\int {\pi }_{2}dt}\ge 0$. Hence, we conclude that $E_d\left(t\right)\ge 0$. That is, the exposed dog population size is a positive quantity.

Positivity of Infectious dog population in prodromal phase: Consider the model equation for prodromal dog population$d{I}_{dp}/dt={\delta }_d{E}_d-\left({\mu }_d+{\rho }_d+{\mathit{c}}_{\mathit{d}}\right)I_{dP}$. Here, the term${\delta }_d{E}_d$is a positive quantity since${ \delta }_d$is a positive parameter and the exposed dog population$E_d$ is already shown positive. Now, the aforementioned ODE can be expressed without loss of generality, after eliminating the positive term${\delta }_d{E}_d$, as an inequality as$d{I}_{dp}/dt\ge -\left({\mu }_d+{\rho }_d+{\mathit{c}}_{\mathit{d}}\right)I_{dP}$. Now, using the variables separation method and integrating, the solution can be obtained as $I_{dp}\left(t\right)=I_{dp}\left(0\right) e^{-\int \left({\mu }_d+{\rho }_d+{\mathit{c}}_{\mathit{d}}\right) dt}.$ Here, the integral constant $I_{dp}\left(0\right)$is the initial population of prodromal infectious dogs and is assumed to be a non-negative quantity $I_{dp}\left(0\right)\ge 0$. Similarly, the exponential term is also a positive quantity i.e., $e^{-\int \left({\mu }_d+{\rho }_d+{\mathit{c}}_{\mathit{d}}\right) dt} \ge 0$since the integrand in the exponent consists of positive parameters only. Hence, we conclude that$I_{dp}\left(t\right)\ge 0$. That is, the prodromal infectious dog population size is a positive quantity.

Positivity of Infectious dog population in furious phase: Consider the model equation for infectious dog population in furious phase$d{I}_{dF}/dt={\rho }_d{I}_{dP}-\left({\mu }_d+{\epsilon }_d+{\mathit{c}}_{\mathit{d}}\right)I_{dF}$. Here, the term ${\rho }_d{I}_{dP}$is positive since $\rho d$is a positive parameter and the prodromal infectious dog population$I_{dP}$ is already shown positive. Now, the aforementioned ODE can be expressed without loss of generality, after eliminating the positive term${\rho }_d{I}_{dP}$, as an inequality as$d{I}_{dF}/dt\ge -\left({\mu }_d+{\epsilon }_d+{\mathit{c}}_{\mathit{d}}\right)I_{dF}$. Now, using the variables separation method and integrating, the solution can be obtained as $I_{dF}\left(t\right)\ge I_{dF}\left(0\right) e^{-\int \left({\mu }_d+{\epsilon }_d+{\mathit{c}}_{\mathit{d}}\right) dt}.$ Here, the integral constant $I_{dF}\left(0\right)$is the initial population of furious infectious dogs and is assumed to be a non-negative quantity $I_{dF}\left(0\right)\ge 0$. Similarly, the exponential term is also positive i.e., $e^{-\int \left({\mu }_d+{\epsilon }_d+{\mathit{c}}_{\mathit{d}}\right) dt} \ge 0$since the integrand in the exponent consists of positive parameters only. Hence, we conclude that$I_{dF}\left(t\right)\ge 0$. That is, the furious infectious dog population size is a positive quantity.

Positivity of Recovered Dog Population: Consider the model equation for the recovered dog population as$d{R}_d/dt=v_d{S}_d+{\theta }_d{E}_d-\left({{\mu }_d+\alpha }_d\right)R_d$. Here, the term$s v_d{S}_d$ and ${\theta }_d{E}_d$ are positive since ${\theta }_d$and $v_d$ are positive parameters and the susceptible and Exposed dog population$s$are already shown positive. Now, the aforementioned ODE can be expressed without loss of generality, after eliminating the positive terms$v_d{S}_d$, ${\theta }_d{E}_d,$as an inequality as$d{R}_d/dt\ge -\left({{\mu }_d+\alpha }_d\right)R_d$. Now, using the variables separation method and integrating, the solution can be obtained as $R_d \left(t\right)\ge R_d \left(0\right) e^{-\int \left({{\mu }_d+\alpha }_d\right) dt}.$ Here, the integral constant $R_d\left(0\right)$is the initial population of recovered dogs and is assumed to be a non-negative quantity $R_d\left(0\right)\ge 0$. Similarly, the exponential term is also positive i.e., $e^{-\int \left({{\mu }_d+\alpha }_d\right) dt} \ge 0$since the integrand in the exponent consists of positive parameters only. Hence, we conclude that$R_d\left(t\right)\ge 0$. That is, the recovered dog population size is a positive quantity.

Positivity of the human population

Theorem 2:

Every solution of the system of model equations representing human populations given in (1) together with the initial conditions exists in the interval$[\begin{array}{cc}0,& \mathrm{\infty })\end{array}$. Also, the solutions are positive.

Proof: 

Here we show that the solutions for the human population equations given in model system (1) exist. Also, we show that they are non-negative i.e., $S_h\left(t\right)\ge 0, E_h\left(t\right)\ge 0, I_h\left(t\right)\ge 0, R_h\left(t\right)\ge 0$. That is, If the initial conditions $S_h\left(0\right), E_h\left(0\right), I_h\left(0\right), and R_h\left(0\right)$ are non-negative then so are the variables $S_h\left(t\right), E_h\left(t\right), I_h\left(t\right), R_h\left(t\right)$ for all $t\ge 0$.

The solutions $\left\{S_h\left(t\right), E_h\left(t\right), I_h\left(t\right), R_h\left(t\right)\right\}$ of the system (1) together with the initial conditions exits and unique on $\left[\begin{array}{cc}0,& k\end{array}\right)$ where $0<k<+\infty$since the right-hand side of the system is completely continuous and locally Lipshitzian on $C$.

We now consider the human population equations given in (1), one by one, and show that the solutions of the human variables are non-negative i.e., $S_h\left(t\right)\ge 0, E_h\left(t\right)\ge 0, I_h\left(t\right)\ge 0, R_h\left(t\right)\ge 0$ for all $t\ge 0$.

Positivity of susceptible human population: Consider the model equation for susceptible human population$d{S}_h/dt={ A}_h-{\beta }_h{\left(I_{dp}+I_{dF}\right)S}_h-{({\mu }_h+v_h)S}_h{+\alpha }_h{R}_h$. It can be expressed without loss of generality, after eliminating the positive term $A_h$appearing on the right-hand side, as an inequality as $d{S}_h/dt\ge -{\left({\beta }_{hP}{I}_{dp}+{{\beta }_{hF}I}_{dF}\right)S}_h-{({\mu }_h+v_h)S}_h$. This inequality can also be expressed as$d{S}_h/S_h\ge {\mathrm{Q}}_{1}dt.$ Here, the function ${\mathrm{Q}}_1$denoting the expression ${\mathrm{Q}}_1=-\left({\beta }_{hP}{I}_{dp}+{{\beta }_{hF}I}_{dF}\right)-({\mu }_h+v_h)$can be negative or zero or positive valued. Now, using the variables separation method and upon integrating, the solution can be obtained as $S_h\left(t\right)\ge S_h\left(0\right) e^{\int Q_{1}dt}.$Here, the integral constant $S_h\left(0\right)$is the initial susceptible human population and is assumed to be a non-negative quantity $S_h\left(0\right)\ge 0$. Similarly, an exponential function is always a non-negative quantity irrespective of the value of the exponent i.e., $e^{\int Q_{1}dt}\ge 0$. Hence, we conclude that$S_h\left(t\right)\ge 0$. That is, the susceptible human population size is a positive quantity.

Positivity of exposed human population: Consider the model equation for exposed human population$d{E}_h/dt\ge -\left({{\mu }_h+\delta }_h+{\theta }_h\right)E_h$.

Alternately, it can be expressed as$d{E}_h/E_h\ge {\mathrm{Q}}_{2}dt$. Here, the function ${\mathrm{Q}}_2$denoting the expression ${\mathrm{Q}}_2=-\left({{\mu }_h+\delta }_h+{\theta }_h\right)$can be negative zero or positive valued. Now, using the variables separation method and applying integration, the solution can be obtained as $E_h\left(t\right)=E_h\left(0\right) e^{\int Q_{2}dt}.$ Here, the integral constant $E_h\left(0\right)$is the initial population of exposed humans and is assumed to be a non-negative quantity $E_h\left(0\right)\ge 0$. Similarly, an exponential function is always a non-negative quantity irrespective of the value of the exponent i.e., $e^{\int Q_{2}dt}\ge 0$. Hence, we conclude that$E_h\left(t\right)\ge 0$. That is, the exposed human population size is a positive quantity.

Positivity of Infectious human population: Consider the model equation for the infectious human population$d{I}_h/dt={\delta }_h{E}_h-\left({\mu }_h+{\epsilon }_h\right)I_h$. Here, the term ${\delta }_h{E}_h$is a positive quantity since ${\delta }_h$is a positive parameter and the exposed human population$E_h$ is already shown positive. Now, the aforementioned ODE can be expressed without loss of generality, after eliminating the positive term${\delta }_h{E}_h$, as an inequality as$d{I}_h/dt\ge -\left({\mu }_h+{\epsilon }_h\right)I_h$. Now, using the variables separation method and integrating, the solution can be obtained as $I_h\left(t\right)\ge I_h\left(0\right) e^{-\int \left({\mu }_h+{\epsilon }_h\right) dt}.$ Here, the integral constant $I_h\left(0\right)$is the initial population of infectious humans and is assumed to be a non-negative quantity$I_h\left(0\right)\ge 0$. Similarly, the exponential term is also a positive quantity i.e., $e^{-\int \left({\mu }_h+{\epsilon }_h\right) dt} \ge 0$since the integrand in the exponent consists of positive parameters only. Hence, we conclude that$I_h\left(t\right)\ge 0$. That is, the infectious human population size is a positive quantity.

Positivity of Recovered human population: Consider the model equation for the recovered dog population as$d{R}_h/dt=v_h{S}_h+{\theta }_h{E}_h-\left({{\mu }_h+\alpha }_h\right)R_h$. Here, the terms $v_h{S}_h$ and${\theta }_h{E}_h$are positive since $v_h, and \theta h$are positive parameters and the susceptible and exposed human population are already shown positive. Now, the aforementioned ODE can be expressed without loss of generality, after eliminating the positive terms $v_h{S}_h, {\theta }_h{E}_h$ as an inequality as$d{R}_h/dt\ge -\left({{\mu }_h+\alpha }_h\right)R_h$. Now, using the variables separation method and integrating, the solution can be obtained as $R_h \left(t\right)\ge R_h \left(0\right) e^{-\int \left({{\mu }_h+\alpha }_h\right) dt}.$ Here, the integral constant $R_h\left(0\right)$is the initial population of recovered humans and is assumed to be a non-negative quantity i.e., $R_h\left(0\right)\ge 0$. Similarly, the exponential term is also positive i.e., $e^{-\int \left({{\mu }_h+\alpha }_h\right) dt} \ge 0$since the integrand in the exponent consists of positive parameters only. Hence, we conclude that$R_h\left(t\right)\ge 0$. That is, the recovered human population size is a positive quantity.

Theorem 3:

The feasible region$\mathsf{\Omega }$defined by$\mathsf{\Omega }={\mathsf{\Omega }}_{\mathrm{d}}\mathrm{x}{\mathsf{\Omega }}_{\mathrm{h}}\subset {\mathrm{R}}_{+}^5\mathrm{x}{\mathrm{R}}_{+}^4$is bounded.

Here,

$$\begin{gathered} {\mathsf{\Omega }}_d=\left\{\left(\begin{array}{ccccc}{S}_d\left(t\right),& E_d\left(t\right),& I_{dP}\left(t\right),& I_{dF}\left(t\right),& R_d\left(t\right)\end{array}\right)\in R_{+}^5 : {\mathrm{N}}_{\mathrm{d}}\left(t\right)\le \frac{A_d}{{\mathsf{\mu }}_{\mathrm{d}}} \right\}{\mathsf{\Omega }}_h \\ =\left\{\left(\begin{array}{cccc}{S}_h\left(t\right),& E_h\left(t\right),& I_h\left(t\right),& R_h\left(t\right)\end{array}\right)\in R_{+}^4 : {\mathrm{N}}_{\mathrm{h}}\left(t\right)\le \frac{A_h}{{\mathsf{\mu }}_{\mathrm{h}}} \right\} \end{gathered}$$

Also, $N_d\left(t\right)=S_d\left(t\right)+E_d\left(t\right)+I_{dP}\left(t\right)+I_{dF}\left(t\right)+R_d\left(t\right)$is the total dog population.

Similarly,$N_h\left(t\right)=S_h\left(t\right)+E_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right)$is the total human population.

Furthermore, the sets $\mathsf{\Omega },{\mathsf{\Omega }}_d,$and ${\mathsf{\Omega }}_h$ are all real-valued.

Proof: 

Here, we show that the model population is bounded.

That is, the dog population is bounded i.e., ${\mathrm{N}}_{\mathrm{d}}\left(t\right)\le \frac{A_d}{{\mathsf{\mu }}_{\mathrm{d}}}$

Also, the human population is bounded i.e., ${\mathrm{N}}_{\mathrm{h}}\left(t\right)\le \frac{A_h}{{\mathsf{\mu }}_{\mathrm{h}}}$

3.2. Boundedness of the Model

The dog population is bounded

We now show that if the total dog population size is given by $N_d\left(t\right)=S_d\left(t\right)+E_d\left(t\right)+I_{dP}\left(t\right)+I_{dF}\left(t\right)+R_d\left(t\right)$ then$\underset{t\to \infty }{\lim }{\mathrm{N}}_{\mathrm{d}}\left(t\right)\le \frac{A_d}{{\mathsf{\mu }}_{\mathrm{d}}}$. In other words, the total size of the dog population given in model system (1) is bounded above.

Consider that $N_d$ denotes the total dog population at any time t.

Therefore

$$N_d=S_d+E_d+I_{dP}+I_{dF}+R_d$$

(3.1)

Upon derivation of both sides of equation (2.1) concerning time, we obtain

$$\frac{\mathrm{d}{\mathrm{N}}_{\mathrm{d}}}{\mathrm{d}\mathrm{t}}=\frac{\mathrm{d}{\mathrm{S}}_{\mathrm{d}}}{\mathrm{d}\mathrm{t}}+\frac{\mathrm{d}{\mathrm{E}}_{\mathrm{d}}}{\mathrm{d}\mathrm{t}}+\frac{\mathrm{d}{\mathrm{I}}_{\mathrm{d}\mathrm{P}}}{\mathrm{d}\mathrm{t}}+\frac{\mathrm{d}{\mathrm{I}}_{\mathrm{d}\mathrm{F}}}{\mathrm{d}\mathrm{t}}+\frac{\mathrm{d}{\mathrm{R}}_{\mathrm{d}}}{\mathrm{d}\mathrm{t}}$$

(3.2)

Now, making use of the system (1) and substituting for the differential terms appearing on the right-hand side of (3.2), the equation reduces to the form as

$${\frac{\mathrm{d}{\mathrm{N}}_{\mathrm{d}}}{\mathrm{d}\mathrm{t}}=\mathrm{A}}_{\mathrm{d}}-{\mathsf{\mu }}_{\mathrm{d}}{\mathrm{N}}_{\mathrm{d}}-{\mathrm{c}}_{\mathrm{d}}{\mathrm{I}}_{\mathrm{d}\mathrm{P}}-({\mathrm{c}}_{\mathrm{d}}+{\mathsf{\epsilon }}_{\mathrm{d}}){\mathrm{I}}_{\mathrm{d}\mathrm{F}}$$

(3.3)

Here, the terms ${\mathrm{c}}_{\mathrm{d}}{\mathrm{I}}_{\mathrm{d}\mathrm{P}}$and $\left({\mathrm{c}}_{\mathrm{d}}+{\mathsf{\epsilon }}_{\mathrm{d}}\right){\mathrm{I}}_{\mathrm{d}\mathrm{F}}$are positive since all the model parameters are assumed to be positive and all the model variables are proved to be positive. Hence, upon removing these two positive terms, the equation (3.3) can be expressed as an inequality

$${\frac{\mathrm{d}{\mathrm{N}}_{\mathrm{d}}}{\mathrm{d}\mathrm{t}}+{\mathsf{\mu }}_{\mathrm{d}}{\mathrm{N}}_{\mathrm{d}}\le \mathrm{A}}_{\mathrm{d}}$$

It is a first-order nonlinear ODE with constant coefficients and its solution is given by

$${\mathrm{N}}_{\mathrm{d}}\le \frac{A_d}{{\mathsf{\mu }}_{\mathrm{d}}}+\left[{\mathrm{N}}_{\mathrm{d}}\left(0\right)-\frac{A_d}{{\mathsf{\mu }}_{\mathrm{d}}}\right]e^{{-\mathsf{\mu }}_{\mathrm{d}}t}$$

(3.4)

As $\mathrm{t}\to \mathrm{\infty }$ in equation (3.4) the population size ${\mathrm{N}}_{\mathrm{d}}\to \frac{{\mathrm{A}}_d}{{\mathsf{\mu }}_{\mathrm{d}}}$ which implies that $0{\le \mathrm{N}}_{\mathrm{d}}\le \frac{{\mathrm{A}}_d}{{\mathsf{\mu }}_{\mathrm{d}}}$, Thus the feasible solution set of the system equation of the model enters and remains in the region:

$${\mathsf{\Omega }}_d=\left\{S_d,E_d,I_{dP},{I_{dF,} R}_d:{\mathrm{N}}_{\mathrm{d}}\le \frac{{\mathrm{A}}_d}{{\mathsf{\mu }}_{\mathrm{d}}}\right\}\in R_{+}^5$$

(3.5)

The human population is bounded

We now show that if the total human population size is given by $N_h\left(t\right)=S_h\left(t\right)+E_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right)$ then$\underset{t\to \infty }{\lim }{\mathrm{N}}_{\mathrm{h}}\left(t\right)\le \frac{A_h}{{\mathsf{\mu }}_{\mathrm{h}}}$. In other words, the total size of the human population given in model system (1) is bounded above. Consider that $N_h$ denotes the total human population at any time t.

Therefore,

$$N_h=S_h+E_h+I_h+R_h$$

(3.6)

Upon derivation of both sides of equation (3.6) concerning time, we obtain

$$\frac{\mathrm{d}{\mathrm{N}}_{\mathrm{h}}}{\mathrm{d}\mathrm{t}}=\frac{\mathrm{d}{\mathrm{S}}_{\mathrm{h}}}{\mathrm{d}\mathrm{t}}+\frac{\mathrm{d}{\mathrm{E}}_{\mathrm{h}}}{\mathrm{d}\mathrm{t}}+\frac{\mathrm{d}{\mathrm{I}}_{\mathrm{h}}}{\mathrm{d}\mathrm{t}}+\frac{\mathrm{d}{\mathrm{R}}_{\mathrm{h}}}{\mathrm{d}\mathrm{t}}$$

(3.7)

Now, making use of the system (1) and substituting for the differential terms appearing on the right-hand side of (3.7), the equation reduces to the form after algebraic simplifications as

$${\frac{\mathrm{d}{\mathrm{N}}_{\mathrm{h}}}{\mathrm{d}\mathrm{t}}=\mathrm{A}}_{\mathrm{h}}-{\mathsf{\mu }}_{\mathrm{h}}{\mathrm{N}}_{\mathrm{h}}-{\mathsf{\epsilon }}_{\mathrm{h}}{\mathrm{I}}_{\mathrm{h}}$$

(3.8)

Here, the term ${\mathsf{\epsilon }}_{\mathrm{h}}{\mathrm{I}}_{\mathrm{h}}$is positive since all the model parameters are assumed to be positive and all the model variables are proved to be positive.

Hence, upon removing this negative term, the equation (3.8) can be expressed as an inequality as

$${\frac{\mathrm{d}{\mathrm{N}}_{\mathrm{h}}}{\mathrm{d}\mathrm{t}}+{\mathsf{\mu }}_{\mathrm{h}}{\mathrm{N}}_{\mathrm{h}}\le \mathrm{A}}_{\mathrm{h}}$$

It is a first-order nonlinear ODE with constant coefficients and its solution is given by

$${\mathrm{N}}_{\mathrm{d}}\le \frac{A_h}{{\mathsf{\mu }}_{\mathrm{h}}}+\left[{\mathrm{N}}_{\mathrm{h}}\left(0\right)-\frac{A_h}{{\mathsf{\mu }}_{\mathrm{h}}}\right]e^{{-\mathsf{\mu }}_{\mathrm{h}}t}$$

(3.9)

As $\mathrm{t}\to \mathrm{\infty }$ in equation (3.9) the population size ${\mathrm{N}}_{\mathrm{h}}\to \frac{{\mathrm{A}}_h}{{\mathsf{\mu }}_{\mathrm{h}}}$ which implies that $0{\le \mathrm{N}}_{\mathrm{h}}\le \frac{{\mathrm{A}}_h}{{\mathsf{\mu }}_{\mathrm{h}}}$, Thus the feasible solution set of the system equation of the model enters and remains in the region:

$${\mathsf{\Omega }}_h=\left\{S_h,E_h,I_h,{ R}_h:{\mathrm{N}}_{\mathrm{h}}\le \frac{{\mathrm{A}}_h}{{\mathsf{\mu }}_{\mathrm{h}}}\right\}\in R_{+}^4$$

(3.10)

Therefore, from equations (3.5) and (3.10) the region $\mathsf{\Omega }={\mathsf{\Omega }}_{+}^5\times {\mathsf{\Omega }}_{+}^4$ is positively invariant and the model (1) is well-posed or biologically and epidemiologically. The proof of Theorem 3. is complete.

3.4. Disease-Free Equilibrium ${\mathcal{E}}_0$

To determine the disease-free equilibrium,${\mathcal{E}}_0$ for the basic model we consider a scenario where none of the compartments are affected by rabies infection, i.e., when $E_d=I_{dP}=I_{dF}=E_h=I_h=0$, and the control parameters $v_d,$ ${\theta }_d,$ $v_h,$ ${\theta }_h$ are zero. We set the right-hand side of system (1) to zero, leading to the following outcome:

$${\mathrm{A}}_d-{{\mu }_{d}S}_d=0$$

$${\mathrm{A}}_h-{\mu }_h{S}_h=0$$

$${\mathcal{E}}_0={(\mathit{s}}_{\mathit{d}}^0,{\mathit{s}}_{\mathit{h}}^0)=\left(\frac{{\mathrm{A}}_d}{{\mu }_d},\frac{{\mathrm{A}}_d}{{\mu }_d}\right)$$

(3.11)

To determine the disease-free equilibrium for the effective model is similar to the basic model but considering the parameters $v_d,$ ${\theta }_d,$ $v_h,$ ${\theta }_h$ are different from zero

$$\left.\begin{array}{c} {\mathrm{A}}_d-{\left({{\beta }_{dP}I}_{dP}+{{\beta }_{dF}I}_{dF}\right)S}_d-\left({\mu }_d+v_d\right)S_d {+\alpha }_d{R}_d=0\\ v_d{S}_d+{\theta }_d{E}_d-\left({{\mu }_d+\alpha }_d\right)R_d=0 \\ {\mathrm{A}}_h-{\left({{\beta }_{hP}I}_{dP}+{{\beta }_{hF}I}_{dF}\right)S}_h-\left({\mu }_h+v_h\right)S_h{+\alpha }_h{R}_h=0 \\ v_h{S}_h+{\theta }_h{E}_h-\left({{\mu }_h+\alpha }_h\right)R_h=0 \end{array}\right\}$$

(3.12)

One (3.12) is algebraically manipulated, the disease-free-equilibrium point can be represented as

$${\mathcal{E}}_e=\left({S_d^e, E}_d^e, I_{dP}^e, I_{dF}^e, R_d^e, S_h^e, E_h^e, I_h^e, R_h^e\right).$$

When we manipulate equation (3.12) algebraically, we get the disease-free equilibrium point for an effective model

$${\mathcal{E}}_e=\left(\frac{A_d\left({{\mu }_d+\alpha }_d\right)}{{\mu }_d\left({{\mu }_d+\alpha }_d+v_d\right)},0,0,0,\frac{A_d{v}_d}{{\mu }_d\left({{\mu }_d+\alpha }_d+v_d\right)},\frac{A_h\left({{\mu }_h+\alpha }_h\right)}{{\mu }_h\left({{\mu }_h+\alpha }_h+v_h\right)},0,0,\frac{A_h{v}_h}{{\mu }_h\left({{\mu }_h+\alpha }_h+v_h\right)}\right)$$

(3.13)

3.5. Basic Reproduction Number $\left({\mathrm{R}}_0\right)$

To determine the basic reproduction number ($R_0$.), we apply the next-generation matrix method [7]. We can obtain the matrices $f_i$ and $V_i$for the new infection terms and the remaining transfer terms, taking into consideration that our infected compartments are $E_d,I_{dP},I_{dF}, E_h,$ and${ I}_h$

$$\begin{array}{c} \\ \frac{d{E}_d}{dt}={({\beta }_{dP}{I}_{dp}+{{\beta }_{dF}I}_{dF})S}_d-\left({\mu }_d+{\delta }_d\right)E_d\\ \frac{d{I}_{dp}}{dt}={\delta }_d{E}_d-\left({\mu }_d+{\rho }_d\right)I_{dP}\\ \frac{d{I}_{dF}}{dt}={\rho }_d{I}_{dP}-\left({\mu }_d+{\epsilon }_d\right)I_{dF} \\ \frac{d{E}_h}{dt}={\left({{\beta }_{hP}I}_{dP}+{{\beta }_{hF}I}_{dF}\right)S}_h-\left({{\mu }_h+\delta }_h\right)E_h \\ \frac{d{I}_h}{dt}={\delta }_h{E}_h-\left({\mu }_h+{\epsilon }_h\right)I_h \end{array}$$

(3.14)

From equation (3.14) take the vectors $f_i$ and $V_i$

$$f_i=\left(\begin{array}{c}{\beta }_d{(I_{dp}+I_{dF})S}_d\\ 0\\ \begin{array}{c}0\\ {\beta }_h{\left(I_{dp}+I_{dF}\right)S}_h\\ 0\end{array}\end{array}\right)\& V_i=\left(\begin{array}{c}\left({\mu }_d+{\delta }_d\right)E_d\\ \left({\mu }_d+{\rho }_d\right)I_{dP}-{\delta }_d{E}_d\\ \begin{array}{c}\left({\mu }_d+{\epsilon }_d\right)I_{dF}-{\rho }_d{I}_{dP}\\ \begin{array}{c}\left({{\mu }_h+\delta }_h\right)E_h\\ \left({\mu }_h+{\epsilon }_h\right)I_h-{\delta }_h{E}_h\end{array}\end{array}\end{array}\right)$$

The Jacobian matrices $f_i$and$V_i$at disease-free equilibrium from equation (3.11), respectively, are given by

$$F=\left(\begin{array}{ccc}\begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}\frac{{ {\beta }_{dP}A}_d}{{\mu }_d}\\ 0\\ \begin{array}{c}0\\ \frac{{{\beta }_{hP }A}_h}{{\mu }_h}\\ 0\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\frac{ {\beta }_{dF}{A}_d}{{\mu }_d}\\ 0\\ \begin{array}{c}0\\ \frac{{ {\beta }_{hF}A}_h}{{\mu }_h}\\ 0\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}\end{array}\end{array}\right)$$

(3.15)

$$V=\left(\begin{array}{ccc}\begin{array}{c}({\mu }_d+{\delta }_d)\\ -{\delta }_d\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ \left({\mu }_d+{\rho }_d\right)\\ \begin{array}{c}-{\rho }_d\\ 0\\ 0\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}0\\ 0\\ \begin{array}{c}\left({\mu }_d+{\epsilon }_d\right)\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ \left({{\mu }_h+\delta }_h\right)\\ -{\delta }_h\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ \left({\mu }_h+{\epsilon }_h\right)\end{array}\end{array}\end{array}\end{array}\right)$$

Given matrix V, its inverse is given by

$$V^{-1}=\left(\begin{array}{ccccc}\frac{1}{({\mu }_d+{\delta }_d)}& 0& 0& 0& 0\\ \frac{{\delta }_d}{({\mu }_d+{\delta }_d)\left({\mu }_d+{\rho }_d\right)}& \frac{1}{\left({\mu }_d+{\rho }_d\right)}& 0& 0& 0\\ \frac{{\rho }_d{\delta }_d}{\left({\mu }_d+{\rho }_d\right)({\mu }_d+{\delta }_d)\left({\mu }_d+{\epsilon }_d\right)}& \frac{{\rho }_d}{\left({\mu }_d+{\rho }_d\right)\left({\mu }_d+{\epsilon }_d\right)} & \frac{1}{\left({\mu }_d+{\epsilon }_d\right)}& 0& 0\\ 0& 0& 0& \frac{1}{\left({{\mu }_h+\delta }_h\right)}& 0\\ 0& 0& 0& \frac{{\delta }_h}{\left({{\mu }_h+\delta }_h\right)\left({\mu }_d+{\epsilon }_d\right)}& \frac{1}{\left({\mu }_d+{\epsilon }_d\right)}\end{array}\right)$$

(3.16)

Take F and${ V}^{-1}$ from equations (3.15) and (3.16), the symbol $\mathsf{\rho }\left(\mathrm{F}{\mathrm{V}}^{-1}\right)$ represents the dominant eigenvalue in magnitude, or spectral radius, of a matrix $\mathrm{F}{\mathrm{V}}^{-1}$, thus basic reproductive number defended by $R_0=\rho \left(F{V}^{-1}\right)$

As a result, the following gives the effective reproduction number:

$${\mathbf{R}}_0=\frac{{\mathit{A}}_{\mathit{d}}{\mathsf{\delta }}_{\mathbf{d}}({\mathsf{\beta }}_{\mathbf{d}\mathbf{P}}\left({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\epsilon }}_{\mathbf{d}}\right)+{\mathsf{\rho }}_{\mathbf{d}}{\mathsf{\beta }}_{\mathbf{d}\mathbf{F}})}{{\mathsf{\mu }}_{\mathbf{d}}\left({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\rho }}_{\mathbf{d}}\right)\left({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\epsilon }}_{\mathbf{d}}\right)({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\delta }}_{\mathbf{d}})}$$

(3.17)

3.5.1. Effective Reproduction Number

The average number of diseases caused by a single infectious individual in a society when intervention techniques are implemented—in this case, vaccinations for humans and dogs both before and after exposure—is known as the effective reproduction number.

By using the same method as ${\mathrm{R}}_0$ but taking into consideration vaccines for dogs and humans both before and after exposure, and culling of infected dogs i.e.,${\theta }_d, {\theta }_h, v_d,v_h c_d\ne 0$, the effective reproduction number $R_e$ of system (1) is calculated.$R_e=\rho \left(F{V}^{-1}\right)$ is the spectral radius, or dominant eigenvalue, of $F{V}^{-1}.$

The Jacobian matrices F and V at effective disease-free equilibrium ${\mathcal{E}}_e$, respectively, are given by

$$F=\left(\begin{array}{ccc}\begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}\frac{{\beta }_{dP}{A}_d\left({{\mu }_d+\alpha }_d\right)}{{\mu }_d\left({{\mu }_d+\alpha }_d+v_d\right)}\\ 0\\ \begin{array}{c}0\\ \frac{{\beta }_{hP}{A}_h\left({{\mu }_h+\alpha }_h\right)}{{\mu }_h\left({{\mu }_h+\alpha }_h+v_h\right)}\\ 0\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\frac{{\beta }_{dF}{A}_d\left({{\mu }_d+\alpha }_d\right)}{{\mu }_d\left({{\mu }_d+\alpha }_d+v_d\right)}\\ 0\\ \begin{array}{c}0\\ \frac{{\beta }_{hF}{A}_h\left({{\mu }_h+\alpha }_h\right)}{{\mu }_h\left({{\mu }_h+\alpha }_h+v_h\right)}\\ 0\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}\end{array}\end{array}\right)$$

(3.18)

$$V=\left(\begin{array}{ccc}\begin{array}{c}({\mu }_d+{\delta }_d+{\theta }_d)\\ -{\delta }_d\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ \left({\mu }_d+{\rho }_d+{\mathit{c}}_{\mathit{d}}\right)\\ \begin{array}{c}-{\rho }_d\\ 0\\ 0\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}0\\ 0\\ \begin{array}{c}\left({\mu }_d+{\epsilon }_d+{\mathit{c}}_{\mathit{d}}\right)\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ \left({{\mu }_h+\delta }_h+{\theta }_h\right)\\ -{\delta }_h\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ \left({\mu }_h+{\epsilon }_h\right)\end{array}\end{array}\end{array}\end{array}\right)$$

The inverse of matrix V is given by

$$V^{-1}=\left(\begin{array}{ccc}\begin{array}{c}\frac{1}{({\mu }_d+{\delta }_d+{\theta }_d)}\\ \frac{{\delta }_d}{({\mu }_d+{\delta }_d+{\theta }_d)\left({\mu }_d+{\rho }_d+c_d\right)}\\ \begin{array}{c}\frac{{\rho }_d{\delta }_d}{\left({\mu }_d+{\rho }_d+c_d\right)({\mu }_d+{\delta }_d+{\theta }_d)\left({\mu }_d+{\epsilon }_d+c_d\right)}\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ \frac{1}{\left({\mu }_d+{\rho }_d+c_d\right)}\\ \begin{array}{c}\frac{{\rho }_d}{\left({\mu }_d+{\rho }_d+c_d\right)\left({\mu }_d+{\epsilon }_d+c_d\right)} \\ 0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}\frac{1}{\left({\mu }_d+{\epsilon }_d+c_d\right)}\\ 0\\ 0\end{array}\end{array}\end{array}\right.$$

$$\left.\begin{array}{cc}\begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ \frac{1}{\left({{\mu }_h+\delta }_h+{\theta }_h\right)}\\ \frac{{\delta }_h}{\left({{\mu }_h+\delta }_h+{\theta }_h\right)\left({\mu }_d+{\epsilon }_d\right)}\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ \frac{1}{\left({\mu }_d+{\epsilon }_d\right)}\end{array}\end{array}\end{array}\right)$$

(3.19)

The symbol $\mathsf{\rho }\left(\mathrm{F}{\mathrm{V}}^{-1}\right)$ represents the dominant eigenvalue in magnitude, or spectral radius, of a matrix $\mathrm{F}{\mathrm{V}}^{-1}$, and take$\mathrm{F}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{V}}^{-1}$ from equation (3.18) and (3.19), thus effective reproductive number defend by $R_e=\rho \left(F{V}^{-1}\right)$Therefore, the effective reproduction number is as follows:

$${\mathbf{R}}_{\mathbf{e}}=\frac{{\mathit{A}}_{\mathit{d}}({\mathsf{\beta }}_{\mathbf{d}\mathbf{P}}\left({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\epsilon }}_{\mathbf{d}}+{\mathbf{c}}_{\mathbf{d}}\right)+{\mathsf{\rho }}_{\mathbf{d}}{\mathsf{\beta }}_{\mathbf{d}\mathbf{F}}){\mathsf{\delta }}_{\mathbf{d}}\left({{\mathsf{\mu }}_{\mathbf{d}}+\mathsf{\alpha }}_{\mathbf{d}}\right)}{{\mathsf{\mu }}_{\mathbf{d}}\left({{\mathsf{\mu }}_{\mathbf{d}}+{\mathbf{v}}_{\mathbf{d}}+\mathsf{\alpha }}_{\mathbf{d}}\right)\left({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\rho }}_{\mathbf{d}}+{\mathbf{c}}_{\mathbf{d}}\right)\left({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\epsilon }}_{\mathbf{d}}+{\mathbf{c}}_{\mathbf{d}}\right)({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\delta }}_{\mathbf{d}}+{\mathit{\theta }}_{\mathit{d}})}$$

(3.20)

When we now replace the parameter values from Table 2 with the expression given in Equations (3.17) and (3.20) for the basic and effective reproductive numbers, we obtain the following results.

$${\mathbf{R}}_0=\frac{{\mathit{A}}_{\mathit{d}}{\mathsf{\delta }}_{\mathbf{d}}({\mathsf{\beta }}_{\mathbf{d}\mathbf{P}}\left({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\epsilon }}_{\mathbf{d}}\right)+{\mathsf{\rho }}_{\mathbf{d}}{\mathsf{\beta }}_{\mathbf{d}\mathbf{F}})}{{\mathsf{\mu }}_{\mathbf{d}}\left({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\rho }}_{\mathbf{d}}\right)\left({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\epsilon }}_{\mathbf{d}}\right)({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\delta }}_{\mathbf{d}})}=2.74$$

The fact that the basic reproductive number ${\mathbf{R}}_0$ is larger than one in the absence of any control measures indicates that the disease will keep spreading in the population.

$${\mathbf{R}}_{\mathbf{e}}=\frac{{\mathit{A}}_{\mathit{d}}({\mathsf{\beta }}_{\mathbf{d}\mathbf{P}}\left({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\epsilon }}_{\mathbf{d}}+{\mathbf{c}}_{\mathbf{d}}\right)+{\mathsf{\rho }}_{\mathbf{d}}{\mathsf{\beta }}_{\mathbf{d}\mathbf{F}}){\mathsf{\delta }}_{\mathbf{d}}\left({{\mathsf{\mu }}_{\mathbf{d}}+\mathsf{\alpha }}_{\mathbf{d}}\right)}{{\mathsf{\mu }}_{\mathbf{d}}\left({{\mathsf{\mu }}_{\mathbf{d}}+{\mathbf{v}}_{\mathbf{d}}+\mathsf{\alpha }}_{\mathbf{d}}\right)\left({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\rho }}_{\mathbf{d}}+{\mathbf{c}}_{\mathbf{d}}\right)\left({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\epsilon }}_{\mathbf{d}}+{\mathbf{c}}_{\mathbf{d}}\right)({\mathsf{\mu }}_{\mathbf{d}}+{\mathsf{\delta }}_{\mathbf{d}}+{\mathit{\theta }}_{\mathit{d}})}=1.78$$

The effective reproductive number ${\mathbf{R}}_{\mathbf{e}}$ is larger than one given the current vaccination coverage, indicating that the disease is still present. This suggests that greater measures should be taken to prevent the spread of rabies.

3.6. Global Stability of Disease-Free Equilibrium Points

We used the method suggested by [7] to investigate the global stability of the disease-free equilibrium point of the system (1).

The structure of our model, as shown in system (1), is as follows.

$$\left.\begin{array}{c}\frac{\mathit{d}\mathit{X}}{\mathit{d}\mathit{t}}={\mathit{B}}_0\left(\mathit{X}-{\mathit{X}}_0\right)+{\mathit{B}}_{1}Y\\ \frac{\mathit{d}\mathit{Y}}{\mathit{d}\mathit{t}}={\mathit{B}}_{2}Y \end{array}\right\}$$

(3.21)

The classifications of susceptible and vaccinated individuals are represented by $X\in {\mathbb{R}}_{+}^m$. Classes of exposed and infectious individuals are represented by $Y\in {\mathbb{R}}_{+}^n$ A vector at DFE point ${\epsilon }_0$ of the vector length as $X$ is represented by Using system (1) as an example, we define:

$$\mathrm{X}=\left[\begin{array}{c}{\mathrm{S}}_{\mathrm{d}}\\ {\mathrm{R}}_{\mathrm{d}}\\ {\mathrm{S}}_{\mathrm{h}}\\ {\mathrm{R}}_{\mathrm{h}}\end{array}\right], \mathrm{Y}=\left[\begin{array}{c}{\mathrm{E}}_{\mathrm{d}}\\ {\mathrm{I}}_{\mathrm{d}\mathrm{P}}\\ {\mathrm{I}}_{\mathrm{d}\mathrm{F}}\\ E_d\\ I_h\end{array}\right] \mathrm{and} X_{{\epsilon }_0}=\left[\begin{array}{c}\frac{A_d\left({{\mu }_d+\alpha }_d\right)}{{\mu }_d\left({{\mu }_d+\alpha }_d+v_d\right)}\\ \frac{A_d{v}_d}{{\mu }_d\left({{\mu }_d+\alpha }_d+v_d\right)}\\ \frac{A_h\left({{\mu }_h+\alpha }_h\right)}{{\mu }_h\left({{\mu }_h+\alpha }_h+v_h\right)}\\ \frac{A_h{v}_h}{{\mu }_h\left({{\mu }_h+\alpha }_h+v_h\right)}\end{array}\right], \mathrm{X}-{\mathrm{X}}_0=\left[\begin{array}{c}{S}_d-\frac{A_d\left({{\mu }_d+\alpha }_d\right)}{{\mu }_d\left({{\mu }_d+\alpha }_d+v_d\right)}\\ R_d-\frac{A_d{v}_d}{{\mu }_d\left({{\mu }_d+\alpha }_d+v_d\right)}\\ S_h-\frac{A_h\left({{\mu }_h+\alpha }_h\right)}{{\mu }_h\left({{\mu }_h+\alpha }_h+v_h\right)}\\ R_d-\frac{A_h{v}_h}{{\mu }_h\left({{\mu }_h+\alpha }_h+v_h\right)}\end{array}\right]$$

To confirm whether the disease-free equilibrium is globally stable, we have to show that:

i.

B₀ should be a matrix whose eigenvalues are real and negative; and

ii.

B_z should be a Metzler matrix.

Using system (1) and the representation in 3.21 the first equation can be rewritten as shown below:

$$\left[\begin{array}{c}{\mathrm{A}}_d{+\alpha }_d{R}_d-{\left({{\beta }_{dP}I}_{dP}+{{\beta }_{dF}I}_{dF}+{\mu }_d+v_d\right)S}_d\\ v_d{S}_d+{\theta }_d{E}_d-\left({{\mu }_d+\alpha }_d\right)R_d\\ {\mathrm{A}}_h{+\alpha }_h{R}_h-{({{\beta }_{hP}I}_{dP}+{{\beta }_{hF}I}_{dF}+{\mu }_h+v_h)S}_h\\ v_h{S}_h+{\theta }_h{E}_h-\left({{\mu }_h+\alpha }_h\right)R_h\end{array}\right]={\mathit{B}}_0\left[\begin{array}{c}{S}_d-\frac{A_d\left({{\mu }_d+\alpha }_d\right)}{{\mu }_d\left({{\mu }_d+\alpha }_d+v_d\right)}\\ R_d-\frac{A_d{v}_d}{{\mu }_d\left({{\mu }_d+\alpha }_d+v_d\right)}\\ S_h-\frac{A_h\left({{\mu }_h+\alpha }_h\right)}{{\mu }_h\left({{\mu }_h+\alpha }_h+v_h\right)}\\ R_d-\frac{A_h{v}_h}{{\mu }_h\left({{\mu }_h+\alpha }_h+v_h\right)}\end{array}\right]+B_1\left[\begin{array}{c}{\mathrm{E}}_{\mathrm{d}}\\ {\mathrm{I}}_{\mathrm{d}\mathrm{P}}\\ {\mathrm{I}}_{\mathrm{d}\mathrm{F}}\\ E_d\\ I_h\end{array}\right]$$

(3.22)

The equation (3.22)’s left side is rewritten as

$$\left[\begin{array}{c}{\mathrm{A}}_d{+\alpha }_d{R}_d-{\left(+{\mu }_d+v_d\right)S}_d\\ v_d{S}_d-\left({{\mu }_d+\alpha }_d\right)R_d\\ {\mathrm{A}}_h{+\alpha }_h{R}_h-{(+{\mu }_h+v_h)S}_h\\ v_h{S}_h-\left({{\mu }_h+\alpha }_h\right)R_h\end{array}\right]+\left[\begin{array}{c}-{\left({{\beta }_{dP}I}_{dP}+{{\beta }_{dF}I}_{dF}\right)S}_d\\ {\theta }_d{E}_d\\ -{({{\beta }_{hP}I}_{dP}+{{\beta }_{hF}I}_{dF})S}_h\\ +{\theta }_h{E}_h\end{array}\right]={\mathit{B}}_0\left[\begin{array}{c}{S}_d-S_d^0\\ R_d-R_d^0\\ S_h-S_h^0\\ R_d-R_h^0\end{array}\right]+B_1\left[\begin{array}{c}{\mathrm{E}}_{\mathrm{d}}\\ {\mathrm{I}}_{\mathrm{d}\mathrm{P}}\\ {\mathrm{I}}_{\mathrm{d}\mathrm{F}}\\ E_d\\ I_h\end{array}\right]$$

(3.23)

Using the state variables $S_d, R_d, S_h, R_h$, and ${\mathrm{E}}_{\mathrm{d}}, {\mathrm{I}}_{\mathrm{d}\mathrm{P}}, {\mathrm{I}}_{\mathrm{d}\mathrm{F}}, and E_h$ of the Jacobian matrix, we obtain the first vector and the second vector on the left side of the equation (3.23).

$${\mathit{B}}_0=\left[\begin{array}{cccc}-\left(+{\mu }_d+v_d\right)& { \alpha }_h& 0& 0\\ v_d& -\left({{\mu }_d+\alpha }_d\right)& 0& 0\\ 0& 0& -({\mu }_h+v_h)& { \alpha }_h\\ 0& 0& v_h& -\left({{\mu }_h+\alpha }_h\right)\end{array}\right],{\mathit{B}}_1=\left[\begin{array}{cccc}0& -{\beta }_{dP}{S}_d& -{\beta }_{dF}{S}_d& 0\\ {\theta }_d& 0& 0& 0\\ 0& -{\beta }_{hP}{S}_h& -{\beta }_{hF}{S}_h& 0\\ 0& 0& 0& {\theta }_h\end{array}\right]$$

(3.23)

Next, show that a matrix ${\mathit{B}}_0$ has real, negative eigenvalues.

$\left|{\mathit{B}}_0-\mathit{I}\mathit{\lambda }\right|=0$ is the characteristic equation of matrix $B_0$.

$$\left|\begin{array}{cccc}-\left({\mu }_d+v_d+\mathit{\lambda }\right)& { \alpha }_h& 0& 0\\ v_d& -\left({{\mu }_d+\alpha }_d+\mathit{\lambda }\right)& 0& 0\\ 0& 0& -({\mu }_h+v_h+\mathit{\lambda })& { \alpha }_h\\ 0& 0& v_h& -\left({{\mu }_h+\alpha }_h+\mathit{\lambda }\right)\end{array}\right|=0$$

$$-\left({\mu }_d+v_d+\mathit{\lambda }\right)\left|\begin{array}{ccc}-\left({{\mu }_d+\alpha }_d+\mathit{\lambda }\right)& 0& 0\\ 0& -\left({\mu }_h+v_h+\mathit{\lambda }\right)& { \alpha }_h\\ 0& v_h& -\left({{\mu }_h+\alpha }_h+\mathit{\lambda }\right)\end{array}\right|=0$$

$$\left({\mu }_d+v_d+\mathit{\lambda }\right)\left({{\mu }_d+\alpha }_d+\mathit{\lambda }\right)\left|\begin{array}{cc}-\left({\mu }_h+v_h+\mathit{\lambda }\right)& { \alpha }_h\\ v_h& -\left({{\mu }_h+\alpha }_h+\mathit{\lambda }\right)\end{array}\right|=0$$

$$\left({\mu }_d+v_d+\mathit{\lambda }\right)\left({{\mu }_d+\alpha }_d+\mathit{\lambda }\right)\left(\left({\mu }_h+v_h+\mathit{\lambda }\right)\left({{\mu }_h+\alpha }_h+\mathit{\lambda }\right)-v_h{ \alpha }_h\right)=0$$

$$\left({\mu }_d+v_d+\mathit{\lambda }\right)\left({{\mu }_d+\alpha }_d+\mathit{\lambda }\right)\left({\mu }_h+\mathit{\lambda }\right)\left(({\mu }_h+v_h+{ \alpha }_h)+\mathit{\lambda }\right)=0$$

As a result, the following Eigen values exist.

$$\left\{\begin{array}{c}{\lambda }_1=-\left({\mu }_d+v_d\right) \\ {\lambda }_2=-\left({{\mu }_d+\alpha }_d\right) \\ {\lambda }_3={-\mu }_h \\ { \lambda }_4=-\left({\mu }_h+v_h+{ \alpha }_h\right) \end{array}\right.$$

(3.24)

The second equation can be remake as follows using system (1) and the representation in 3.21:

$$\left[\begin{array}{c}{({\beta }_{dP}{I}_{dp}+{{\beta }_{dF}I}_{dF})S}_d-\left({\mu }_d+{\delta }_d\right)E_d\\ {\delta }_d{E}_d-\left({\mu }_d+{\rho }_d\right)I_{dP}\\ {\rho }_d{I}_{dP}-\left({\mu }_d+{\epsilon }_d\right)I_{dF}\\ {\left({{\beta }_{hP}I}_{dP}+{{\beta }_{hF}I}_{dF}\right)S}_h-\left({{\mu }_h+\delta }_h\right)E_h\\ {\delta }_h{E}_h-\left({\mu }_h+{\epsilon }_h\right)I_h\end{array}\right]=B_2\left[\begin{array}{c}{\mathrm{E}}_{\mathrm{d}}\\ {\mathrm{I}}_{\mathrm{d}\mathrm{P}}\\ {\mathrm{I}}_{\mathrm{d}\mathrm{F}}\\ E_d\\ I_h\end{array}\right]$$

(3.25)

Using the Jacobian matrix’s elements ${\mathrm{E}}_{\mathrm{d}}, {\mathrm{I}}_{\mathrm{d}\mathrm{P}}, {\mathrm{I}}_{\mathrm{d}\mathrm{F}}, E_{h, }$ $I_h$ the left side of equation (3.25) becomes:

$$B_2=\left[\begin{array}{ccccc}-\left({\mu }_d+{\delta }_d\right)& {\beta }_{dP}{S}_d& {\beta }_{dF}{S}_d& 0& 0\\ {\delta }_d& -\left({\mu }_d+{\rho }_d\right)& 0& 0& 0\\ 0& {\rho }_d& -\left({\mu }_d+{\epsilon }_d\right)& 0& 0\\ 0& {\beta }_{hP}{S}_h& {\beta }_{hF}{S}_h& -\left({{\mu }_h+\delta }_h\right)& 0\\ 0& 0& 0& {\delta }_h& -\left({\mu }_h+{\epsilon }_h\right)\end{array}\right]$$

(3.26)

Thus, $B_2$ is a Metzler matrix

Therefore, from equation (3.24) all Eigenvalues of matrix $B_0$ are $-\left({\mu }_d+v_d\right),$ $-\left({{\mu }_d+\alpha }_d\right),$ ${-\mu }_h$ and $-\left({\mu }_h+v_h+{ \alpha }_h\right)$ which are real and negative. Additionally, matrix $B_2$ off-diagonal elements are non-negative according to equation (3.26) since all of the parameters are positive, showing that the matrix is a Metzler matrix. This further proves that the region Ω is globally asymptotically stable for the disease-free equilibrium points of the system (1). This leads us to the critical theorem that follows.

Theorem 4:The disease-free equilibrium point is globally asymptotically stable in the region Ω.if${\mathrm{R}}_{\mathrm{e}}<1$and unstable in the region Ω. if${\mathrm{R}}_{\mathrm{e}}>1$

4. Sensitivity Analysis

We can determine how much a state variable changes in response to a changing parameter by using sensitivity indices. Sensitivity analysis is commonly used to figure out how sensitive model predictions are to parameter values because errors in data collection and expected parameter values occur frequently. Therefore, we use it to identify the characteristics that have a major impact on $R_e$ and should be the focus of attention for intervention strategies [8, 9]. The relationship between the parameters and $R_e$ is inversely proportional if the result is negative. To decrease the size of the effect of modifying that parameter, we can in this case take the modulus of the sensitivity index. A positive sensitivity index, on the other hand, indicates that $R_e$ varies in direct proportion to the parameter [9, 10]. The equation (2.20) provides the explicit expression of ${\mathrm{R}}_{\mathrm{e}}=\frac{{\mathrm{A}}_{\mathrm{d}}({\mathsf{\beta }}_{\mathrm{d}\mathrm{P}}\left({\mathsf{\mu }}_{\mathrm{d}}+{\mathsf{\epsilon }}_{\mathrm{d}}+{\mathrm{c}}_{\mathrm{d}}\right)+{\mathsf{\rho }}_{\mathrm{d}}{\mathsf{\beta }}_{\mathrm{d}\mathrm{F}}){\mathsf{\delta }}_{\mathrm{d}}\left({{\mathsf{\mu }}_{\mathrm{d}}+\mathsf{\alpha }}_{\mathrm{d}}\right)}{{\mathsf{\mu }}_{\mathrm{d}}\left({{\mathsf{\mu }}_{\mathrm{d}}+{\mathrm{v}}_{\mathrm{d}}+\mathsf{\alpha }}_{\mathrm{d}}\right)\left({\mathsf{\mu }}_{\mathrm{d}}+{\mathsf{\rho }}_{\mathrm{d}}+{\mathrm{c}}_{\mathrm{d}}\right)\left({\mathsf{\mu }}_{\mathrm{d}}+{\mathsf{\epsilon }}_{\mathrm{d}}+{\mathrm{c}}_{\mathrm{d}}\right)({\mathsf{\mu }}_{\mathrm{d}}+{\mathsf{\delta }}_{\mathrm{d}}+{\mathsf{\theta }}_{\mathrm{d}})}$. Since $R_e$only depends on eleven parameters $A_d,{\beta }_{dP},{\beta }_{dF}$, ${\delta }_d, {\rho }_d, v_d, {\theta }_d, cd {\epsilon }_d, { \mu }_d,$ and ${\alpha }_d$. we can use the normalized forward sensitivity index to calculate an analytical equation for$R_e$ sensitivity to each of these parameters [9, 13].

The normalized sensitivity index $S_{\omega }^{R_e}$ is defined as

$$S_{\omega }^{R_e}=\frac{\partial R_e}{\partial \omega }\times \frac{\omega }{R_e}$$

Thus, using the values in Table 2, the normalized sensitivity indices for the eleven parameters are derived. We use literature data and certain assumptions for the parameter values.

Based on our model, the $R_e$ sensitivity index concerning $A_d$ is as follows:

$${\mathrm{S}}_{{\mathrm{A}}_{\mathrm{d}}}^{{\mathrm{R}}_{\mathrm{e}}}=\frac{\partial {\mathrm{R}}_{\mathrm{e}}}{\partial {\mathrm{A}}_{\mathrm{d}}}\times \frac{{\mathrm{A}}_{\mathrm{d}}}{{\mathrm{R}}_{\mathrm{e}}}=+1$$

Similarly, the sensitivity index concerning${\mathsf{\beta }}_{\mathrm{d}\mathrm{P}},{{\mathsf{\beta }}_{\mathrm{d}\mathrm{F}},\delta }_d,{\rho }_d,v_d,{\theta }_d,c_d{\epsilon }_d,{\mu }_d,$ and ${\alpha }_d$ is given by: In the same way, concerning ${\mathsf{\beta }}_{\mathrm{d}\mathrm{P}}, {{\mathsf{\beta }}_{\mathrm{d}\mathrm{F}}, \delta }_d, {\rho }_d, v_d, {\theta }_d,{ c}_d {\epsilon }_d, { \mu }_d,$ and ${\alpha }_d$, the sensitivity index is provided by:

$${\mathrm{S}}_{{\mathsf{\beta }}_{\mathrm{d}\mathrm{P}}}^{{\mathrm{R}}_{\mathrm{e}}}=\frac{\partial {\mathrm{R}}_{\mathrm{e}}}{\partial {\mathsf{\beta }}_{\mathrm{d}\mathrm{P}}}\times \frac{{\mathsf{\beta }}_{\mathrm{d}\mathrm{P}}}{{\mathrm{R}}_{\mathrm{e}}}=0.001747639115$$

$$S_{{\mathsf{\beta }}_{\mathrm{d}\mathrm{F}}}^{R_e}=\frac{\partial R_e}{\partial {\mathsf{\beta }}_{\mathrm{d}\mathrm{F}} }\times \frac{{\mathsf{\beta }}_{\mathrm{d}\mathrm{F}} }{R_e}=0.9982523609$$

$$S_{{\mathsf{\delta }}_{\mathrm{d}}}^{R_e}=\frac{\partial R_e}{\partial {\mathsf{\delta }}_{\mathrm{d}}}\times \frac{{\mathsf{\delta }}_{\mathrm{d}}}{R_e}=0.6056304706$$

$$S_{{\alpha }_d}^{R_e}=\frac{\partial R_e}{\partial {\alpha }_d}\times \frac{{\alpha }_d}{R_e}=0.1812729130$$

$${ S}_{{ v}_d}^{R_e}=\frac{\partial R_e}{\partial v_d}\times \frac{{ v}_d}{R_e}=-0.1914241960$$

$$S_{c_d}^{R_e}=\frac{\partial R_e}{\partial c_d }\times \frac{c_d }{R_e}=-0.3144890128$$

$${ S}_{{\mu }_d}^{R_e}=\frac{\partial R_e}{\partial {\mu }_d}\times \frac{{\mu }_d}{R_e}=-1.220608197$$

$$S_{{\theta }_d}^{R_e}=\frac{\partial R_e}{\partial {\theta }_d }\times \frac{{\theta }_d }{R_e}=-0.4731488052$$

$${ S}_{{ \epsilon }_d}^{R_e}=\frac{\partial R_e}{\partial {\epsilon }_d}\times \frac{{ \epsilon }_d}{R_e}=-0.8081706290$$

$$S_{{\mathsf{\rho }}_{\mathrm{d}}}^{R_e}=\frac{\partial R_e}{\partial {\mathsf{\rho }}_{\mathrm{d}} }\times \frac{{\mathsf{\rho }}_{\mathrm{d}} }{R_e}=0.2209374584$$

4.1. Interpretation of Sensitivity Indices

We observe the sensitivity index values’ sign. The parameters ($A_d,{\beta }_{dP},{\beta }_{dF},{ \delta }_d,$and${ \alpha }_d$) which have positive indices, indicate that they have a significant effect on the spread of the disease within the community. As their values increase, the community’s average number of secondary case infections increases due to an increase in the effective reproduction number ($R_e$) i.e., the parameters ($A_d,{\beta }_{dP},{\beta }_{dF},{ \delta }_d,$and${ \alpha }_d$) directly proportional to $R_e$. In addition, the parameters (${ v}_d,{\theta }_d,c_d,{\rho }_d$, and ${\epsilon }_d$) with negative sensitivity indices have an impact on decreasing the disease burden in the community when their values increase while the others remain unchanged. Also, as their values increase, the effective reproduction number ($R_e$) decreases, minimizing the disease’s widespread within the community i.e., inversely proportional to $R_e$ [13]. The relationship between the parameters ($A_d,{\beta }_{dP},{\beta }_{dF},{ \delta }_d, {\theta }_d,{ \delta }_d,{\epsilon }_d$) and the effective reproductive number is numerically simulated and is displayed in the figures below.

Figure 2. Shows a direct relationship between the effective reproduction number (R_e) and Furious phase dog-to-dog transmission rate ${\beta }_{dF}$. Also, Figure 3 shows a direct relationship between the effective reproduction number $(R_e$) and the Prodromal phase dog-to-dog transmission rate. Both the infection rates ${\beta }_{dF},$ and ${\beta }_{dP}$ values varied as shown in Figure 2 and Figure 3. The result agrees with reality in the sense that as the rate of infection rabies increases, the number of individuals that will be infected in the population also increases thereby increasing the effective reproduction number.

Figure 4. Shows that the effective reproduction numbe (R_e) is an increasing function of of the recruitment rate $A_d$. Here the values of the recruitment rate $A_d$were varied as shown in the figure. This shows that direct relationship between the annual birth rate of the dog population and the effective reproductive number. In Figure 6, the effective reproduction number ${ R}_e$ varies directly with the rate of transfer from Exposed to Prodromal stage${\delta }_d$. Values of ${\delta }_d$were varied as shown in the figure.

Figure 4. Effect of parameter $A_d$ on $R_e$

Figure 4. Effect of parameter $A_d$ on $R_e$

Figure 5. Effect of parameter ${\delta }_d$ on $R_e$

Figure 5. Effect of parameter ${\delta }_d$ on $R_e$

Figure 6. Effect of parameter ${\theta }_d$ on $R_e$

Figure 6. Effect of parameter ${\theta }_d$ on $R_e$

Figure 6 shows an inverse relationship between the effective reproduction number ${ R}_e$ and the recovery rate ${\theta }_d$due to post-exposed treatment. This is true because if proper treatment is given to Exposed individual dogs, fewer people will be infected thereby reducing the effective reproduction number${ R}_e$.

Figure 7. Effect of parameter ${\epsilon }_d$ on $R_e$

Figure 7. Effect of parameter ${\epsilon }_d$ on $R_e$

In Figure 7. The effective reproduction number ${ R}_e$ is a decreasing function of rabies-induced death ${\epsilon }_d$. Values of ${\epsilon }_d$were varied as shown in the figure. This also agrees with reality because as more infectious individuals die as a result of rabies infection, the infection rate will also be reduced.

5. Analysis of the Optimal Control Model

This section does a detailed study of the time-dependent rabies model, given by the system (1). Pontryagin’s Maximum Principle [14], which has been widely used in mathematical models of biological processes including optimal forms the basis of the analysis [15, 16]. The following goal or cost-functional is used to decrease the populations of infectious dogs in the prodromal phase,${ I}_{dP}$, and infectious dogs in the furious phase ${ I}_{dF}$, as well as exposed humans $E_h$and infectious humans${ I}_h$. Our goal is to minimize the objective function by identifying the most efficient controls, such as $u_i\left(t\right),(i=1, 2, 3, 4):$.

Our objective function becomes

$$J=\underset{u\in U}{\min }{\int }_0^T\left({ A_{1}E}_d+{A_{2}E}_h+A_3{I}_{dP}+{A_{4}I}_{dF}+{A_{5}I}_h+\frac{1}{2} \sum _{i=1}^4{B}_{i }{u}_i^2 \left(t\right)\right)dt$$

(4.1)

subject to the model (1),

$$\left.\begin{array}{c} \frac{d{S}_d}{dt}={\mathrm{A}}_d-{\left({{\beta }_{dP}I}_{dP}+{{\beta }_{dF}I}_{dF}\right)S}_d-\left({\mu }_d+u_1\right)S_d {+\alpha }_d{R}_d\\ \frac{d{E}_d}{dt}={({\beta }_{dP}{I}_{dp}+{{\beta }_{dF}I}_{dF})S}_d-\left({\mu }_d+{\delta }_d+u_2\right)E_d\\ \frac{d{I}_{dp}}{dt}={\delta }_d{E}_d-\left({\mu }_d+{\rho }_d+{\mathit{c}}_{\mathit{d}}\right)I_{dP}\\ \frac{d{I}_{dF}}{dt}={\rho }_d{I}_{dP}-\left({\mu }_d+{\epsilon }_d+{\mathit{c}}_{\mathit{d}}\right)I_{dF}\\ \frac{d{R}_d}{dt}=u_1{S}_d+u_2{E}_d-\left({{\mu }_d+\alpha }_d\right)R_d \\ \frac{d{S}_h}{dt}={\mathrm{A}}_h-{({{\beta }_{hP}I}_{dP}+{{\beta }_{hF}I}_{dF})S}_h-\left({\mu }_h+u_3\right)S_h{+\alpha }_h{R}_h\\ \frac{d{E}_h}{dt}={({{\beta }_{hP}I}_{dP}+{{\beta }_{hF}I}_{dF})S}_h-\left({{\mu }_h+\delta }_h+u_4\right)E_h\\ \frac{d{I}_h}{dt}={\delta }_h{E}_h-\left({\mu }_h+{\epsilon }_h\right)I_h \\ \frac{d{R}_h}{dt}=u_3{S}_h+u_4{E}_h-\left({{\mu }_h+\alpha }_h\right)R_h \end{array}\right\}$$

where $A_1$, and A₂, are the weight constants of the exposed dog and human class, $A_3,$and $A_4$ are the weight constants of prodromal and Furious phase infectious dog classes, and $A_5$ are the weight constants of the infectious human class. Similarly, $B_1{u}_{ 1}^2, {B_{2}u}_{ 2}^2, {B_{3}u}_{ 3}^2,$and ${B_{4}u}_{ 4}^2$ describe the cost associated with rabies vaccination and treatment. In this work, as in other studies [17, 18], the cost control functions take a quadratic form. We want to find the optimal control $u_i^{*}=u_1^{*},u_2^{*}, u_3^{*}, u_4^{*}$such as

$$J\left(u_1^{*},u_2^{*}, u_3^{*}, u_4^{*} \right)=\min J\{\left(u_1, u_2, u_3, u_4\right) | u_i\in U \mathrm{for} i=1, 2, 3, 4\}$$

(4.2)

subject to the control set provided by the dynamical system in (1)

$U=\{u_1, u_2, u_3, u_4\}$ such that, $0\le u_i\le 1$ for $i=1, 2, 3, 4\};\mathrm{L}\mathrm{e}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{g}\mathrm{u}\mathrm{e}\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e} {\forall }_t\in \left[0, T\right],$ the non-empty control set.

5.1. Characterization of the Optimal Control

The Pontryagin’s Maximum Principle [14] is used to try and obtain the necessary conditions for the optimal control of rabies governed by the non-autonomous system. This principle becomes the state system (1), with the objective functional (4.1) and (4.2), into a pointwise minimizing problem concerning the controls $u_1, u_2, u_3 and u_4$.

The Hamiltonian equation with state variables $S_d=S_d^{\mathrm{*}},E_d=E_d^{\mathrm{*}},$ $I_{dP}=I_{dP}^{\mathrm{*}}$, $I_{dF}=I_{dF}^{\mathrm{*}}$, $R_d=R_d^{\mathrm{*}}$, $S_h=S_h^{\mathrm{*}},E_h=E_h^{\mathrm{*}},$ $I_h=I_h^{\mathrm{*}},$ and $R_h=R_h^{\mathrm{*}}$ is formed as

$H=Iintegrand+\left(adjoins\right)\times \left(RHS\right)$ of the

$$\begin{array}{c} H=A_1{E}_d^{*}+A_2{E}_h^{*}+A_3{I}_{dP}^{*}+A_4{I}_{dF}^{*}+A_5{I}_h^{*} \\ +\frac{1}{2}\left[ B_1{u}_{ 1}^2+{B_{2}u}_{ 2}^2+{B_{3}u}_{ 3}^2+{B_{4}u}_{ 4}^2\right] \\ +{\lambda }_1\left[A_d-\left({\beta }_{dP}{I}_{dP}^{*}+{\beta }_{dF}{I}_{dF}^{*}\right)S_d^{*}-\left({\mu }_d+u_1\right)S_d^{*} {+\alpha }_d{R}_d^{*}\right] 4.3 \\ +{\lambda }_2\left[\left({\beta }_{dP}{I}_{dP}^{*}+{\beta }_{dF}{I}_{dF}^{*}\right)S_d^{*}-\left({\mu }_d+{\delta }_d+u_2\right)E_d^{*}\right] \\ +{\lambda }_3\left[{\delta }_d{E}_d^{*}-\left({\mu }_d+{\rho }_d+c_d\right)I_{dP}^{*} \right] \\ +{\lambda }_4\left[{\rho }_d{I}_{dP}^{*}-\left({\mu }_d+{\epsilon }_d+c_d\right)I_{dF}^{*} \right] \\ +{\lambda }_5\left[u_1{S}_d^{*}+u_2{E}_d^{*}-\left({{\mu }_d+\alpha }_d\right)R_d^{*} \right] \\ +{\lambda }_6[A_h-\left({{\beta }_{hP}I}_{dP}^{*}+{\beta }_{hF}{I}_{dF}^{*}\right)S_h^{*}-{(\mu }_h+u_3)S_h^{*}{+\alpha }_h{R}_h^{*} ] \\ +{\lambda }_7\left[\left({{\beta }_{hP}I}_{dP}^{*}+{\beta }_{hF}{I}_{dF}^{*}\right)S_h^{*}-\left({{\mu }_h+\delta }_h+u_4\right)E_h^{*} \right] \\ +{\lambda }_8\left[ {\delta }_h{E}_h^{*}-\left({\mu }_h+{\epsilon }_h\right)I_h^{*}\right] \\ +{\lambda }_9\left[u_3{S}_h^{*}+u_4{E}_h^{*}-\left({{\mu }_h+\alpha }_h\right)R_h^{*} \right] \end{array}$$

where the adjoint (co-state) variables are ${\mathsf{\lambda }}_{\mathrm{i}},$i = 1, 2, 3, 5, 6, 7, 8, 9. The necessary conditions for the optimal control by the following result.

Theorem4.2:Given an optimal control quadruple${\mathrm{u}}_1^{\mathrm{*}},{\mathrm{u}}_2^{\mathrm{*}},{\mathrm{u}}_3^{\mathrm{*}},{\mathrm{u}}_4^{\mathrm{*}}$that minimizes objective functional (4.1) over the control set U subject to the state system (1), then there exist adjoint variables${\mathsf{\lambda }}_1,{\mathsf{\lambda }}_{2,}{\mathsf{\lambda }}_{3,}{\mathsf{\lambda }}_4,{\mathsf{\lambda }}_5,{\mathsf{\lambda }}_6,{\mathsf{\lambda }}_7,{\mathsf{\lambda }}_8,{\mathsf{\lambda }}_9$satisfying

$$\frac{d{\lambda }_1}{dt}={\lambda }_1\left[\left({\beta }_{dP}{I}_{dP}^{*}+{\beta }_{dF}{I}_{dF}^{*}\right)+\left({\mu }_d+u_1\right)\right]-{\lambda }_2\left({\beta }_{dP}{I}_{dP}^{*}+{\beta }_{dF}{I}_{dF}^{*}\right)-{\lambda }_5{u}_1$$

$$\frac{d{\lambda }_2}{dt}={\lambda }_2\left({\mu }_d+{\delta }_d+u_2\right)-{\lambda }_3{\delta }_d-{\lambda }_5{u}_2-A_1$$

$$\frac{d{\lambda }_3}{dt}={\lambda }_1{\beta }_{dP}{S}_d^{*}+{\lambda }_3\left({\mu }_d+{\rho }_d+c_d\right)+{\lambda }_6{\beta }_{hP}{S}_h^{*}-{\lambda }_2{\beta }_{dP}{S}_d^{*}$$

$$-{\lambda }_4{\rho }_d-{\lambda }_7{\beta }_{hP}{S}_h^{*}-A_3 4.4$$

$$\frac{d{\lambda }_4}{dt}={\lambda }_1{\beta }_{dF}{S}_d^{*}+{\lambda }_4\left({\mu }_d+{\epsilon }_d+c_d\right)+{\lambda }_6{\beta }_{hF}{S}_h^{*}-{\lambda }_2{\beta }_{dF}{S}_d^{*}-{\lambda }_7{\beta }_{hF}{S}_h^{*}-A_{4 }$$

$$\frac{ d{\lambda }_5}{dt}={\lambda }_5\left({{\mu }_d+\alpha }_d\right)-{\lambda }_1{\alpha }_d$$

$$\frac{d{\lambda }_6}{dt}={\lambda }_6\left({\beta }_{hP}{I}_{dP}^{*}+{\beta }_{hF}{I}_{dF}^{*}\right)+{\lambda }_6\left({\mu }_h+u_3\right)-{\lambda }_7\left({\beta }_{hP}{I}_{dP}^{*}+{\beta }_{hF}{I}_{dF}^{*}\right)-{\lambda }_9{u}_3$$

$$\frac{d{\lambda }_7}{dt}={\lambda }_7\left({{\mu }_h+\delta }_h+u_4\right)-{\lambda }_8{\delta }_h-{\lambda }_9{u}_4-A_2$$

$$\frac{d{\lambda }_8}{dt}={\lambda }_8\left({\mu }_h+{\epsilon }_h\right)-A_5$$

$$\frac{d{\lambda }_9}{dt}=-\frac{\partial H}{\partial R_h^{*}}=-{\lambda }_6{\alpha }_h+{\lambda }_9\left({{\mu }_h+\alpha }_h\right)$$

With transversely conditions, ${\lambda }_i\left(T\right)=0$, $i=1, 2, 3, 4, 5, 6, 7, 8$, 9 i.e.,

$${\lambda }_1\left(T\right)=0, {\lambda }_2\left(T\right)=0, {\lambda }_3\left(T\right)=0, { \lambda }_4\left(T\right)=0, {\lambda }_5\left(T\right)=0, {\lambda }_6\left(T\right)=0,$$

$${\lambda }_7\left(T\right)=0, {\lambda }_8\left(T\right)=0, {\lambda }_9\left(T\right)=0,4.5$$

And

$$u_1^{*}=min\left\{ max\left\{0, \frac{({\lambda }_1-{\lambda }_5)S_d^{*}}{B_1}\right\}, 1\right\}$$

$$u_2^{*}=min\left\{ max\left\{0, \frac{({\lambda }_2-{\lambda }_5)E_d^{*}}{B_2}\right\}, 1\right\}$$

$$u_3^{*}=min\left\{ max\left\{0, \frac{{(\lambda }_6-{\lambda }_9{)S}_h^{*}}{B_3}\right\}, 1\right\}$$

$$u_4^{*}=min\left\{ max\left\{0, \frac{({\lambda }_7-{\lambda }_9{)E}_h^{*}}{B_4}\right\}, 1\right\}$$

Proof: 

By considering partial derivatives of the Hamiltonian H given by (4.3) with respect to the corresponding state variables, one can get the adjoint equations governed by the non-autonomous system (4.4).

$$\frac{d{\lambda }_1}{dt}=-\frac{\partial H}{\partial S_d^{*}}, \frac{d{\lambda }_2}{dt}=-\frac{\partial H}{\partial E_d^{*}}, \frac{d{\lambda }_3}{dt}=-\frac{\partial H}{\partial I_{dP}^{*}}, \frac{d{\lambda }_4}{dt}=-\frac{\partial H}{\partial I_{dF}^{*}}, \frac{d{\lambda }_5}{dt}=-\frac{\partial H}{\partial R_d^{*}}$$

$$\frac{d{\lambda }_6}{dt}=-\frac{\partial H}{\partial S_h^{*}}, \frac{d{\lambda }_7}{dt}=-\frac{\partial H}{\partial E_h^{*}} , \frac{d{\lambda }_8}{dt}=-\frac{\partial H}{\partial I_h^{*}} ,\frac{d{\lambda }_9}{dt}=-\frac{\partial H}{\partial R_h^{*}}$$

$$\frac{d{\lambda }_1}{dt}=-\frac{\partial H}{\partial S_d^{*}}={\lambda }_1\left[\left({\beta }_{dP}{I}_{dP}^{*}+{\beta }_{dF}{I}_{dF}^{*}\right)+\left({\mu }_d+u_1\right)\right]-{\lambda }_2\left({\beta }_{dP}{I}_{dP}^{*}+{\beta }_{dF}{I}_{dF}^{*}\right)-{\lambda }_5{u}_1$$

$$\frac{d{\lambda }_2}{dt}=-\frac{\partial H}{\partial E_d^{*}}={\lambda }_2\left({\mu }_d+{\delta }_d+u_2\right)-{\lambda }_3{\delta }_d-{\lambda }_5{u}_2-A_1$$

$$\frac{d{\lambda }_3}{dt}=-\frac{\partial H}{\partial I_{dP}^{*}}={\lambda }_1{\beta }_{dP}{S}_d^{*}+{\lambda }_3\left({\mu }_d+{\rho }_d+c_d\right)+{\lambda }_6{\beta }_{hP}{S}_h^{*}-{\lambda }_2{\beta }_{dP}{S}_d^{*}$$

$$-{\lambda }_4{\rho }_d-{\lambda }_7{\beta }_{hP}{S}_h^{*}-A_3$$

$$\frac{d{\lambda }_4}{dt}=-\frac{\partial H}{\partial I_{dF}^{*}}={\lambda }_1{\beta }_{dF}{S}_d^{*}+{\lambda }_4\left({\mu }_d+{\epsilon }_d+c_d\right)+{\lambda }_6{\beta }_{hF}{S}_h^{*}-{\lambda }_2{\beta }_{dF}{S}_d^{*}-{\lambda }_7{\beta }_{hF}{S}_h^{*}-A_4$$

$$\frac{d{\lambda }_5}{dt}=-\frac{\partial H}{\partial R_d^{*}}={\lambda }_5\left({{\mu }_d+\alpha }_d\right)-{\lambda }_1{\alpha }_d$$

$$\frac{d{\lambda }_6}{dt}=-\frac{\partial H}{\partial S_h^{*}}={\lambda }_6\left({\beta }_{hP}{I}_{dP}^{*}+{\beta }_{hF}{I}_{dF}^{*}\right)+{\lambda }_6\left({\mu }_h+u_3\right)-{\lambda }_7\left({\beta }_{hP}{I}_{dP}^{*}+{\beta }_{hF}{I}_{dF}^{*}\right)-{\lambda }_9{u}_3$$

$$\frac{d{\lambda }_7}{dt}=-\frac{\partial H}{\partial E_h^{*}}={\lambda }_7\left({{\mu }_h+\delta }_h+u_4\right)-{\lambda }_8{\delta }_h-{\lambda }_9{u}_4-A_2$$

$$\frac{d{\lambda }_8}{dt}=-\frac{\partial H}{\partial I_h^{*}}={\lambda }_8\left({\mu }_h+{\epsilon }_h\right)-A_5$$

$$\frac{d{\lambda }_9}{dt}=-\frac{\partial H}{\partial R_h^{*}}=-{\lambda }_6{\alpha }_h+{\lambda }_9\left({{\mu }_h+\alpha }_h\right)$$

Under terminal conditions or transversality (4.5). Additionally, the following partial differential equations must be solved to determine the optimal control characterization given by (4.6):

$$\frac{\partial H}{\partial u_i}=0, i=1, 2, 3, 4. \mathrm{The\; optimality\; condition\; gives\; us}$$

$$\frac{\partial H}{\partial u_1}{|}_{u_1=u_1^{*}}=B_1{u}_1^{*}-{\lambda }_1{S}_d^{*}+{\lambda }_5{S}_d^{*}=0 \mathrm{i}.\mathrm{e}., u_1^{*}=\frac{\left({\lambda }_1-{\lambda }_5\right)S_d^{*}}{B_1}$$

$$\frac{\partial H}{\partial u_2}{|}_{u_2=u_2^{*}}=B_2{u}_2^{*}-{\lambda }_2{E}_d^{*}+{\lambda }_5{E}_d^{*}=0 \mathrm{i}.\mathrm{e}., u_2^{*}=\frac{({\lambda }_2-{\lambda }_5)E_d^{*}}{B_2}$$

$$\frac{\partial H}{\partial u_3}{|}_{u_3=u_3^{*}}=B_3{u}_3^{*}-{\lambda }_6{S}_h^{*}+{\lambda }_9{S}_h^{*}=0 \mathrm{i}.\mathrm{e}., u_3^{*}=\frac{{(\lambda }_6-{\lambda }_9{)S}_h^{*}}{B_3}$$

$$\frac{\partial H}{\partial u_4}{|}_{u_4=u_4^{*}}=B_4{u}_4^{*}-{\lambda }_7{E}_h^{*}+{\lambda }_9{E}_h^{*}=0 \mathrm{i}.\mathrm{e}., u_4^{*}=\frac{({\lambda }_7-{\lambda }_9{)E}_h^{*}}{B_4}$$

By standard control arguments involving bounds on the control, then

$$u_i^{*}=\left\{\begin{array}{c} 0, if {\phi }_i^{*} \le 0 \\ {\phi }_i^{*} if 0\le {\phi }_i^{*} \le 1\\ 1 if {\phi }_i^{*} \ge 1 \end{array}\right.$$

for $i=1, 2, 3, 4$ and where

$${\phi }_1^{*}=\frac{\left({\lambda }_1-{\lambda }_5\right)S_d^{*}}{B_1}$$

$${\phi }_2^{*}=\frac{({\lambda }_2-{\lambda }_5)E_d^{*}}{B_2}$$

$${\phi }_3^{*}=\frac{{(\lambda }_6-{\lambda }_9{)S}_h^{*}}{B_3}$$

$${\phi }_4^{*}=\frac{({\lambda }_7-{\lambda }_9{)E}_h^{*}}{B_4}$$

The proof is now complete.

6. Simulations and Cost-Effectiveness Analysis

The numerical simulations of the optimality system and the cost-effectiveness analysis of some of the control strategy combinations under consideration are the primary focus of this part.

6.1. Numerical Simulations

The state system (1) together with the adjoint equations (4.4) containing the initial conditions at t = 0, the terminal conditions (4.5), and the optimal control’s characterization (4.6) form the optimality system. Using an iterative approach and a fourth-order forward-backward Runge-Kutta scheme, this optimality system is solved. with the initial conditions at t = 0 for the controls throughout the simulated time, the state system (1) is solved forward in time. Because of the terminal conditions (4.5), the adjoint system is solved backward in time using the state equations’ current iteration responses. The specifics of the numerical process for solving this kind of optimality system with various time orientations are provided [19]. The parameter values from Table 2 are used to show how different combinations of the optimal control intervention options influence the population’s risk of developing rabies. with initial conditions : $S_d=1000000,E_d=8000,$ $I_{dP}=150,$ $I_{dF}=250$, $R_d=50000,$ $S_h=\mathrm{5,000,000}$ $E_h=100,$ $I_h=38,$ and $R_h=2500,$. The weight constants values are chosen as $A_1=A_2={ A}_3={ A}_4=A_5=B_1=B_3=1$ and $B_2=B_4=4$

6.1.1. Intervention 1: Optimal Use of Pre-Exposed Prophylaxis and Treating the Exposed Dogs

This intervention method combines the control effort to treat the exposed dogs, with the control effort to increase immunity in susceptible dogs (i.e., $u_1, u_2\ne 0)$. The magnitudes of infected prodromal phase dogs, infected furious phase dogs, and infectious, and exposed humans decrease more rapidly when controls are in use than in the case without controls, as Figure 6(a-d) shows.

Figure 6. Simulations of the Optimal Control showing the effect of control Strategy 1.

Figure 6. Simulations of the Optimal Control showing the effect of control Strategy 1.

6.1.2. Intervention 2: Optimal Use of Pre-Exposed Prophylaxis Dogs and Pre-Exposed Prophylaxis Humans

Figure 7 shows this intervention method, which combines the optimal use of control effort to increase susceptible dog’s immunity and control effort to increase susceptible human immunity (i.e., $u_1\ne u_3\ne 0)$. The magnitudes of infected prodromal phase dogs, infectious furious phase dogs, and infectious, and exposed humans decrease more quickly in the presence of controls than in the absence of controls, as Figure 7(a-d) shows.

Figure 7. Simulations of the Optimal Control showing the effect of control Strategy 2.

Figure 7. Simulations of the Optimal Control showing the effect of control Strategy 2.

6.1.3. Intervention 3: Optimal Use of Dogs Per-Exposed Prophylaxis and Human Post-Exposed Prophylaxis

This strategy illustrates combines the optimal use of the control effort to increase the immunity of susceptible dogs (per-exposed prophylaxis for dogs) and the control effort at treating the exposed humans (post-exposed prophylaxis for humans) i.e., $u_1, u_4\ne 0$. As expected, the numbers of infectious prodromal phase dogs (Figure 8a), infectious furious phase dogs (Figure 8b), and infectious and Exposed humans (Figure 8c & 8d) diminish more rapidly with controls than the case without controls.

6.1.4. Intervention 4: Optimal Use of Dogs Per-Exposed Prophylaxis, Human’s Per-Exposed Prophylaxis and Human’s Post-Exposed Prophylaxis

This strategy shows the control effort to increase the immunity of susceptible dogs, the control effort increasing the immunity of susceptible humans, and the control effort treating the exposed humans $(u_1, u_3, u_4\ne 0)$. on the spread of rabies dynamics in the population. Figure 9(a) and 9(b) display that the controls decrease the size of infectious in prodromal and Furious phase dogs more rapidly than when controls are not applied. Similarly, the numbers of infectious humans (Figure 9c) and infectious-exposed humans (Figure 9d) decreased more rapidly with controls than the case without controls.

6.1.5. Intervention 4: Implementing All Optimal Control Strategies

This uses the four controls $(u_1\left(t\right), { u}_2\left(t\right), u_3\left(t\right),$and $u_4\left(t\right)$ to minimize the spread of rabies governed by the model (1). The sizes of infectious in prodromal (Figure 9a), Furious phase dogs (Figure 9b), infectious humans (Figure 9c), and exposed humans (Figure 9d), rapidly decrease when controls are in use than the case when controls are not used intervention strategy.

7. Cost-Effectiveness Analysis

We can determine the most effective and least expensive method based on the optimality system’s simulation result using the parameter values listed in Table 2. We used a technique known as the incremental cost-effectiveness ratio (ICER) to arrive at this strategy.

$$ICER=\frac{Change in total cost}{ change in control benefits}$$

By using this method, we were able to compare exceed techniques with more than unity and evaluate which intervention was more effective than the next less effective one. The ratio of the difference in the total number of infections averted to the difference in the averted costs between the two methods was the definition of this technique [19]. We computed the overall cost saved and the total number of infections saved based on the simulation result of the optimal control problem. Based on the total number of infections averted, the control measures are arranged in increasing order in Table 4. The total number of infections prevented was calculated as the difference between the total number of people infected with rabies with controls and the total number of people infected with rabies without controls. The cost function represented by ${\frac{1}{2}B}_1{u}_{ 1}^2, {\frac{1}{2}{B}_{2}u}_{ 2}^2, {{\frac{1}{2}B}_{3}u}_{ 3}^2, and {{\frac{1}{2}B}_{4}u}_{ 4}^2$over was used to determine the total number of infections prevented [19]. Table 3 shows the overall number of diseases prevented and the total cost of all interventions. However, as a single control is unsuccessful in completely eradicating the rabies virus from the human population, we did not take into consideration a strategy that applies only one single control.

we implement the following combinations of the controls.

Strategy 1: The control effort to increase the immunity of susceptible dogs and the control effort aimed at treating the exposed dogs (i.e., $u_1\ne u_2\ne 0)$.

Strategy 2: The control effort to increase the immunity of susceptible dogs and the control effort aimed at increasing the immunity of susceptible humans (i.e., $u_1\ne u_3\ne 0)$.

Strategy 3: The control effort to increase the immunity of susceptible dogs and the control effort aimed at treating the exposed humans (i.e., $u_1\ne u_4\ne 0)$.

Strategy 4: The control effort to increase the immunity of susceptible dogs, the control effort aimed at increasing the immunity of susceptible humans and the control effort aimed at treating the exposed humans (i.e., $u_1\ne u_3\ne u_4\ne 0)$.

Strategy 5: Implementing all controls (i.e., $u_1\ne u_2\ne u_3\ne u_4\ne 0)$

$$ICER\left(2\right)=\frac{\mathrm{44,140,000}}{\mathrm{1,046,800}}=42.16660 , ICER\left(3\right)=\frac{-5876000}{832380}=-7.05928$$

$ICER \left(2\right)>ICER \left(3\right)$, as can be shown. This indicates that strategy 2 is less successful and more expensive than approach 3. As a result, strategy 2 is removed from the list of alternative interventions, and strategies 3 and 4’s ICER are recalculated, as seen in Table 7.

$$ICER\left(3\right)=\frac{\mathrm{38,264,000}}{\mathrm{1,879,180}}=20.36207, ICER\left(4\right)=\frac{-1229000}{174060}=-7.06078$$

Given that ICER(3)$>$ ICER(4), approach 3 is dominant over strategy 4, more expensive, and less successful. As a result, strategy 3 is removed from the list of alternative interventions, and strategies 4 and 1’s ICER are revised as shown in Table 8.

$$ICER\left(4\right)=\frac{\mathrm{37,035,000}}{\mathrm{2,053,240}}=18.03734$$

$ICER\left(4\right)>ICER\left(1\right),$ as could have been observed. This indicates that compared to approach 1, strategy 4 is less effective and more expensive. As a result, strategy 4 is removed from the list of alternative interventions, and strategies 1 and 5’s ICER are recalculated, as shown in Table 9.

$ICER\left(1\right)=\frac{\mathrm{8,649,100}}{\mathrm{2,348,380}}=3.68301$, $ICER\left(5\right)=\frac{-\mathrm{1,977,900}}{\mathrm{280,490}}=-7.06078$Table 9 shows that ICER (1) is greater than ICER (5). This suggests that compared to strategy 5, approach 1 is less effective and more expensive. Thus, out of all the techniques examined in this study, intervention strategy 5 is thought to be the most cost-effective (the most economical).

8. Conclusion and Recommendations

8.1. Conclusion

The dynamics of rabies infection in a human population and among dogs are being studied using a nonlinear mathematical model that is proposed in this study. Because the behaviors of the prodromal and furious phases of infection differ, the model assumes that the infected dog’s compartment was divided into two. Based on these different behaviors, the transmission rate of dogs was divided into two for each population, i.e., ${\beta }_{dP},$ ${\beta }_{dF},$ ${\beta }_{hP}$, and ${\beta }_{hF}$. Additionally, to simulate the effect of those parameters, we solved the model for different values of the most significant model parameters. Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 show the results. The parameters $A_d,{\beta }_{dP},{\beta }_{dF},{ \delta }_d,{ \alpha }_d,$and${ \alpha }_d$ have a positive relationship with the effective reproduction number $\left(R_e\right)$ based on a sensitivity analysis of the parameter. On the other hand, the parameters ${\mu }_d, v_d,{\theta }_d{, c}_d,{\rho }_d$, and ${\epsilon }_d$ are negative signs, then these are inversely proportional to the effective number $\left(R_e\right)$. By adding time-dependent controls to the base model, the model is expanded to an optimum control problem, and the optimality conditions are obtained by applying Pontryagin’s Maximum Principle. Lastly, to determine the effectiveness of different control combinations, numerical simulations of the resulting control problem are performed. The simulation of the control problem shows that the most effective strategy for controlling the disease is the one that makes use of all relevant control measures. Therefore, preventing the spread of rabies requires a multipronged strategy.

8.2. Recommendations

From the results of this study, we suggest the following points to control/minimize/ Dog rabies disease:

• ✓ The dynamics of rabies transmission should be investigated more thoroughly to develop the optimal control measures.
• ✓ Workshops, seminars, and trainings are examples of educational initiatives that should be conducted to raise awareness regarding rabies transmission and prevention strategies.
• ✓ Dog owners should be encouraged by media coverage to crate their pets rather than allow them free reign on the streets.
• ✓ The government and policymakers should come up with strategies to control and reduce the number of stray dogs. A plan that confines owned dogs to specific locations to restrict the annual crop of newborn puppies and lower the frequency of dog-to-dog transmission.
• ✓ The primary focus of strategy development should be on implementing the most effective control methods, such as vaccinating dog populations, decreasing the annual number of newborn puppies, and removing stray dogs, to reduce the human population and rapidly eliminate diseases.