Cost Effectiveness Analysis of Optimal Rabies Control Strategies

Preprint

Article

Cost Effectiveness Analysis of Optimal Rabies Control Strategies

Altmetrics

Downloads

170

Views

Comments

Demsis Dejene Hailemichael^*

Geremew Kenassa Edessa,Purnachandra Rao Koya

Demsis Dejene Hailemichael^*

Geremew Kenassa Edessa,Purnachandra Rao Koya

This version is not peer-reviewed

Submitted:

27 March 2024

Posted:

28 March 2024

You are already at the latest version

Alerts

Abstract

A mathematical model that incorporates human attacks on dogs is presented to explore the dynamics of rabies transmission. The model divided the infection rate into two categories: dog-to-dog transmission rates during the prodromal phase (βdP) and the furious phase (βdF). It has been determined that the model is well-posed and that all of the solutions are positive, as well as the feasibility and positivity of the model's solutions. Both the basic reproduction number (R0) and the effective reproductive number (Re) are computed, and it is demonstrated that the model has a unique disease-free equilibrium that is globally asymptotically stable whenever Re < 1. To identify the model's most sensitive parameters and the ones that intervention efforts should focus on, sensitivity analysis of the effective reproduction number is carried out. We can construct an optimal control system using the state equations, the adjoint equations, and the optimality condition that determines the controls by applying Pontryagin's Maximum/Minimum principle. To help decision-makers in allocating funds for rabies interventions, cost-effectiveness analyses are conducted. To determine the degree to which the intervention strategies are advantageous and cost-effective, this study performs a cost-effective analysis of one or all possible combinations of the optimal rabies control strategies (pre-exposed and post-exposed vaccination for both dog and human populations) for the four different transmission settings. After arranging the techniques to increase effectiveness (total infections prevented), a cost-effective analysis using the Incremental Cost Effectiveness Ratio (ICER) was conducted. The result supports the function that the four intervention options are performing to either completely eradicate or significantly reduce rabies disease among the populations. The dynamical behavior of the system is studied through simulations using the ode45 numerical technique from MATLAB.

Keywords:

Subject: Computer Science and Mathematics - Mathematical and Computational Biology

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 (

β_{d P}

) and the furious phase (

β_{d F}

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_{d P} I_{d F} 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_{d P}

(iv) infected dogs in the furious phase

I_{d F}

(v) recovered dogs

R_{d}

due to pre-exposure vaccination. This

S_{d} E_{d} I_{d P} I_{d F} 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} a n d 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.

\begin{matrix} \frac{d S_{d}}{d t} = A_{d} - {({β_{d P} I}_{d P} + {β_{d F} I}_{d F}) S}_{d} - (μ_{d} + v_{d}) S_{d} {+ α}_{d} R_{d} \\ \frac{d E_{d}}{d t} = {(β_{d P} I_{d p} + {β_{d F} I}_{d F}) S}_{d} - (μ_{d} + δ_{d} + θ_{d}) E_{d} \\ \frac{d I_{d p}}{d t} = δ_{d} E_{d} - (μ_{d} + ρ_{d} + c_{d}) I_{d P} \\ \frac{d I_{d F}}{d t} = ρ_{d} I_{d P} - (μ_{d} + ε_{d} + c_{d}) I_{d F} \\ \frac{d R_{d}}{d t} = v_{d} S_{d} + θ_{d} E_{d} - ({μ_{d} + α}_{d}) R_{d} \\ \frac{d S_{h}}{d t} = A_{h} - {({β_{h P} I}_{d P} + {β_{h F} I}_{d F}) S}_{h} - (μ_{h} + v_{h}) S_{h} {+ α}_{h} R_{h} \\ \frac{d E_{h}}{d t} = {({β_{h P} I}_{d P} + {β_{h F} I}_{d F}) S}_{h} - ({μ_{h} + δ}_{h} + θ_{h}) E_{h} \\ \frac{d I_{h}}{d t} = δ_{h} E_{h} - (μ_{h} + ε_{h}) I_{h} \\ \frac{d R_{h}}{d t} = v_{h} S_{h} + θ_{h} E_{h} - ({μ_{h} + α}_{h}) R_{h} \end{matrix}\}

(1)

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

Variables and Parameters

Description

S_{d}

E_{d}

I_{d p}

I_{d F}

R_{d}

S_{h}

E_{h}

I_{h}

R_{h}

A_{d}

β_{d P}

β_{d F}

δ_{d}

ρ_{d}

v_{d}

θ_{d}

ε_{d}

μ_{d}

α_{d}

A_{h}

β_{h P}

β_{h F}

δ_{h}

v_{h}

θ_{h}

ε_{h}

μ_{h}

α_{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{matrix} 0, & \infty) \end{matrix}

. 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} (t) > 0, E_{d} (t) \geq 0, I_{d P} (t) \geq 0, I_{d F} (t) \geq 0, R_{d} (t) \geq 0

. That is, if the initial conditions

S_{d} (0), E_{d} (0), I_{d P} (0), I_{d F} (0), R_{d} (0)

are non-negative then so are the variables

S_{d} (t), E_{d} (t), I_{d P} (t), I_{d F} (t), R_{d} (t)

for all

t \geq 0

The solutions

\{S_{d} (t), E_{d} (t), I_{d P} (t), I_{d F} (t), R_{d} (t)\}

of the system (1) together with the initial conditions exits and unique on

[\begin{matrix} 0, & k \end{matrix})

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} (t) > 0, E_{d} (t) \geq 0, I_{d P} (t) \geq 0, I_{d F} (t) \geq 0, R_{d} (t) \geq 0

for all

t \geq 0

Positivity of susceptible dog population: Consider the model equation for susceptible dog population

d S_{d} / d t = A_{d} - {(β_{d P} I_{d p} + β_{d F} I_{d F}) S}_{d} - (μ_{d} + v_{d}) S_{d} {+ α}_{d} R_{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} / d t \geq - {(β_{d P} I_{d p} + β_{d F} I_{d F}) S}_{d} - (μ_{d} + v_{d}) S_{d}

.This inequality can also be expressed as

d S_{d} / S_{d} \geq π_{1} d t .

Here, the function

π_{1}

denoting the expression

π_{1} = - (β_{d P} I_{d p} + β_{d F} I_{d F}) - (μ_{d} + v_{d})

can be negative zero or positive valued. Now, using the variables separation method and upon integrating, the solution of

d S_{d} / S_{d} \geq π_{1} d t

can be obtained as

S_{d} (t) \geq S_{d} (0) e^{\int π_{1} d t} .

Here, the integral constant

S_{d} (0)

is the initial susceptible dog population and is assumed to be a non-negative quantity

S_{d} (0) \geq 0

. Similarly, an exponential function is always a non-negative quantity irrespective of the value of the exponent i.e.,

e^{\int π_{1} d t} \geq 0

. Hence, we conclude that

S_{d} (t) \geq 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} / d t \geq - (μ_{d} + δ_{d} + θ_{d}) E_{d}

. Alternately, it can be expressed as

d E_{d} / E_{d} \geq π_{2} d t

. Here, the function

π_{2}

denoting the expression

π_{2} = - (μ_{d} + δ_{d} + θ_{d})

can be negative zero or positive valued. Now, using the variables separation method and applying integration, the solution of

d E_{d} / E_{d} \geq π_{2} d t

can be obtained as

E_{d} (t) = E_{d} (0) e^{\int π_{2} d t} .

Here, the integral constant

E_{d} (0)

is the initial population of exposed dogs and is assumed to be a non-negative quantity

E_{d} (0) \geq 0

. Similarly, an exponential function is always a non-negative quantity irrespective of the value of the exponent i.e.,

e^{\int π_{2} d t} \geq 0

. Hence, we conclude that

E_{d} (t) \geq 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_{d p} / d t = δ_{d} E_{d} - (μ_{d} + ρ_{d} + c_{d}) I_{d P}

. Here, the term

δ_{d} E_{d}

is a positive quantity since

δ_{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

δ_{d} E_{d}

, as an inequality as

d I_{d p} / d t \geq - (μ_{d} + ρ_{d} + c_{d}) I_{d P}

. Now, using the variables separation method and integrating, the solution can be obtained as

I_{d p} (t) = I_{d p} (0) e^{- \int (μ_{d} + ρ_{d} + c_{d}) d t} .

Here, the integral constant

I_{d p} (0)

is the initial population of prodromal infectious dogs and is assumed to be a non-negative quantity

I_{d p} (0) \geq 0

. Similarly, the exponential term is also a positive quantity i.e.,

e^{- \int (μ_{d} + ρ_{d} + c_{d}) d t} \geq 0

since the integrand in the exponent consists of positive parameters only. Hence, we conclude that

I_{d p} (t) \geq 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_{d F} / d t = ρ_{d} I_{d P} - (μ_{d} + ε_{d} + c_{d}) I_{d F}

. Here, the term

ρ_{d} I_{d P}

is positive since

ρ d

is a positive parameter and the prodromal infectious dog population

I_{d P}

is already shown positive. Now, the aforementioned ODE can be expressed without loss of generality, after eliminating the positive term

ρ_{d} I_{d P}

, as an inequality as

d I_{d F} / d t \geq - (μ_{d} + ε_{d} + c_{d}) I_{d F}

. Now, using the variables separation method and integrating, the solution can be obtained as

I_{d F} (t) \geq I_{d F} (0) e^{- \int (μ_{d} + ε_{d} + c_{d}) d t} .

Here, the integral constant

I_{d F} (0)

is the initial population of furious infectious dogs and is assumed to be a non-negative quantity

I_{d F} (0) \geq 0

. Similarly, the exponential term is also positive i.e.,

e^{- \int (μ_{d} + ε_{d} + c_{d}) d t} \geq 0

since the integrand in the exponent consists of positive parameters only. Hence, we conclude that

I_{d F} (t) \geq 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} / d t = v_{d} S_{d} + θ_{d} E_{d} - ({μ_{d} + α}_{d}) R_{d}

. Here, the term

s v_{d} S_{d}

and

θ_{d} E_{d}

are positive since

θ_{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}

θ_{d} E_{d},

as an inequality as

d R_{d} / d t \geq - ({μ_{d} + α}_{d}) R_{d}

. Now, using the variables separation method and integrating, the solution can be obtained as

R_{d} (t) \geq R_{d} (0) e^{- \int ({μ_{d} + α}_{d}) d t} .

Here, the integral constant

R_{d} (0)

is the initial population of recovered dogs and is assumed to be a non-negative quantity

R_{d} (0) \geq 0

. Similarly, the exponential term is also positive i.e.,

e^{- \int ({μ_{d} + α}_{d}) d t} \geq 0

since the integrand in the exponent consists of positive parameters only. Hence, we conclude that

R_{d} (t) \geq 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{matrix} 0, & \infty) \end{matrix}

. 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} (t) \geq 0, E_{h} (t) \geq 0, I_{h} (t) \geq 0, R_{h} (t) \geq 0

. That is, If the initial conditions

S_{h} (0), E_{h} (0), I_{h} (0), a n d R_{h} (0)

are non-negative then so are the variables

S_{h} (t), E_{h} (t), I_{h} (t), R_{h} (t)

for all

t \geq 0

The solutions

\{S_{h} (t), E_{h} (t), I_{h} (t), R_{h} (t)\}

of the system (1) together with the initial conditions exits and unique on

[\begin{matrix} 0, & k \end{matrix})

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} (t) \geq 0, E_{h} (t) \geq 0, I_{h} (t) \geq 0, R_{h} (t) \geq 0

for all

t \geq 0

Positivity of susceptible human population: Consider the model equation for susceptible human population

d S_{h} / d t = A_{h} - β_{h} {(I_{d p} + I_{d F}) S}_{h} - {(μ_{h} + v_{h}) S}_{h} {+ α}_{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} / d t \geq - {(β_{h P} I_{d p} + {β_{h F} I}_{d F}) S}_{h} - {(μ_{h} + v_{h}) S}_{h}

. This inequality can also be expressed as

d S_{h} / S_{h} \geq Q_{1} d t .

Here, the function

Q_{1}

denoting the expression

Q_{1} = - (β_{h P} I_{d p} + {β_{h F} I}_{d F}) - (μ_{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} (t) \geq S_{h} (0) e^{\int Q_{1} d t} .

Here, the integral constant

S_{h} (0)

is the initial susceptible human population and is assumed to be a non-negative quantity

S_{h} (0) \geq 0

. Similarly, an exponential function is always a non-negative quantity irrespective of the value of the exponent i.e.,

e^{\int Q_{1} d t} \geq 0

. Hence, we conclude that

S_{h} (t) \geq 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} / d t \geq - ({μ_{h} + δ}_{h} + θ_{h}) E_{h}

Alternately, it can be expressed as

d E_{h} / E_{h} \geq Q_{2} d t

. Here, the function

Q_{2}

denoting the expression

Q_{2} = - ({μ_{h} + δ}_{h} + θ_{h})

can be negative zero or positive valued. Now, using the variables separation method and applying integration, the solution can be obtained as

E_{h} (t) = E_{h} (0) e^{\int Q_{2} d t} .

Here, the integral constant

E_{h} (0)

is the initial population of exposed humans and is assumed to be a non-negative quantity

E_{h} (0) \geq 0

. Similarly, an exponential function is always a non-negative quantity irrespective of the value of the exponent i.e.,

e^{\int Q_{2} d t} \geq 0

. Hence, we conclude that

E_{h} (t) \geq 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} / d t = δ_{h} E_{h} - (μ_{h} + ε_{h}) I_{h}

. Here, the term

δ_{h} E_{h}

is a positive quantity since

δ_{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

δ_{h} E_{h}

, as an inequality as

d I_{h} / d t \geq - (μ_{h} + ε_{h}) I_{h}

. Now, using the variables separation method and integrating, the solution can be obtained as

I_{h} (t) \geq I_{h} (0) e^{- \int (μ_{h} + ε_{h}) d t} .

Here, the integral constant

I_{h} (0)

is the initial population of infectious humans and is assumed to be a non-negative quantity

I_{h} (0) \geq 0

. Similarly, the exponential term is also a positive quantity i.e.,

e^{- \int (μ_{h} + ε_{h}) d t} \geq 0

since the integrand in the exponent consists of positive parameters only. Hence, we conclude that

I_{h} (t) \geq 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} / d t = v_{h} S_{h} + θ_{h} E_{h} - ({μ_{h} + α}_{h}) R_{h}

. Here, the terms

v_{h} S_{h}

and

θ_{h} E_{h}

are positive since

v_{h}, a n d θ 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}, θ_{h} E_{h}

as an inequality as

d R_{h} / d t \geq - ({μ_{h} + α}_{h}) R_{h}

. Now, using the variables separation method and integrating, the solution can be obtained as

R_{h} (t) \geq R_{h} (0) e^{- \int ({μ_{h} + α}_{h}) d t} .

Here, the integral constant

R_{h} (0)

is the initial population of recovered humans and is assumed to be a non-negative quantity i.e.,

R_{h} (0) \geq 0

. Similarly, the exponential term is also positive i.e.,

e^{- \int ({μ_{h} + α}_{h}) d t} \geq 0

since the integrand in the exponent consists of positive parameters only. Hence, we conclude that

R_{h} (t) \geq 0

. That is, the recovered human population size is a positive quantity.

Theorem 3:

The feasible region

Ω

defined by

Ω = Ω_{d} x Ω_{h} \subset R_{+}^{5} x R_{+}^{4}

is bounded.

Here,

Ω_{d} = \{(\begin{matrix} S_{d} (t), & E_{d} (t), & I_{d P} (t), & I_{d F} (t), & R_{d} (t) \end{matrix}) \in R_{+}^{5} : N_{d} (t) \leq \frac{A_{d}}{μ_{d}}\} Ω_{h} = \{(\begin{matrix} S_{h} (t), & E_{h} (t), & I_{h} (t), & R_{h} (t) \end{matrix}) \in R_{+}^{4} : N_{h} (t) \leq \frac{A_{h}}{μ_{h}}\}

Also,

N_{d} (t) = S_{d} (t) + E_{d} (t) + I_{d P} (t) + I_{d F} (t) + R_{d} (t)

is the total dog population.

Similarly,

N_{h} (t) = S_{h} (t) + E_{h} (t) + I_{h} (t) + R_{h} (t)

is the total human population.

Furthermore, the sets

Ω, Ω_{d},

and

Ω_{h}

are all real-valued.

Proof:

Here, we show that the model population is bounded.

That is, the dog population is bounded i.e.,

N_{d} (t) \leq \frac{A_{d}}{μ_{d}}

Also, the human population is bounded i.e.,

N_{h} (t) \leq \frac{A_{h}}{μ_{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} (t) = S_{d} (t) + E_{d} (t) + I_{d P} (t) + I_{d F} (t) + R_{d} (t)

then

\lim_{t \to \infty} N_{d} (t) \leq \frac{A_{d}}{μ_{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_{d P} + I_{d F} + R_{d}

(3.1)

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

\frac{d N_{d}}{d t} = \frac{d S_{d}}{d t} + \frac{d E_{d}}{d t} + \frac{d I_{d P}}{d t} + \frac{d I_{d F}}{d t} + \frac{d R_{d}}{d 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{d N_{d}}{d t} = A}_{d} - μ_{d} N_{d} - c_{d} I_{d P} - (c_{d} + ε_{d}) I_{d F}

(3.3)

Here, the terms

c_{d} I_{d P}

and

(c_{d} + ε_{d}) I_{d 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{d N_{d}}{d t} + μ_{d} N_{d} \leq A}_{d}

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

N_{d} \leq \frac{A_{d}}{μ_{d}} + [N_{d} (0) - \frac{A_{d}}{μ_{d}}] e^{{- μ}_{d} t}

(3.4)

t \to \infty

in equation (3.4) the population size

N_{d} \to \frac{A_{d}}{μ_{d}}

which implies that

0 {\leq N}_{d} \leq \frac{A_{d}}{μ_{d}}

, Thus the feasible solution set of the system equation of the model enters and remains in the region:

Ω_{d} = \{S_{d}, E_{d}, I_{d P}, {I_{d F,} R}_{d} : N_{d} \leq \frac{A_{d}}{μ_{d}}\} \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} (t) = S_{h} (t) + E_{h} (t) + I_{h} (t) + R_{h} (t)

then

\lim_{t \to \infty} N_{h} (t) \leq \frac{A_{h}}{μ_{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{d N_{h}}{d t} = \frac{d S_{h}}{d t} + \frac{d E_{h}}{d t} + \frac{d I_{h}}{d t} + \frac{d R_{h}}{d 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{d N_{h}}{d t} = A}_{h} - μ_{h} N_{h} - ε_{h} I_{h}

(3.8)

Here, the term

ε_{h} I_{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{d N_{h}}{d t} + μ_{h} N_{h} \leq A}_{h}

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

N_{d} \leq \frac{A_{h}}{μ_{h}} + [N_{h} (0) - \frac{A_{h}}{μ_{h}}] e^{{- μ}_{h} t}

(3.9)

t \to \infty

in equation (3.9) the population size

N_{h} \to \frac{A_{h}}{μ_{h}}

which implies that

0 {\leq N}_{h} \leq \frac{A_{h}}{μ_{h}}

, Thus the feasible solution set of the system equation of the model enters and remains in the region:

Ω_{h} = \{S_{h}, E_{h}, I_{h}, R_{h} : N_{h} \leq \frac{A_{h}}{μ_{h}}\} \in R_{+}^{4}

(3.10)

Therefore, from equations (3.5) and (3.10) the region

Ω = Ω_{+}^{5} \times Ω_{+}^{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 $E_{0}$

To determine the disease-free equilibrium,

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_{d P} = I_{d F} = E_{h} = I_{h} = 0

, and the control parameters

v_{d},

θ_{d},

v_{h},

θ_{h}

are zero. We set the right-hand side of system (1) to zero, leading to the following outcome:

A_{d} - {μ_{d} S}_{d} = 0

A_{h} - μ_{h} S_{h} = 0

E_{0} = {(s}_{d}^{0}, s_{h}^{0}) = (\frac{A_{d}}{μ_{d}}, \frac{A_{d}}{μ_{d}})

(3.11)

To determine the disease-free equilibrium for the effective model is similar to the basic model but considering the parameters

v_{d},

θ_{d},

v_{h},

θ_{h}

are different from zero

\begin{matrix} A_{d} - {({β_{d P} I}_{d P} + {β_{d F} I}_{d F}) S}_{d} - (μ_{d} + v_{d}) S_{d} {+ α}_{d} R_{d} = 0 \\ v_{d} S_{d} + θ_{d} E_{d} - ({μ_{d} + α}_{d}) R_{d} = 0 \\ A_{h} - {({β_{h P} I}_{d P} + {β_{h F} I}_{d F}) S}_{h} - (μ_{h} + v_{h}) S_{h} {+ α}_{h} R_{h} = 0 \\ v_{h} S_{h} + θ_{h} E_{h} - ({μ_{h} + α}_{h}) R_{h} = 0 \end{matrix}\}

(3.12)

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

E_{e} = ({S_{d}^{e}, E}_{d}^{e}, I_{d P}^{e}, I_{d F}^{e}, R_{d}^{e}, S_{h}^{e}, E_{h}^{e}, I_{h}^{e}, R_{h}^{e}) .

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

E_{e} = (\frac{A_{d} ({μ_{d} + α}_{d})}{μ_{d} ({μ_{d} + α}_{d} + v_{d})}, 0, 0, 0, \frac{A_{d} v_{d}}{μ_{d} ({μ_{d} + α}_{d} + v_{d})}, \frac{A_{h} ({μ_{h} + α}_{h})}{μ_{h} ({μ_{h} + α}_{h} + v_{h})}, 0, 0, \frac{A_{h} v_{h}}{μ_{h} ({μ_{h} + α}_{h} + v_{h})})

(3.13)

3.5. Basic Reproduction Number $(R_{0})$

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_{d P}, I_{d F}, E_{h},

and

I_{h}

\begin{matrix} \frac{d E_{d}}{d t} = {(β_{d P} I_{d p} + {β_{d F} I}_{d F}) S}_{d} - (μ_{d} + δ_{d}) E_{d} \\ \frac{d I_{d p}}{d t} = δ_{d} E_{d} - (μ_{d} + ρ_{d}) I_{d P} \\ \frac{d I_{d F}}{d t} = ρ_{d} I_{d P} - (μ_{d} + ε_{d}) I_{d F} \\ \frac{d E_{h}}{d t} = {({β_{h P} I}_{d P} + {β_{h F} I}_{d F}) S}_{h} - ({μ_{h} + δ}_{h}) E_{h} \\ \frac{d I_{h}}{d t} = δ_{h} E_{h} - (μ_{h} + ε_{h}) I_{h} \end{matrix}

(3.14)

From equation (3.14) take the vectors

f_{i}

and

V_{i}

f_{i} = (\begin{matrix} β_{d} {(I_{d p} + I_{d F}) S}_{d} \\ 0 \\ \begin{matrix} 0 \\ β_{h} {(I_{d p} + I_{d F}) S}_{h} \\ 0 \end{matrix} \end{matrix}) & V_{i} = (\begin{matrix} (μ_{d} + δ_{d}) E_{d} \\ (μ_{d} + ρ_{d}) I_{d P} - δ_{d} E_{d} \\ \begin{matrix} (μ_{d} + ε_{d}) I_{d F} - ρ_{d} I_{d P} \\ \begin{matrix} ({μ_{h} + δ}_{h}) E_{h} \\ (μ_{h} + ε_{h}) I_{h} - δ_{h} E_{h} \end{matrix} \end{matrix} \end{matrix})

The Jacobian matrices

f_{i}

and

V_{i}

at disease-free equilibrium from equation (3.11), respectively, are given by

F = (\begin{matrix} \begin{matrix} 0 \\ 0 \\ \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} \end{matrix} & \begin{matrix} \frac{{β_{d P} A}_{d}}{μ_{d}} \\ 0 \\ \begin{matrix} 0 \\ \frac{{β_{h P} A}_{h}}{μ_{h}} \\ 0 \end{matrix} \end{matrix} & \begin{matrix} \begin{matrix} \frac{β_{d F} A_{d}}{μ_{d}} \\ 0 \\ \begin{matrix} 0 \\ \frac{{β_{h F} A}_{h}}{μ_{h}} \\ 0 \end{matrix} \end{matrix} & \begin{matrix} 0 \\ 0 \\ \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} \end{matrix} & \begin{matrix} 0 \\ 0 \\ \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} \end{matrix} \end{matrix} \end{matrix})

(3.15)

V = (\begin{matrix} \begin{matrix} (μ_{d} + δ_{d}) \\ - δ_{d} \\ \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} \end{matrix} & \begin{matrix} 0 \\ (μ_{d} + ρ_{d}) \\ \begin{matrix} - ρ_{d} \\ 0 \\ 0 \end{matrix} \end{matrix} & \begin{matrix} \begin{matrix} 0 \\ 0 \\ \begin{matrix} (μ_{d} + ε_{d}) \\ 0 \\ 0 \end{matrix} \end{matrix} & \begin{matrix} 0 \\ 0 \\ \begin{matrix} 0 \\ ({μ_{h} + δ}_{h}) \\ - δ_{h} \end{matrix} \end{matrix} & \begin{matrix} 0 \\ 0 \\ \begin{matrix} 0 \\ 0 \\ (μ_{h} + ε_{h}) \end{matrix} \end{matrix} \end{matrix} \end{matrix})

Given matrix V, its inverse is given by

V^{- 1} = (\begin{matrix} \frac{1}{(μ_{d} + δ_{d})} & 0 & 0 & 0 & 0 \\ \frac{δ_{d}}{(μ_{d} + δ_{d}) (μ_{d} + ρ_{d})} & \frac{1}{(μ_{d} + ρ_{d})} & 0 & 0 & 0 \\ \frac{ρ_{d} δ_{d}}{(μ_{d} + ρ_{d}) (μ_{d} + δ_{d}) (μ_{d} + ε_{d})} & \frac{ρ_{d}}{(μ_{d} + ρ_{d}) (μ_{d} + ε_{d})} & \frac{1}{(μ_{d} + ε_{d})} & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{({μ_{h} + δ}_{h})} & 0 \\ 0 & 0 & 0 & \frac{δ_{h}}{({μ_{h} + δ}_{h}) (μ_{d} + ε_{d})} & \frac{1}{(μ_{d} + ε_{d})} \end{matrix})

(3.16)

Take F and

V^{- 1}

from equations (3.15) and (3.16), the symbol

ρ (F V^{- 1})

represents the dominant eigenvalue in magnitude, or spectral radius, of a matrix

F V^{- 1}

, thus basic reproductive number defended by

R_{0} = ρ (F V^{- 1})

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

R_{0} = \frac{A_{d} δ_{d} (β_{d P} (μ_{d} + ε_{d}) + ρ_{d} β_{d F})}{μ_{d} (μ_{d} + ρ_{d}) (μ_{d} + ε_{d}) (μ_{d} + δ_{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

R_{0}

but taking into consideration vaccines for dogs and humans both before and after exposure, and culling of infected dogs i.e.,

θ_{d}, θ_{h}, v_{d}, v_{h} c_{d} \neq 0

, the effective reproduction number

R_{e}

of system (1) is calculated.

R_{e} = ρ (F V^{- 1})

is the spectral radius, or dominant eigenvalue, of

F V^{- 1} .

The Jacobian matrices F and V at effective disease-free equilibrium

E_{e}

, respectively, are given by

F = (\begin{matrix} \begin{matrix} 0 \\ 0 \\ \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} \end{matrix} & \begin{matrix} \frac{β_{d P} A_{d} ({μ_{d} + α}_{d})}{μ_{d} ({μ_{d} + α}_{d} + v_{d})} \\ 0 \\ \begin{matrix} 0 \\ \frac{β_{h P} A_{h} ({μ_{h} + α}_{h})}{μ_{h} ({μ_{h} + α}_{h} + v_{h})} \\ 0 \end{matrix} \end{matrix} & \begin{matrix} \begin{matrix} \frac{β_{d F} A_{d} ({μ_{d} + α}_{d})}{μ_{d} ({μ_{d} + α}_{d} + v_{d})} \\ 0 \\ \begin{matrix} 0 \\ \frac{β_{h F} A_{h} ({μ_{h} + α}_{h})}{μ_{h} ({μ_{h} + α}_{h} + v_{h})} \\ 0 \end{matrix} \end{matrix} & \begin{matrix} 0 \\ 0 \\ \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} \end{matrix} & \begin{matrix} 0 \\ 0 \\ \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} \end{matrix} \end{matrix} \end{matrix})

(3.18)

V = (\begin{matrix} \begin{matrix} (μ_{d} + δ_{d} + θ_{d}) \\ - δ_{d} \\ \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} \end{matrix} & \begin{matrix} 0 \\ (μ_{d} + ρ_{d} + c_{d}) \\ \begin{matrix} - ρ_{d} \\ 0 \\ 0 \end{matrix} \end{matrix} & \begin{matrix} \begin{matrix} 0 \\ 0 \\ \begin{matrix} (μ_{d} + ε_{d} + c_{d}) \\ 0 \\ 0 \end{matrix} \end{matrix} & \begin{matrix} 0 \\ 0 \\ \begin{matrix} 0 \\ ({μ_{h} + δ}_{h} + θ_{h}) \\ - δ_{h} \end{matrix} \end{matrix} & \begin{matrix} 0 \\ 0 \\ \begin{matrix} 0 \\ 0 \\ (μ_{h} + ε_{h}) \end{matrix} \end{matrix} \end{matrix} \end{matrix})

The inverse of matrix V is given by

V^{- 1} = (\begin{matrix} \begin{matrix} \frac{1}{(μ_{d} + δ_{d} + θ_{d})} \\ \frac{δ_{d}}{(μ_{d} + δ_{d} + θ_{d}) (μ_{d} + ρ_{d} + c_{d})} \\ \begin{matrix} \frac{ρ_{d} δ_{d}}{(μ_{d} + ρ_{d} + c_{d}) (μ_{d} + δ_{d} + θ_{d}) (μ_{d} + ε_{d} + c_{d})} \\ 0 \\ 0 \end{matrix} \end{matrix} & \begin{matrix} 0 \\ \frac{1}{(μ_{d} + ρ_{d} + c_{d})} \\ \begin{matrix} \frac{ρ_{d}}{(μ_{d} + ρ_{d} + c_{d}) (μ_{d} + ε_{d} + c_{d})} \\ 0 \\ 0 \end{matrix} \end{matrix} & \begin{matrix} 0 \\ 0 \\ \begin{matrix} \frac{1}{(μ_{d} + ε_{d} + c_{d})} \\ 0 \\ 0 \end{matrix} \end{matrix} \end{matrix}

\begin{matrix} \begin{matrix} 0 \\ 0 \\ \begin{matrix} 0 \\ \frac{1}{({μ_{h} + δ}_{h} + θ_{h})} \\ \frac{δ_{h}}{({μ_{h} + δ}_{h} + θ_{h}) (μ_{d} + ε_{d})} \end{matrix} \end{matrix} & \begin{matrix} 0 \\ 0 \\ \begin{matrix} 0 \\ 0 \\ \frac{1}{(μ_{d} + ε_{d})} \end{matrix} \end{matrix} \end{matrix})

(3.19)

The symbol

ρ (F V^{- 1})

represents the dominant eigenvalue in magnitude, or spectral radius, of a matrix

F V^{- 1}

, and take

F a n d V^{- 1}

from equation (3.18) and (3.19), thus effective reproductive number defend by

R_{e} = ρ (F V^{- 1})

Therefore, the effective reproduction number is as follows:

R_{e} = \frac{A_{d} (β_{d P} (μ_{d} + ε_{d} + c_{d}) + ρ_{d} β_{d F}) δ_{d} ({μ_{d} + α}_{d})}{μ_{d} ({μ_{d} + v_{d} + α}_{d}) (μ_{d} + ρ_{d} + c_{d}) (μ_{d} + ε_{d} + c_{d}) (μ_{d} + δ_{d} + θ_{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.

R_{0} = \frac{A_{d} δ_{d} (β_{d P} (μ_{d} + ε_{d}) + ρ_{d} β_{d F})}{μ_{d} (μ_{d} + ρ_{d}) (μ_{d} + ε_{d}) (μ_{d} + δ_{d})} = 2.74

The fact that the basic reproductive number

R_{0}

is larger than one in the absence of any control measures indicates that the disease will keep spreading in the population.

R_{e} = \frac{A_{d} (β_{d P} (μ_{d} + ε_{d} + c_{d}) + ρ_{d} β_{d F}) δ_{d} ({μ_{d} + α}_{d})}{μ_{d} ({μ_{d} + v_{d} + α}_{d}) (μ_{d} + ρ_{d} + c_{d}) (μ_{d} + ε_{d} + c_{d}) (μ_{d} + δ_{d} + θ_{d})} = 1.78

The effective reproductive number

R_{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.

\begin{matrix} \frac{d X}{d t} = B_{0} (X - X_{0}) + B_{1} Y \\ \frac{d Y}{d t} = B_{2} Y \end{matrix}\}

(3.21)

The classifications of susceptible and vaccinated individuals are represented by

X \in R_{+}^{m}

. Classes of exposed and infectious individuals are represented by

Y \in R_{+}^{n}

A vector at DFE point

ε_{0}

of the vector length as

X

is represented by Using system (1) as an example, we define:

X = [\begin{matrix} S_{d} \\ R_{d} \\ S_{h} \\ R_{h} \end{matrix}], Y = [\begin{matrix} E_{d} \\ I_{d P} \\ I_{d F} \\ E_{d} \\ I_{h} \end{matrix}] and X_{ε_{0}} = [\begin{matrix} \frac{A_{d} ({μ_{d} + α}_{d})}{μ_{d} ({μ_{d} + α}_{d} + v_{d})} \\ \frac{A_{d} v_{d}}{μ_{d} ({μ_{d} + α}_{d} + v_{d})} \\ \frac{A_{h} ({μ_{h} + α}_{h})}{μ_{h} ({μ_{h} + α}_{h} + v_{h})} \\ \frac{A_{h} v_{h}}{μ_{h} ({μ_{h} + α}_{h} + v_{h})} \end{matrix}], X - X_{0} = [\begin{matrix} S_{d} - \frac{A_{d} ({μ_{d} + α}_{d})}{μ_{d} ({μ_{d} + α}_{d} + v_{d})} \\ R_{d} - \frac{A_{d} v_{d}}{μ_{d} ({μ_{d} + α}_{d} + v_{d})} \\ S_{h} - \frac{A_{h} ({μ_{h} + α}_{h})}{μ_{h} ({μ_{h} + α}_{h} + v_{h})} \\ R_{d} - \frac{A_{h} v_{h}}{μ_{h} ({μ_{h} + α}_{h} + v_{h})} \end{matrix}]

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:

[\begin{matrix} A_{d} {+ α}_{d} R_{d} - {({β_{d P} I}_{d P} + {β_{d F} I}_{d F} + μ_{d} + v_{d}) S}_{d} \\ v_{d} S_{d} + θ_{d} E_{d} - ({μ_{d} + α}_{d}) R_{d} \\ A_{h} {+ α}_{h} R_{h} - {({β_{h P} I}_{d P} + {β_{h F} I}_{d F} + μ_{h} + v_{h}) S}_{h} \\ v_{h} S_{h} + θ_{h} E_{h} - ({μ_{h} + α}_{h}) R_{h} \end{matrix}] = B_{0} [\begin{matrix} S_{d} - \frac{A_{d} ({μ_{d} + α}_{d})}{μ_{d} ({μ_{d} + α}_{d} + v_{d})} \\ R_{d} - \frac{A_{d} v_{d}}{μ_{d} ({μ_{d} + α}_{d} + v_{d})} \\ S_{h} - \frac{A_{h} ({μ_{h} + α}_{h})}{μ_{h} ({μ_{h} + α}_{h} + v_{h})} \\ R_{d} - \frac{A_{h} v_{h}}{μ_{h} ({μ_{h} + α}_{h} + v_{h})} \end{matrix}] + B_{1} [\begin{matrix} E_{d} \\ I_{d P} \\ I_{d F} \\ E_{d} \\ I_{h} \end{matrix}]

(3.22)

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

[\begin{matrix} A_{d} {+ α}_{d} R_{d} - {(+ μ_{d} + v_{d}) S}_{d} \\ v_{d} S_{d} - ({μ_{d} + α}_{d}) R_{d} \\ A_{h} {+ α}_{h} R_{h} - {(+ μ_{h} + v_{h}) S}_{h} \\ v_{h} S_{h} - ({μ_{h} + α}_{h}) R_{h} \end{matrix}] + [\begin{matrix} - {({β_{d P} I}_{d P} + {β_{d F} I}_{d F}) S}_{d} \\ θ_{d} E_{d} \\ - {({β_{h P} I}_{d P} + {β_{h F} I}_{d F}) S}_{h} \\ + θ_{h} E_{h} \end{matrix}] = B_{0} [\begin{matrix} S_{d} - S_{d}^{0} \\ R_{d} - R_{d}^{0} \\ S_{h} - S_{h}^{0} \\ R_{d} - R_{h}^{0} \end{matrix}] + B_{1} [\begin{matrix} E_{d} \\ I_{d P} \\ I_{d F} \\ E_{d} \\ I_{h} \end{matrix}]

(3.23)

Using the state variables

S_{d}, R_{d}, S_{h}, R_{h}

, and

E_{d}, I_{d P}, I_{d F}, a n d E_{h}

of the Jacobian matrix, we obtain the first vector and the second vector on the left side of the equation (3.23).

B_{0} = [\begin{matrix} - (+ μ_{d} + v_{d}) & α_{h} & 0 & 0 \\ v_{d} & - ({μ_{d} + α}_{d}) & 0 & 0 \\ 0 & 0 & - (μ_{h} + v_{h}) & α_{h} \\ 0 & 0 & v_{h} & - ({μ_{h} + α}_{h}) \end{matrix}], B_{1} = [\begin{matrix} 0 & - β_{d P} S_{d} & - β_{d F} S_{d} & 0 \\ θ_{d} & 0 & 0 & 0 \\ 0 & - β_{h P} S_{h} & - β_{h F} S_{h} & 0 \\ 0 & 0 & 0 & θ_{h} \end{matrix}]

(3.23)

Next, show that a matrix

B_{0}

has real, negative eigenvalues.

|B_{0} - I λ| = 0

is the characteristic equation of matrix

B_{0}

|\begin{matrix} - (μ_{d} + v_{d} + λ) & α_{h} & 0 & 0 \\ v_{d} & - ({μ_{d} + α}_{d} + λ) & 0 & 0 \\ 0 & 0 & - (μ_{h} + v_{h} + λ) & α_{h} \\ 0 & 0 & v_{h} & - ({μ_{h} + α}_{h} + λ) \end{matrix}| = 0

- (μ_{d} + v_{d} + λ) |\begin{matrix} - ({μ_{d} + α}_{d} + λ) & 0 & 0 \\ 0 & - (μ_{h} + v_{h} + λ) & α_{h} \\ 0 & v_{h} & - ({μ_{h} + α}_{h} + λ) \end{matrix}| = 0

(μ_{d} + v_{d} + λ) ({μ_{d} + α}_{d} + λ) |\begin{matrix} - (μ_{h} + v_{h} + λ) & α_{h} \\ v_{h} & - ({μ_{h} + α}_{h} + λ) \end{matrix}| = 0

(μ_{d} + v_{d} + λ) ({μ_{d} + α}_{d} + λ) ((μ_{h} + v_{h} + λ) ({μ_{h} + α}_{h} + λ) - v_{h} α_{h}) = 0

(μ_{d} + v_{d} + λ) ({μ_{d} + α}_{d} + λ) (μ_{h} + λ) ((μ_{h} + v_{h} + α_{h}) + λ) = 0

As a result, the following Eigen values exist.

\{\begin{matrix} λ_{1} = - (μ_{d} + v_{d}) \\ λ_{2} = - ({μ_{d} + α}_{d}) \\ λ_{3} = {- μ}_{h} \\ λ_{4} = - (μ_{h} + v_{h} + α_{h}) \end{matrix}

(3.24)

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

[\begin{matrix} {(β_{d P} I_{d p} + {β_{d F} I}_{d F}) S}_{d} - (μ_{d} + δ_{d}) E_{d} \\ δ_{d} E_{d} - (μ_{d} + ρ_{d}) I_{d P} \\ ρ_{d} I_{d P} - (μ_{d} + ε_{d}) I_{d F} \\ {({β_{h P} I}_{d P} + {β_{h F} I}_{d F}) S}_{h} - ({μ_{h} + δ}_{h}) E_{h} \\ δ_{h} E_{h} - (μ_{h} + ε_{h}) I_{h} \end{matrix}] = B_{2} [\begin{matrix} E_{d} \\ I_{d P} \\ I_{d F} \\ E_{d} \\ I_{h} \end{matrix}]

(3.25)

Using the Jacobian matrix’s elements

E_{d}, I_{d P}, I_{d F}, E_{h,}

I_{h}

the left side of equation (3.25) becomes:

B_{2} = [\begin{matrix} - (μ_{d} + δ_{d}) & β_{d P} S_{d} & β_{d F} S_{d} & 0 & 0 \\ δ_{d} & - (μ_{d} + ρ_{d}) & 0 & 0 & 0 \\ 0 & ρ_{d} & - (μ_{d} + ε_{d}) & 0 & 0 \\ 0 & β_{h P} S_{h} & β_{h F} S_{h} & - ({μ_{h} + δ}_{h}) & 0 \\ 0 & 0 & 0 & δ_{h} & - (μ_{h} + ε_{h}) \end{matrix}]

(3.26)

Thus,

B_{2}

is a Metzler matrix

Therefore, from equation (3.24) all Eigenvalues of matrix

B_{0}

are

- (μ_{d} + v_{d}),

- ({μ_{d} + α}_{d}),

{- μ}_{h}

and

- (μ_{h} + v_{h} + α_{h})

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

R_{e} < 1

and unstable in the region Ω. if

R_{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

R_{e} = \frac{A_{d} (β_{d P} (μ_{d} + ε_{d} + c_{d}) + ρ_{d} β_{d F}) δ_{d} ({μ_{d} + α}_{d})}{μ_{d} ({μ_{d} + v_{d} + α}_{d}) (μ_{d} + ρ_{d} + c_{d}) (μ_{d} + ε_{d} + c_{d}) (μ_{d} + δ_{d} + θ_{d})}

. Since

R_{e}

only depends on eleven parameters

A_{d}, β_{d P}, β_{d F}

δ_{d}, ρ_{d}, v_{d}, θ_{d}, c d ε_{d}, μ_{d},

and

α_{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_{ω}^{R_{e}}

is defined as

S_{ω}^{R_{e}} = \frac{\partial R_{e}}{\partial ω} \times \frac{ω}{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:

S_{A_{d}}^{R_{e}} = \frac{\partial R_{e}}{\partial A_{d}} \times \frac{A_{d}}{R_{e}} = + 1

Similarly, the sensitivity index concerning

β_{d P}, {β_{d F}, δ}_{d}, ρ_{d}, v_{d}, θ_{d}, c_{d} ε_{d}, μ_{d},

and

α_{d}

is given by: In the same way, concerning

β_{d P}, {β_{d F}, δ}_{d}, ρ_{d}, v_{d}, θ_{d}, c_{d} ε_{d}, μ_{d},

and

α_{d}

, the sensitivity index is provided by:

S_{β_{d P}}^{R_{e}} = \frac{\partial R_{e}}{\partial β_{d P}} \times \frac{β_{d P}}{R_{e}} = 0.001747639115

S_{β_{d F}}^{R_{e}} = \frac{\partial R_{e}}{\partial β_{d F}} \times \frac{β_{d F}}{R_{e}} = 0.9982523609

S_{δ_{d}}^{R_{e}} = \frac{\partial R_{e}}{\partial δ_{d}} \times \frac{δ_{d}}{R_{e}} = 0.6056304706

S_{α_{d}}^{R_{e}} = \frac{\partial R_{e}}{\partial α_{d}} \times \frac{α_{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_{μ_{d}}^{R_{e}} = \frac{\partial R_{e}}{\partial μ_{d}} \times \frac{μ_{d}}{R_{e}} = - 1.220608197

S_{θ_{d}}^{R_{e}} = \frac{\partial R_{e}}{\partial θ_{d}} \times \frac{θ_{d}}{R_{e}} = - 0.4731488052

S_{ε_{d}}^{R_{e}} = \frac{\partial R_{e}}{\partial ε_{d}} \times \frac{ε_{d}}{R_{e}} = - 0.8081706290

S_{ρ_{d}}^{R_{e}} = \frac{\partial R_{e}}{\partial ρ_{d}} \times \frac{ρ_{d}}{R_{e}} = 0.2209374584

4.1. Interpretation of Sensitivity Indices

We observe the sensitivity index values’ sign. The parameters (

A_{d}, β_{d P}, β_{d F}, δ_{d},

and

α_{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}, β_{d P}, β_{d F}, δ_{d},

and

α_{d}

) directly proportional to

R_{e}

. In addition, the parameters (

v_{d}, θ_{d}, c_{d}, ρ_{d}

, and

ε_{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}, β_{d P}, β_{d F}, δ_{d}, θ_{d}, δ_{d}, ε_{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

β_{d F}

. 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

β_{d F},

and

β_{d P}

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

δ_{d}

. Values of

δ_{d}

were varied as shown in the figure.

Figure 4. Effect of parameter

A_{d}

R_{e}

Figure 4. Effect of parameter

A_{d}

R_{e}

Figure 5. Effect of parameter

δ_{d}

R_{e}

Figure 5. Effect of parameter

δ_{d}

R_{e}

Figure 6. Effect of parameter

θ_{d}

R_{e}

Figure 6. Effect of parameter

θ_{d}

R_{e}

Figure 6 shows an inverse relationship between the effective reproduction number

R_{e}

and the recovery rate

θ_{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

ε_{d}

R_{e}

Figure 7. Effect of parameter

ε_{d}

R_{e}

In Figure 7. The effective reproduction number

R_{e}

is a decreasing function of rabies-induced death

ε_{d}

. Values of

ε_{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_{d P}

, and infectious dogs in the furious phase

I_{d F}

, 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} (t), (i = 1, 2, 3, 4) :

Our objective function becomes

J = \min_{u \in U} \int_{0}^{T} ({A_{1} E}_{d} + {A_{2} E}_{h} + A_{3} I_{d P} + {A_{4} I}_{d F} + {A_{5} I}_{h} + \frac{1}{2} \sum_{i = 1}^{4} B_{i} u_{i}^{2} (t)) d t

(4.1)

subject to the model (1),

\begin{matrix} \frac{d S_{d}}{d t} = A_{d} - {({β_{d P} I}_{d P} + {β_{d F} I}_{d F}) S}_{d} - (μ_{d} + u_{1}) S_{d} {+ α}_{d} R_{d} \\ \frac{d E_{d}}{d t} = {(β_{d P} I_{d p} + {β_{d F} I}_{d F}) S}_{d} - (μ_{d} + δ_{d} + u_{2}) E_{d} \\ \frac{d I_{d p}}{d t} = δ_{d} E_{d} - (μ_{d} + ρ_{d} + c_{d}) I_{d P} \\ \frac{d I_{d F}}{d t} = ρ_{d} I_{d P} - (μ_{d} + ε_{d} + c_{d}) I_{d F} \\ \frac{d R_{d}}{d t} = u_{1} S_{d} + u_{2} E_{d} - ({μ_{d} + α}_{d}) R_{d} \\ \frac{d S_{h}}{d t} = A_{h} - {({β_{h P} I}_{d P} + {β_{h F} I}_{d F}) S}_{h} - (μ_{h} + u_{3}) S_{h} {+ α}_{h} R_{h} \\ \frac{d E_{h}}{d t} = {({β_{h P} I}_{d P} + {β_{h F} I}_{d F}) S}_{h} - ({μ_{h} + δ}_{h} + u_{4}) E_{h} \\ \frac{d I_{h}}{d t} = δ_{h} E_{h} - (μ_{h} + ε_{h}) I_{h} \\ \frac{d R_{h}}{d t} = u_{3} S_{h} + u_{4} E_{h} - ({μ_{h} + α}_{h}) R_{h} \end{matrix}\}

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 (u_{1}^{*}, u_{2}^{*}, u_{3}^{*}, u_{4}^{*}) = \min J {(u_{1}, u_{2}, u_{3}, u_{4}) | u_{i} \in U 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 \leq u_{i} \leq 1

for

i = 1, 2, 3, 4}; L e b e s g u e m e a s u r a b l e \forall_{t} \in [0, T],

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} a n d u_{4}

The Hamiltonian equation with state variables

S_{d} = S_{d}^{*}, E_{d} = E_{d}^{*},

I_{d P} = I_{d P}^{*}

I_{d F} = I_{d F}^{*}

R_{d} = R_{d}^{*}

S_{h} = S_{h}^{*}, E_{h} = E_{h}^{*},

I_{h} = I_{h}^{*},

and

R_{h} = R_{h}^{*}

is formed as

H = I i n t e g r a n d + (a d j o i n s) \times (R H S)

of the

\begin{matrix} H = A_{1} E_{d}^{*} + A_{2} E_{h}^{*} + A_{3} I_{d P}^{*} + A_{4} I_{d F}^{*} + A_{5} I_{h}^{*} \\ + \frac{1}{2} [B_{1} u_{1}^{2} + {B_{2} u}_{2}^{2} + {B_{3} u}_{3}^{2} + {B_{4} u}_{4}^{2}] \\ + λ_{1} [A_{d} - (β_{d P} I_{d P}^{*} + β_{d F} I_{d F}^{*}) S_{d}^{*} - (μ_{d} + u_{1}) S_{d}^{*} {+ α}_{d} R_{d}^{*}] 4.3 \\ + λ_{2} [(β_{d P} I_{d P}^{*} + β_{d F} I_{d F}^{*}) S_{d}^{*} - (μ_{d} + δ_{d} + u_{2}) E_{d}^{*}] \\ + λ_{3} [δ_{d} E_{d}^{*} - (μ_{d} + ρ_{d} + c_{d}) I_{d P}^{*}] \\ + λ_{4} [ρ_{d} I_{d P}^{*} - (μ_{d} + ε_{d} + c_{d}) I_{d F}^{*}] \\ + λ_{5} [u_{1} S_{d}^{*} + u_{2} E_{d}^{*} - ({μ_{d} + α}_{d}) R_{d}^{*}] \\ + λ_{6} [A_{h} - ({β_{h P} I}_{d P}^{*} + β_{h F} I_{d F}^{*}) S_{h}^{*} - {(μ}_{h} + u_{3}) S_{h}^{*} {+ α}_{h} R_{h}^{*}] \\ + λ_{7} [({β_{h P} I}_{d P}^{*} + β_{h F} I_{d F}^{*}) S_{h}^{*} - ({μ_{h} + δ}_{h} + u_{4}) E_{h}^{*}] \\ + λ_{8} [δ_{h} E_{h}^{*} - (μ_{h} + ε_{h}) I_{h}^{*}] \\ + λ_{9} [u_{3} S_{h}^{*} + u_{4} E_{h}^{*} - ({μ_{h} + α}_{h}) R_{h}^{*}] \end{matrix}

where the adjoint (co-state) variables are

λ_{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

u_{1}^{*}, u_{2}^{*}, u_{3}^{*}, u_{4}^{*}

that minimizes objective functional (4.1) over the control set U subject to the state system (1), then there exist adjoint variables

λ_{1}, λ_{2,} λ_{3,} λ_{4}, λ_{5}, λ_{6}, λ_{7}, λ_{8}, λ_{9}

satisfying

\frac{d λ_{1}}{d t} = λ_{1} [(β_{d P} I_{d P}^{*} + β_{d F} I_{d F}^{*}) + (μ_{d} + u_{1})] - λ_{2} (β_{d P} I_{d P}^{*} + β_{d F} I_{d F}^{*}) - λ_{5} u_{1}

\frac{d λ_{2}}{d t} = λ_{2} (μ_{d} + δ_{d} + u_{2}) - λ_{3} δ_{d} - λ_{5} u_{2} - A_{1}

\frac{d λ_{3}}{d t} = λ_{1} β_{d P} S_{d}^{*} + λ_{3} (μ_{d} + ρ_{d} + c_{d}) + λ_{6} β_{h P} S_{h}^{*} - λ_{2} β_{d P} S_{d}^{*}

- λ_{4} ρ_{d} - λ_{7} β_{h P} S_{h}^{*} - A_{3} 4.4

\frac{d λ_{4}}{d t} = λ_{1} β_{d F} S_{d}^{*} + λ_{4} (μ_{d} + ε_{d} + c_{d}) + λ_{6} β_{h F} S_{h}^{*} - λ_{2} β_{d F} S_{d}^{*} - λ_{7} β_{h F} S_{h}^{*} - A_{4}

\frac{d λ_{5}}{d t} = λ_{5} ({μ_{d} + α}_{d}) - λ_{1} α_{d}

\frac{d λ_{6}}{d t} = λ_{6} (β_{h P} I_{d P}^{*} + β_{h F} I_{d F}^{*}) + λ_{6} (μ_{h} + u_{3}) - λ_{7} (β_{h P} I_{d P}^{*} + β_{h F} I_{d F}^{*}) - λ_{9} u_{3}

\frac{d λ_{7}}{d t} = λ_{7} ({μ_{h} + δ}_{h} + u_{4}) - λ_{8} δ_{h} - λ_{9} u_{4} - A_{2}

\frac{d λ_{8}}{d t} = λ_{8} (μ_{h} + ε_{h}) - A_{5}

\frac{d λ_{9}}{d t} = - \frac{\partial H}{\partial R_{h}^{*}} = - λ_{6} α_{h} + λ_{9} ({μ_{h} + α}_{h})

With transversely conditions,

λ_{i} (T) = 0

i = 1, 2, 3, 4, 5, 6, 7, 8

, 9 i.e.,

λ_{1} (T) = 0, λ_{2} (T) = 0, λ_{3} (T) = 0, λ_{4} (T) = 0, λ_{5} (T) = 0, λ_{6} (T) = 0,

λ_{7} (T) = 0, λ_{8} (T) = 0, λ_{9} (T) = 0, 4.5

And

u_{1}^{*} = m i n \{m a x \{0, \frac{(λ_{1} - λ_{5}) S_{d}^{*}}{B_{1}}\}, 1\}

u_{2}^{*} = m i n \{m a x \{0, \frac{(λ_{2} - λ_{5}) E_{d}^{*}}{B_{2}}\}, 1\}

u_{3}^{*} = m i n \{m a x \{0, \frac{{(λ}_{6} - λ_{9} {) S}_{h}^{*}}{B_{3}}\}, 1\}

u_{4}^{*} = m i n \{m a x \{0, \frac{(λ_{7} - λ_{9} {) E}_{h}^{*}}{B_{4}}\}, 1\}

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 λ_{1}}{d t} = - \frac{\partial H}{\partial S_{d}^{*}}, \frac{d λ_{2}}{d t} = - \frac{\partial H}{\partial E_{d}^{*}}, \frac{d λ_{3}}{d t} = - \frac{\partial H}{\partial I_{d P}^{*}}, \frac{d λ_{4}}{d t} = - \frac{\partial H}{\partial I_{d F}^{*}}, \frac{d λ_{5}}{d t} = - \frac{\partial H}{\partial R_{d}^{*}}

\frac{d λ_{6}}{d t} = - \frac{\partial H}{\partial S_{h}^{*}}, \frac{d λ_{7}}{d t} = - \frac{\partial H}{\partial E_{h}^{*}}, \frac{d λ_{8}}{d t} = - \frac{\partial H}{\partial I_{h}^{*}}, \frac{d λ_{9}}{d t} = - \frac{\partial H}{\partial R_{h}^{*}}

\frac{d λ_{1}}{d t} = - \frac{\partial H}{\partial S_{d}^{*}} = λ_{1} [(β_{d P} I_{d P}^{*} + β_{d F} I_{d F}^{*}) + (μ_{d} + u_{1})] - λ_{2} (β_{d P} I_{d P}^{*} + β_{d F} I_{d F}^{*}) - λ_{5} u_{1}

\frac{d λ_{2}}{d t} = - \frac{\partial H}{\partial E_{d}^{*}} = λ_{2} (μ_{d} + δ_{d} + u_{2}) - λ_{3} δ_{d} - λ_{5} u_{2} - A_{1}

\frac{d λ_{3}}{d t} = - \frac{\partial H}{\partial I_{d P}^{*}} = λ_{1} β_{d P} S_{d}^{*} + λ_{3} (μ_{d} + ρ_{d} + c_{d}) + λ_{6} β_{h P} S_{h}^{*} - λ_{2} β_{d P} S_{d}^{*}

- λ_{4} ρ_{d} - λ_{7} β_{h P} S_{h}^{*} - A_{3}

\frac{d λ_{4}}{d t} = - \frac{\partial H}{\partial I_{d F}^{*}} = λ_{1} β_{d F} S_{d}^{*} + λ_{4} (μ_{d} + ε_{d} + c_{d}) + λ_{6} β_{h F} S_{h}^{*} - λ_{2} β_{d F} S_{d}^{*} - λ_{7} β_{h F} S_{h}^{*} - A_{4}

\frac{d λ_{5}}{d t} = - \frac{\partial H}{\partial R_{d}^{*}} = λ_{5} ({μ_{d} + α}_{d}) - λ_{1} α_{d}

\frac{d λ_{6}}{d t} = - \frac{\partial H}{\partial S_{h}^{*}} = λ_{6} (β_{h P} I_{d P}^{*} + β_{h F} I_{d F}^{*}) + λ_{6} (μ_{h} + u_{3}) - λ_{7} (β_{h P} I_{d P}^{*} + β_{h F} I_{d F}^{*}) - λ_{9} u_{3}

\frac{d λ_{7}}{d t} = - \frac{\partial H}{\partial E_{h}^{*}} = λ_{7} ({μ_{h} + δ}_{h} + u_{4}) - λ_{8} δ_{h} - λ_{9} u_{4} - A_{2}

\frac{d λ_{8}}{d t} = - \frac{\partial H}{\partial I_{h}^{*}} = λ_{8} (μ_{h} + ε_{h}) - A_{5}

\frac{d λ_{9}}{d t} = - \frac{\partial H}{\partial R_{h}^{*}} = - λ_{6} α_{h} + λ_{9} ({μ_{h} + α}_{h})

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 . The optimality condition gives us

\frac{\partial H}{\partial u_{1}} |_{u_{1} = u_{1}^{*}} = B_{1} u_{1}^{*} - λ_{1} S_{d}^{*} + λ_{5} S_{d}^{*} = 0 i . e ., u_{1}^{*} = \frac{(λ_{1} - λ_{5}) S_{d}^{*}}{B_{1}}

\frac{\partial H}{\partial u_{2}} |_{u_{2} = u_{2}^{*}} = B_{2} u_{2}^{*} - λ_{2} E_{d}^{*} + λ_{5} E_{d}^{*} = 0 i . e ., u_{2}^{*} = \frac{(λ_{2} - λ_{5}) E_{d}^{*}}{B_{2}}

\frac{\partial H}{\partial u_{3}} |_{u_{3} = u_{3}^{*}} = B_{3} u_{3}^{*} - λ_{6} S_{h}^{*} + λ_{9} S_{h}^{*} = 0 i . e ., u_{3}^{*} = \frac{{(λ}_{6} - λ_{9} {) S}_{h}^{*}}{B_{3}}

\frac{\partial H}{\partial u_{4}} |_{u_{4} = u_{4}^{*}} = B_{4} u_{4}^{*} - λ_{7} E_{h}^{*} + λ_{9} E_{h}^{*} = 0 i . e ., u_{4}^{*} = \frac{(λ_{7} - λ_{9} {) E}_{h}^{*}}{B_{4}}

By standard control arguments involving bounds on the control, then

u_{i}^{*} = \{\begin{matrix} 0, i f φ_{i}^{*} \leq 0 \\ φ_{i}^{*} i f 0 \leq φ_{i}^{*} \leq 1 \\ 1 i f φ_{i}^{*} \geq 1 \end{matrix}

for

i = 1, 2, 3, 4

and where

φ_{1}^{*} = \frac{(λ_{1} - λ_{5}) S_{d}^{*}}{B_{1}}

φ_{2}^{*} = \frac{(λ_{2} - λ_{5}) E_{d}^{*}}{B_{2}}

φ_{3}^{*} = \frac{{(λ}_{6} - λ_{9} {) S}_{h}^{*}}{B_{3}}

φ_{4}^{*} = \frac{(λ_{7} - λ_{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_{d P} = 150,

I_{d F} = 250

R_{d} = 50000,

S_{h} = 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} \neq 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.

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} \neq u_{3} \neq 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.

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} \neq 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} \neq 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} (t), u_{2} (t), u_{3} (t),

and

u_{4} (t)

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.

I C E R = \frac{C h a n g e i n t o t a l c o s t}{c h a n g e i n c o n t r o l b e n e f i t s}

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}, a n d {{\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} \neq u_{2} \neq 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} \neq u_{3} \neq 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} \neq u_{4} \neq 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} \neq u_{3} \neq u_{4} \neq 0)

Strategy 5: Implementing all controls (i.e.,

u_{1} \neq u_{2} \neq u_{3} \neq u_{4} \neq 0)

I C E R (2) = \frac{44,140,000}{1,046,800} = 42.16660, I C E R (3) = \frac{- 5876000}{832380} = - 7.05928

I C E R (2) > I C E R (3)

, 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.

I C E R (3) = \frac{38,264,000}{1,879,180} = 20.36207, I C E R (4) = \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.

I C E R (4) = \frac{37,035,000}{2,053,240} = 18.03734, I C E R (1) = \frac{- 28385900}{295140} = - 96.17775

I C E R (4) > I C E R (1),

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.

I C E R (1) = \frac{8,649,100}{2,348,380} = 3.68301

I C E R (5) = \frac{- 1,977,900}{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.,

β_{d P},

β_{d F},

β_{h P}

, and

β_{h F}

. 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}, β_{d P}, β_{d F}, δ_{d}, α_{d},

and

α_{d}

have a positive relationship with the effective reproduction number

(R_{e})

based on a sensitivity analysis of the parameter. On the other hand, the parameters

μ_{d}, v_{d}, θ_{d} {, c}_{d}, ρ_{d}

, and

ε_{d}

are negative signs, then these are inversely proportional to the effective number

(R_{e})

. 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.

Data Availability Statement

The data used to support the findings of this manuscript are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References

Asfaw, G. B., Abagero, A., Addissie, A., Yalew, A. W., Watere, S. H., Desta, G. B., ... & Deressa, S. G. (2024). Epidemiology of suspected rabies cases in Ethiopia: 2018–2022. One Health Advances, 2(1), 3. [CrossRef]
Shigui Ruan, (2017). Modeling the transmission dynamics and control of rabies in China. Mathematical. Biosciences Volume 286, April, Pages 65-93.
Musaili, J. S., & Chepkwony, I. (2020). A Mathematical Model of Rabies Transmission Dynamics in Dogs Incorporating Public Health Education as a Control Strategy-A Case Study of Makueni County. Journal of Advances in Mathematics and Computer Science, 35(1), 1-11. [CrossRef]
Eze, O. C., Mbah, G. E., Nnaji, D. U., & Onyiaji, N. E. (2020). Mathematical modelling of transmission dynamics of rabies virus. International Journal of Mathematics Trends and Technology-IJMTT, 66.
Thongtha, A., & Modnak, C. (2021). A mathematical modeling of rabies with vaccination and culling. International Journal of Biomathematics, 14(06), 2150039. [CrossRef]
Pantha, B., Giri, S., Joshi, H. R., & Vaidya, N. K. (2021). Modeling transmission dynamics of rabies in Nepal. Infectious Disease Modelling, 6, 284-301. [CrossRef]
Iggidr, A., Mbang, J., Sallet, G., & Tewa, J. J. (2007, August). Multi-compartment models. In Conference Publications (Vol. 2007, No. Special, pp. 506-519). Conference Publications.
Edward, S., & Nyerere, N. (2015). A mathematical model for the dynamics of cholera with control measures. Applied and computational Mathematics, 4(2), 53-63. [CrossRef]
Josephine, E. A. (2009). The basic reproduction number: Bifurcation and stability. A PGD thesis submitted to the African Institute of Mathematical Sciences.
Hailemichael, D. D., Edessa, G. K., & Koya, P. R. (2023). Mathematical Modeling of Dog Rabies Transmission Dynamics Using Optimal Control Analysis. Contemporary Mathematics, 296-319.
Asamoah, J. K. K., Oduro, F. T., Bonyah, E., & Seidu, B. (2017). Modelling of rabies transmission dynamics using optimal control analysis. Journal of Applied Mathematics, 2017.
Ega, T. T., Luboobi, L., & Kuznetsov, D. (2015). Modeling the dynamics of rabies transmission with vaccination and stability analysis.
Hailemichael, D. D., Edessa, G. K., & Koya, P. R. (2022). Effect of vaccination and culling on the dynamics of rabies transmission from stray dogs to domestic dogs. Journal of Applied Mathematics, 2022. [CrossRef]
Pontryagin, L.S. (2018). Mathematical theory of optimal processes. Routledge.
Olaniyi, S. , Okosun, K. O., Adesanya, S. O., & Lebelo, R. S. (2020). Modelling malaria dynamics with partial immunity and protected travellers: optimal control and cost-effectiveness analysis. Journal of biological dynamics, 14(1), 90-115. [CrossRef]
Olaniyi, S.J.A.M. (2018). Dynamics of Zika virus model with nonlinear incidence and optimal control strategies. Appl. Math. Inf. Sci, 12(5), 969-982. [CrossRef]
Modnak, C. (2017). Mathematical modelling of an avian influenza: optimal control study for intervention strategies. Appl. Math. Inf. Sci, 11(4), 1049-1057. [CrossRef]
Hugo, A., Makinde, O. D., Kumar, S., & Chibwana, F. F. (2017). Optimal control and cost effectiveness analysis for Newcastle disease eco-epidemiological model in Tanzania. Journal of Biological Dynamics, 11(1), 190-209. [CrossRef]
Lenhart, S., & Workman, J. T. (2007). Optimal control applied to biological models. Chapman and Hall/CRC.
Renalda, E. K., Kuznetsov, D., & Kreppel, K. (2020). Desirable Dog-Rabies Control Methods in an Urban setting in Africa-a Mathematical Model.
Sharma, P. M., Dobhal, Y., Sharma, V., Kumar, S., & Gupta, S. (2015). Rabies: Understanding Neuropathology and Management. Int. J. Heal. Sci. Res, 5, 320-332.

Figure 1. Model Diagram.

Figure 2. Effect of parameter

β_{d F}

R_{e}

Figure 2. Effect of parameter

β_{d F}

R_{e}

Figure 3. Effect of parameter

β_{d P}

R_{e}

Figure 3. Effect of parameter

β_{d P}

R_{e}

Figure 8. Simulations of the Optimal Control showing the effect of control Strategy 3.

Figure 9. Simulations of the Optimal Control showing the effect of control Strategy 4.

Figure 10. Simulations of the Optimal Control showing the effect of control Strategy 5.

Table 2. Parameters and their values for model (1) (unit: year^-1).

Parameter	Value	Interpretation	Source
$A_{d}$	15900	Dogs annual crop of newborn puppies	Assumed
$β_{d P}$	$1.28 \times 10^{- 8}$	Prodromal dog-to-dog transmission rate	Assumed
$β_{d F}$	$1.10 \times 10^{- 5}$	Furious dog-to-dog transmission rate	Assumed
$δ_{d}$	0.17	The incubation period of dog populations	[12]
$ρ_{d}$	0.821	Rate of prodromal to furious stage	Assumed
$v_{d}$	0.25	Pre-exposure prophylaxis for dogs	[11]
$θ_{d}$	0.2	Post-exposure prophylaxis for dogs	[11]
$ε_{d}$	1	The death rate of dogs due to rabies	[11]
$μ_{d}$	0.056	The natural death rate of dogs	[11]
$c_{d}$	0.01792	Rabid dog culling rate	[20]
$α_{d}$	1	Loss of immunity in dogs	[11]
$A_{h}$	112980	Human annual birth	[12]
$β_{h P}$	$9.79 \times 10^{- 10}$	Prodromal dog-to-human transmission rate	Assumed
$β_{h F}$	$9.86 \times 10^{- 8}$	Furious dog-to-human transmission rate	Assumed
$δ_{h}$	0.1667	The incubation period of human populations	[12]
$v_{h}$	0.2	Pre-exposure prophylaxis for humans	Assumed
$θ_{h}$	0.2	Post-exposure prophylaxis for humans	[11]
$ε_{h}$	1	The death rate of humans due to rabies	[11]
$μ_{h}$	0.0074	Natural death rate of humans	[11]

Table 3. Total number of infections saved and cost averted for all strategies.

Strategy	Total infection Averted	Total cost ($)
1	2,348,380	8,649,100
2	1,046,800	44,140,000
3	1,879,180	38,264,000
4	2,053,240	37,035,000
5	2,628,870	6,671,200

Table 6. Increasing order of infections averted.

Intervention	Total infection	Total cost
Strategy	Averted	($)	ICER
2	1,046,800	44,140,000	42.16660
3	1,879,180	38,264,000	-7.05928
4	2,053,240	37,035,000	-
1	2,348,380	8,649,100	-
5	2,628,870	6,671,200	-

Table 7. Comparison between intervention strategies 3 and 4.

Intervention	Total infection	Total cost
Strategy	Averted	($)	ICER
3	1,879,180	38,264,000	20.36207
4	2,053,240	37,035,000	-7.06078
1	2,348,380	8,649,100	-
5	2,628,870	6,671,200	-

Table 8. Comparison between interventions strategy 4 and 1.

Intervention	Total infection	Total cost
Strategy	Averted	($)	ICER
4	2,053,240	37,035,000	18.03734
1	2,348,380	8,649,100	-96.17775
5	2,628,870	6,671,200	-

Table 9. Comparison between interventions strategy 1 and 5.

Intervention	Total infection	Total cost
Strategy	Averted	($)	ICER
1	2,348,380	8,649,100	3.68301
5	2,628,870	6,671,200	-7.051588

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Cost Effectiveness Analysis of Optimal Rabies Control Strategies

Abstract

1. Introduction

2. Description and Model Formulation

3. Model Analysis

3.1. Positivity

3.2. Boundedness of the Model

3.4. Disease-Free Equilibrium E 0

3.5. Basic Reproduction Number ( R 0 )

3.5.1. Effective Reproduction Number

3.6. Global Stability of Disease-Free Equilibrium Points

4. Sensitivity Analysis

4.1. Interpretation of Sensitivity Indices

5. Analysis of the Optimal Control Model

5.1. Characterization of the Optimal Control

6. Simulations and Cost-Effectiveness Analysis

6.1. Numerical Simulations

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

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

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

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

6.1.5. Intervention 4: Implementing All Optimal Control Strategies

7. Cost-Effectiveness Analysis

8. Conclusion and Recommendations

8.1. Conclusion

8.2. Recommendations

Data Availability Statement

Conflicts of Interest

References

MDPI Initiatives

Important Links

Subscribe

3.4. Disease-Free Equilibrium $E_{0}$

3.5. Basic Reproduction Number $(R_{0})$