Modeling Anthrax with Optimal Control and Cost Effectiveness Analysis

Shaibu Osman^1*, Dominic Otoo², Oluwole Daniel Makinde³

¹Department of Basic Sciences, School of Basic and Biomedical Sciences, University of Health and Allied Sciences, Ho, Volta Region, Ghana.
²Department of Mathematics and Statistics, University of Energy and Natural Resources, Sunyani, Ghana.
³Faculty of Military Science, Stellenbosch University, Saldanha, South Africa.

DOI: 10.4236/am.2020.113020                

Abstract

Anthrax is an infection caused by bacteria and it affects both human and animal populations. The disease can be categorized under zoonotic diseases and humans can contract infections through contact with infected animals, ingest contaminated dairy and animal products. In this paper, we developed a mathematical model for anthrax transmission dynamics in both human and animal populations with optimal control. The qualitative solution of the model behaviour was analyzed by determining R_hv, equilibrium points and sensitivity analysis. A vaccination class was incorporated into the model with waning immunity. Local and global stability of the model’s equilibria was found to be locally asymptotically stable whenever R_hv < 1 and unstable otherwise. Analysis of parameter contribution was conducted to determine the contribution of each to R_hv. It was revealed that reducing animal and human interaction rate, would decrease R_hv. We extended the model to optimal control in order to find the best control strategy in reducing anthrax infections. It showed that the effective strategy in combating the anthrax epidemics is vaccination of animals and prevention of humans.

Keywords

Anthrax Model, Stability Analysis, Sensitivity, Optimal Control, Cost Effectiveness

1. Introduction

Anthrax is an infectious disease that is caused by the bacteria Bacillus anthracis. The disease affects wild, domestic animals and humans. Susceptible individuals contract the anthrax disease if they interact with infected animals or consumed contaminated dairy and animal products. Anthrax is a zoonotic disease found naturally in the soil and it affects both animal and human populations worldwide. Individuals can contract the disease from direct contact with animals or from contaminated animal products [1] [2].

Mathematical models are capable of describing a natural phenomenon and the transmission dynamics of infectious diseases. These models play a key role in combating infectious diseases in epidemiology. These models can explicitly explain the transmission mechanism and dynamics of disease [3] [4]. Complex models for transmission dynamics of diseases such as periodic orbits, Hoff bifurcations and multiple equilibrium have been proposed and worked on for some time now. They give a concise qualitative illustration of the disease dynamics and better analysis and implications for disease prediction of [5] [6].

[1] investigated the effects of constant vaccination on anthrax model but never considered effects of optimal control. From the theoretical results of their study under constant vaccination, the transmission of the disease model is similar to dynamics without vaccination. Moreover, [7] developed a model by considering vaccine impact and concluded by establishing the optimal vaccine coverage threshold required for disease eradication. However, [8] used optimal control in the study of a nonlinear SIR epidemic model with a vaccination strategy.

[9] considered the application of optimal control to investigate the impact of chemo-therapy on malaria disease with infection immigrants and [10] applied optimal control methods associated with preventing exogenous reinfection based on a exogenous reinfection tuberculosis model. Authors in [11], investigated the essential role of three basic controls: personal protection, treatment, and mosquito reduction strategies in combating an infectious disease.

However, [12] formulated a general epidemic model of a vector-borne disease

Figure 1. Cases of anthrax outbreaks in Ghana between 2005-2016.

consisting two vertebrate host species and one insect vector species. The qualitative analysis of the study revealed that model exhibited a multiple endemic equilibrium. The spread of anthrax has an impact on the life’s of people and cost of treatment is of concern in every human endeavour. Many studies in literature have been performed to establish the role of vaccination and treatment on transmission of diseases [13].

Some studies reveal the complex nature of the anthrax disease. The transmission dynamics of the disease varies depending on the prevailing conditions of the country. There has not been any mathematical model developed to explain the transmission dynamics of anthrax with optimal control and cost effectiveness in Ghana. We therefore developed a model to investigate the dynamics of anthrax with optimal control and cost effectiveness analysis using the anthrax data between 2005 to 2016 as shown in Figure 1.

2. Model Formulation and Description

The model divides human and animal populations at any time, (t) into compartments with respect to their disease status as shown in Figure 2. Total animal

population, $N_v\left(t\right)$, is subdivided into Susceptible animals, $\left(S_v\right)$, Infectious animals, $\left(I_v\right)$, Vaccinated animals, $\left(V_v\right)$, and Recovered animals, $\left(R_v\right)$.

Hence total animal population;

$N_v\left(t\right)\mathrm{<mo>\; =\; </mo>}{S}_v\mathrm{<mo\; stretchy="false">\; (\; </mo>}t\mathrm{<mo\; stretchy="false">\; )\; </mo>}+V_v\mathrm{<mo\; stretchy="false">\; (\; </mo>}t\mathrm{<mo\; stretchy="false">\; )\; </mo>}+I_v\mathrm{<mo\; stretchy="false">\; (\; </mo>}t\mathrm{<mo\; stretchy="false">\; )\; </mo>}+R_v\mathrm{<mo\; stretchy="false">\; (\; </mo>}t\mathrm{<mo\; stretchy="false">\; )\; </mo>.}$ (1)

Total human population denoted by $N_h$, is subdivided into Susceptible humans, $\left(S_h\right)$, Infectious humans, $\left(I_h\right)$, and Recovered, $\left(R_h\right)$.

Hence, total human population:

$N_h\left(t\right)=S_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right).$ (2)

Susceptible humans to are recruited at a rate ${\Lambda }_h$. Humans acquire anthrax through direct contact, inhalation and ingestion of contaminated animal products at a rate $\left(I_h+I_v\right)\beta$. Infected humans recover from anthrax disease at a rate $\gamma$. Humans infected with anthrax die at a rate ${\delta }_h$ and the recovered humans may loose immunity and return to the susceptible class at a rate ${\sigma }_h$. The entire human population has a natural death rate of ${\mu }_h$.

Susceptible animals are recruited at a rate ${\Lambda }_v$, only a fraction of the animals are successfully vaccinated at a rate $u_3$, where $u_3\in \left[\mathrm{0,1}\right]$. The disease can be obtained through contacts with infectious animals and humans at a rate $\left(I_h+I_v\right)\lambda$. Natural death rate of all animal compartments is ${\mu }_v$ and death rate as a result of anthrax infections is ${\delta }_v$. Animal recovery rate is $\alpha$ and vaccinated animals may move to the infected class at a rate $b{\beta }_m^{\mathrm{<mo>\; *\; </mo>}}\lambda$ due to waning effect. Where $\left(1-b\right)\in \left[\mathrm{0,1}\right]$ is the efficacy of the vaccine. Rate at which animals may loose immunity and join susceptible class is $\tau$. Where ${\beta }_m^{\mathrm{<mo>\; *\; </mo>}}=I_h+I_v$

The system of ordinary differential equations obtained from the model in Figure 2 are as follows:

Figure 2. Flow diagram for the anthrax disease transmission dynamics.

$\begin{array}{l}\frac{\text{d}{S}_h}{\text{d}t}={\Lambda }_h+{\sigma }_h{R}_h-{\beta }_m^{\mathrm{<mo>\; *\; </mo>}}\beta S_h-{\mu }_h{S}_h\hfill \\ \frac{\text{d}{I}_h}{\text{d}t}={\beta }_m^{\mathrm{<mo>\; *\; </mo>}}\beta S_h-\gamma I_h-\left({\mu }_h+{\delta }_h\right)I_h\hfill \\ \frac{\text{d}{R}_h}{\text{d}t}=\gamma I_h-\left({\sigma }_h+{\mu }_h\right)R_h\hfill \\ \frac{\text{d}{S}_v}{\text{d}t}=\left(1-u_3\right){\Lambda }_v-{\beta }_m^{\mathrm{<mo>\; *\; </mo>}}\lambda S_v-{\mu }_v{S}_v+{\sigma }_v{R}_v+\tau V_v\hfill \\ \frac{\text{d}{I}_v}{\text{d}t}={\beta }_m^{\mathrm{<mo>\; *\; </mo>}}\lambda S_v+b{\beta }_m^{\mathrm{<mo>\; *\; </mo>}}\lambda V_v-\alpha I_v-\left({\mu }_v+{\delta }_v\right)I_v\hfill \\ \frac{\text{d}{R}_v}{\text{d}t}=\alpha I_v-\left({\sigma }_v+{\mu }_v\right)R_v\hfill \\ \frac{\text{d}{V}_v}{\text{d}t}=u_3{\Lambda }_v-\left(\tau +{\mu }_v\right)V_v-b\lambda {\beta }_m^{\mathrm{<mo>\; *\; </mo>}}{V}_v\hfill \end{array}\}$ (3)

3. Model Analysis

3.1. Positivity and Boundedness of Solutions

In epidemiological models, conditions under which the system should have non-negative solutions is paramount. The model would be biologically meaningful if all the solutions with positive initial conditions remain positive at every point in time.

Theorem 1. Let $\begin{array}{l}\Pi =\{\left(S_h\left(t\right)\mathrm{,}{I}_h\left(t\right)\mathrm{,}{R}_h\left(t\right)\mathrm{,}{S}_v\left(t\right)\mathrm{,}{I}_v\left(t\right)\mathrm{,}{V}_v\left(t\right)\mathrm{,}{R}_v\left(t\right)\right)\in {ℝ}_{+}^7\mathrm{<mo>\; :\; </mo>}\\ \left(S_h\left(0\right)\mathrm{,}{I}_h\left(0\right)\mathrm{,}{R}_h\left(0\right)\mathrm{,}{S}_v\left(0\right)\mathrm{,}{I}_v\left(0\right)\mathrm{,}{V}_v\left(0\right)\mathrm{,}{R}_v\left(0\right)\right)>0\}\end{array}$, then the solution of $\left\{\left(S_h\left(t\right)\mathrm{,}{I}_h\left(t\right)\mathrm{,}{R}_h\left(t\right)\mathrm{,}{S}_v\left(t\right)\mathrm{,}{I}_v\left(t\right)\mathrm{,}{V}_v\left(t\right)\mathrm{,}{R}_v\left(t\right)\right)\right\}$ are non-negative for all time $t\ge 0$.

Hence, if $S_h\left(0\right)\mathrm{,}{I}_h\left(0\right)\mathrm{,}{R}_h\left(0\right)\mathrm{,}{S}_v\left(0\right)\mathrm{,}{I}_v\left(0\right)\mathrm{,}{V}_v\left(0\right)\mathrm{,}{R}_v\left(0\right)$, then $S_h\left(t\right)\mathrm{,}{I}_h\left(t\right)\mathrm{,}{R}_h\left(t\right)\mathrm{,}{S}_v\left(t\right)\mathrm{,}{I}_v\left(t\right)\mathrm{,}{V}_v\left(t\right)\mathrm{,}{R}_v\left(t\right)$ are also non-negative for all $t>0$.

Human total population at any given time:

$N_h\left(t\right)=S_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right).$ (4)

$\frac{\text{d}{N}_h}{\text{d}t}=\frac{\text{d}{S}_h}{\text{d}t}+\frac{\text{d}{I}_h}{\text{d}t}+\frac{\text{d}{R}_h}{\text{d}t}.$ (5)

$\frac{\text{d}{N}_h}{\text{d}t}={\Lambda }_h-{\mu }_h{S}_h-\left({\mu }_h+{\delta }_h\right)I_h-\left({\sigma }_h+{\mu }_h\right)R_h\mathrm{.}$

In the absence of mortality due to Anthrax infections, the above equation becomes;

$\frac{\text{d}{N}_h}{\text{d}t}\le {\Lambda }_h-{\mu }_h{N}_h\mathrm{.}$ (6)

Solving the equation and as $t\to \infty$, the population size, $N_h\to \frac{{\Lambda }_h}{{\mu }_h}$.

$0\le N_h\le \frac{{\Lambda }_h}{{\mu }_h}$ and $N_h\left(t\right)\le \frac{{\Lambda }_h}{{\mu }_h}$.

Also, if

$N_h\left(0\right)\le \frac{{\Lambda }_h}{{\mu }_h}$, then $N_h\left(t\right)\le \frac{{\Lambda }_h}{{\mu }_h}$.

${\Pi }_h=\left\{\left(S_h\mathrm{,}{I}_h\mathrm{,}{R}_h\right)\in {ℝ}_3^3\mathrm{<mo>\; :\; </mo>}{S}_h+I_h+R_h\le \frac{{\Lambda }_h}{{\mu }_h}\right\}$ (7)

However, for total animal (vector) population at any given time:

$N_v\left(t\right)=S_v\left(t\right)+V_v\left(t\right)+I_v\left(t\right)+R_v\left(t\right).$ (8)

$\frac{\text{d}{N}_v}{\text{d}t}=\frac{\text{d}{S}_v}{\text{d}t}+\frac{\text{d}{V}_v}{\text{d}t}+\frac{\text{d}{I}_v}{\text{d}t}+\frac{\text{d}{R}_v}{\text{d}t}.$ (9)

$\frac{\text{d}{N}_v}{\text{d}t}={\Lambda }_v-{\mu }_v{N}_v-{\delta }_v{I}_v\mathrm{.}$

$\frac{\text{d}{N}_v}{\text{d}t}\le {\Lambda }_v-{\mu }_v{N}_v\mathrm{.}$ (10)

Solving the differential equation and as $t\to \infty$, the population size, $N_v\to \frac{{\Lambda }_v}{{\mu }_v}$.

$0\le N_v\le \frac{{\Lambda }_v}{{\mu }_v}$ and $N_v\left(t\right)\le \frac{{\Lambda }_v}{{\mu }_v}$.

Also, if

$N_v\left(0\right)\le \frac{{\Lambda }_v}{{\mu }_v}$, then $N_v\left(t\right)\le \frac{{\Lambda }_v}{{\mu }_v}$.

${\Pi }_v=\left\{\left(S_v,I_v,R_v,V_v\right)\in {ℝ}_{+}^4:S_v+I_v+R_v+V_v\le \frac{{\Lambda }_v}{{\mu }_v}\right\}.$ (11)

Hence, feasible region for the dynamical system in (1) is given by:

$\Pi ={\Pi }_h\times {\Pi }_v\subset {ℝ}_{+}^3\times {ℝ}_4^4.$ (12)

${\Pi }_h=\left\{\left(S_h\mathrm{,}{I}_h\mathrm{,}{R}_h\right)\in {ℝ}_3^3\mathrm{<mo>\; :\; </mo>}{S}_h+I_h+R_h\le \frac{{\Lambda }_h}{{\mu }_h}\right\}$ (13)

${\Pi }_v=\left\{\left(S_v,I_v,R_v,V_v\right)\in {ℝ}_{+}^4:S_v+I_v+R_v+V_v\le \frac{{\Lambda }_v}{{\mu }_v}\right\}.$ (14)

where $\Pi$ is positively invariant.

3.2. Disease Free Equilibrium

The Disease Free Equilibrium (DFE) of the dynamical system in (3) only exists in the absence of infections. We can then compute DFE by setting the system of equations in (3) to zero and establish DFE points.

At DFE, there are no infections and recovery. Hence, $R_v^{*}=I_h^{*}=I_v^{*}=R_h^{*}=0$.

${\xi }_0=\left(S_h^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}_h^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{R}_h^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{S}_v^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}_v^{\mathrm{<mo>\; *\; </mo>}},R_v^{\mathrm{<mo>\; *\; </mo>}},V_v^{\mathrm{<mo>\; *\; </mo>}}\right)\mathrm{.}$

$\begin{array}{l}\frac{\text{d}{S}_h}{\text{d}t}={\Lambda }_h+{\sigma }_h{R}_h-{\beta }_m^{*}\beta S_h-{\mu }_h{S}_h=0\hfill \\ \frac{\text{d}{V}_v}{\text{d}t}=u_3{\Lambda }_v-\left(\tau +{\mu }_v\right)V_v-b\lambda {\beta }_m^{*}{V}_v=0\hfill \\ \frac{\text{d}{S}_v}{\text{d}t}=\left(1-u_3\right){\Lambda }_v-{\beta }_m^{*}\lambda S_v-{\mu }_v{S}_v+{\sigma }_v{R}_v+\tau V_v=0\hfill \end{array}\}$ (15)

The DFE point is given by;

${\xi }_0=\left(\frac{{\Lambda }_h}{{\mu }_h}\mathrm{,0,0,}\frac{{\Lambda }_v\left(\tau +{\mu }_v\left(1-u_3\right)\right)}{{\mu }_v\left(\tau +{\mu }_v\right)}\mathrm{,0,0,}\frac{u_3{\Lambda }_v}{\tau +{\mu }_v}\right)\mathrm{.}$ (16)

3.3. The Basic Reproductive Number

Using the next generation matrix approach in [4], the reproductive rate can be established. The reproductive rate combines the biology of infections with the social and behaviour of the factors influencing contact rate.

The basic reproductive rate refers to the number of secondary cases one infectious individual will produce in a completely susceptible population [14] [15]. This is a threshold parameter that governs the spread of a disease.

Definition 2. The spectral radius of a matrix A is defined as the maximum of the absolute values of the eigenvalues of the matrix A: $\rho \left(A\right)=\sup \left\{\left|\lambda \right|:\lambda \in \rho \left(A\right)\right\}$, where $\rho \left(A\right)$ denotes the set of eigenvalues of the matrix A.

Considering only the infective compartments in the system of differential equations in (3):

$\begin{array}{l}\frac{\text{d}{I}_h}{\text{d}t}=\left(I_h+I_v\right)\beta S_h-\gamma I_h-\left({\mu }_h+{\delta }_h\right)I_h\hfill \\ \frac{\text{d}{I}_v}{\text{d}t}=\left(I_h+I_v\right)\lambda S_v+\left(I_h+I_v\right)b\lambda V_v-\alpha I_v-\left({\mu }_v+{\delta }_v\right)I_v\hfill \end{array}\}$ (17)

Let f be the number of new infection coming into the system and v be the number of infectives that are leaving the system either by death or birth.

$f=\left[\begin{array}{c}\left(I_h+I_v\right)\beta S_h\\ \left(I_h+I_v\right)\lambda S_v+\left(I_h+I_v\right)b\lambda V_v\end{array}\right]$,

$v=\left[\begin{array}{c}\gamma I_h+\left({\mu }_h+{\delta }_h\right)I_h\\ \alpha I_v+\left({\mu }_v+{\delta }_v\right)I_v\end{array}\right]$.

The Jacobian matrix of f and v at DFE is obtained by F and V as follows:

$F=\left[\begin{array}{cc}\beta S_h^{\mathrm{<mo>\; *\; </mo>}}& \beta S_h^{\mathrm{<mo>\; *\; </mo>}}\\ \lambda S_v^{\mathrm{<mo>\; *\; </mo>}}+b\lambda V_v^{\mathrm{<mo>\; *\; </mo>}}& \lambda S_v^{\mathrm{<mo>\; *\; </mo>}}+b\lambda V_v^{\mathrm{<mo>\; *\; </mo>}}\end{array}\right]$,

$V=\left[\begin{array}{cc}\gamma +\left({\mu }_h+{\delta }_h\right)& 0\\ 0& \alpha +\left({\mu }_v+{\delta }_v\right)\end{array}\right]$.

Computing the product of $F{V}^{-1}$.

$F{V}^{-1}=\left[\begin{array}{cc}\frac{\beta S_h^{\mathrm{<mo>\; *\; </mo>}}}{\gamma +\left({\mu }_h+{\delta }_h\right)}& \frac{\beta S_h^{\mathrm{<mo>\; *\; </mo>}}}{\alpha +\left({\mu }_v+{\delta }_v\right)}\\ \frac{\lambda S_v^{\mathrm{<mo>\; *\; </mo>}}+b\lambda V_v^{\mathrm{<mo>\; *\; </mo>}}}{\gamma +\left({\mu }_h+{\delta }_h\right)}& \frac{\lambda S_v^{\mathrm{<mo>\; *\; </mo>}}+b\lambda V_v^{\mathrm{<mo>\; *\; </mo>}}}{\alpha +\left({\mu }_v+{\delta }_v\right)}\end{array}\right]$

By computing the eigenvalues of $F{V}^{-1}$ and selecting the dominant eigenvalue. The eigenvalues are as follows;

0 and $\left[\left(\frac{\lambda S_v^{\mathrm{<mo>\; *\; </mo>}}+b\lambda V_v^{\mathrm{<mo>\; *\; </mo>}}}{\alpha +\left({\mu }_v+{\delta }_v\right)}\right)+\left(\frac{\beta S_h^{\mathrm{<mo>\; *\; </mo>}}}{\gamma +\left({\mu }_h+{\delta }_h\right)}\right)\right]\mathrm{.}$

The dominant eigenvalue is $\left(\frac{\lambda S_v^{\mathrm{<mo>\; *\; </mo>}}+b\lambda V_v^{\mathrm{<mo>\; *\; </mo>}}}{\alpha +\left({\mu }_v+{\delta }_v\right)}\right)+\left(\frac{\beta S_h^{\mathrm{<mo>\; *\; </mo>}}}{\gamma +\left({\mu }_h+{\delta }_h\right)}\right)$.

This implies that the basic reproductive rate is given by;

$R_{hv}=\left[\left(\frac{\beta S_h^{\mathrm{<mo>\; *\; </mo>}}}{\gamma +\left({\mu }_h+{\delta }_h\right)}\right)+\left(\frac{\lambda S_v^{\mathrm{<mo>\; *\; </mo>}}+b\lambda V_v^{\mathrm{<mo>\; *\; </mo>}}}{\left[\alpha +\left({\mu }_v+{\delta }_v\right)\right]}\right)\right]\mathrm{.}$ (18)

At DFE, it becomes;

$R_{hv}=\left(\frac{\beta {\Lambda }_h}{{\mu }_h\left[\gamma +\left({\mu }_h+{\delta }_h\right)\right]}\right)+\frac{b\lambda {\Lambda }_v\left(\tau +{\mu }_v\left(1-2u_3\right)\right)}{{\mu }_v\left(\tau +{\mu }_v\right)\left[\alpha +\left({\mu }_v+{\delta }_v\right)\right]}$ (19)

where;

$R_h\mathrm{<mo>\; =\; </mo>}\left(\frac{\beta {\Lambda }_h}{{\mu }_h\left[\gamma +\left({\mu }_h+{\delta }_h\right)\right]}\right),R_v=\frac{b\lambda {\Lambda }_v\left(\tau +{\mu }_v\left(1-2u_3\right)\right)}{{\mu }_v\left(\tau +{\mu }_v\right)\left[\alpha +\left({\mu }_v+{\delta }_v\right)\right]}.$

3.4. Global Stability of the Disease-Free Equilibrium

Proposition 3. The disease-free equilibrium (DFE) of model (3) is locally asymptotically stable if $R_{hv}<1$ and unstable if $R_{hv}>1$.

Theorem 4. If $R_{hv}\le 1$, the disease-free equilibrium is globally asymptotically stable in the interior of $\Omega$.

Proof. Considering the Lyapunov function below, □

$P\left(t\right)=\left(\alpha +{\mu }_v+{\delta }_v\right)I_h+\left(\gamma +{\mu }_h+{\delta }_h\right)I_v$ (20)

Taking time derivative of $P\left(t\right)$ along solutions of the differential equations in (3);

$\begin{array}{c}\frac{\text{d}P\left(t\right)}{\text{d}t}=\left(\alpha +{\mu }_v+{\delta }_v\right)\frac{\text{d}{I}_h}{\text{d}t}+\left(\gamma +{\mu }_h+{\delta }_h\right)\frac{\text{d}{I}_v}{\text{d}t}\\ =\left(\beta S_h\left(I_h+I_v\right)-\gamma I_h-\left({\mu }_h+{\delta }_h\right)I_h\right)+\lambda S_v\left(I_h+I_v\right)\\ +b\left(I_h+I_v\right)\lambda V_v-\alpha I_v-\left({\mu }_v+{\delta }_v\right)I_v\\ \le \left(\alpha +{\mu }_v+{\delta }_v\right)\frac{\beta {\Lambda }_h{I}_h}{{\mu }_h}+\left(\alpha +{\mu }_v+{\delta }_v\right)\frac{\beta {\Lambda }_h{I}_v}{{\mu }_h}\\ -\left(\alpha +{\mu }_v+{\delta }_v\right)\left(\gamma +{\mu }_h+{\delta }_h\right)I_h\end{array}$

$\begin{array}{l}+I_h\left(\gamma +{\mu }_h+{\delta }_h\right)\left(\frac{\lambda {\Lambda }_v\left(\tau +{\mu }_v\left(I-u_3\right)\right)}{{\mu }_v\left(\tau +{\mu }_v\right)}\right)\\ +I_v\left(\gamma +{\mu }_h+\delta \right)\left(\frac{\lambda {\Lambda }_v\left(\tau +{\mu }_v\left(I-u_3\right)\right)}{{\mu }_v\left(\tau +{\mu }_v\right)}\right)\\ +I_h\left(\gamma +{\mu }_h+{\delta }_h\right)\left(\frac{b{u}_3\lambda {\Lambda }_v}{\tau +{\mu }_v}\right)+I_v\left(\gamma +{\mu }_h+\delta \right)\left(\frac{b{u}_3\lambda {\Lambda }_v}{\tau +{\mu }_v}\right)\\ -I_v\left(\gamma +{\mu }_h+\delta \right)\left(\alpha +{\mu }_v+{\delta }_v\right)\end{array}$

$\begin{array}{l}\le -I_h\left(\gamma +{\mu }_h+{\delta }_h\right)\left(\alpha +{\mu }_v+{\delta }_v\right)\left(1-R_{hv}\right)\\ -I_v\left(\gamma +{\mu }_h+{\delta }_h\right)\left(\alpha +{\mu }_v+{\delta }_v\right)\left(1-R_{hv}\right)\\ =-\left(I_h+I_v\right)\left(\gamma +{\mu }_h+{\delta }_h\right)\left(\alpha +{\mu }_v+{\delta }_v\right)\left(1-R_{hv}\right)\end{array}$ (21)

Time derivative of P along the solutions of differential equations in (3) yields:

$\left(\frac{\text{d}P\left(t\right)}{\text{d}t}\right)\le 0$, if and only if $R_{hv}<0$

$\left(\frac{\text{d}P\left(t\right)}{\text{d}t}\right)=1$, if and only if $I_h+I_v=0$ or $R_{hv}=1$.

The highest compact invariant set in $S_h\mathrm{,}{I}_h\mathrm{,}{I}_v\in \Omega$, $\frac{\text{d}P\left(t\right)}{\text{d}t}=0$, if $R_{hv}\le \text{1}$, is the singleton, ${\xi }_0$.

It shows that ${\xi }_0$ is globally asymptotically stable in $\Omega$. By LaSalle’s invariant principle [16].

3.5. Endemic Equilibrium

Considering the system in 3, at equilibrium, ${\beta }_m^{*}=I_h+I_v=0$. This corresponds to the DFE or the relation:

$H_0{\beta }_m^{\mathrm{<mo>\; *\; </mo>3}}+H_1{\beta }_m^{\mathrm{<mo>\; *\; </mo>2}}+H_2{\beta }_m^{\mathrm{<mo>\; *\; </mo>}}+H_3=0.$ (22)

Remark. The system of equations in (3) would have an EE, $E^{\mathrm{<mo>\; *\; </mo>}}$, if $R_{hv}>1$. This is satisfied by cases (2, 4, 6) in Table 1. The system would have more than one EE point if $R_{hv}>1$. This is satisfied by case (8) as shown in Table 1. The system in Equation (2.3), would have more than one EE point if $R_{hv}<1$, as satisfied by case (3, 5, 7).

From Table 1, multiple EE exists when $\left(R_{hv}\right)$ is less than unity. This indicates the tendency of backward bifurcation. This is a situation where the DFE and EE coexists. The existence of backward bifurcation has significant implication

Table 1. Possible positive real roots of $P\left({\beta }_m^{\mathrm{<mo>\; *\; </mo>}}\right)$ for $R_{hv}>1$ and $R_{hv}<1$.

for epidemiological control measures because an epidemic may persist at steady state even when $\left(R_{hv}<1\right)$.

3.6. Global Stability of Endemic Equilibrium

In this section, the global behaviour of the system in Equation (3) is analysed. Considering the non-linear Lyapunov function.

Theorem 5. The system in Equation (3), is said to have a unique endemic equilibrium if $R_{hv}>1$, and it is globally asymptotically stable.

The EE can only exists if and only if $R_{hv}>1$. So by letting $R_{hv}>1$, it implies that the EE exists.

Considering the non-linear Lyapunov function bellow;

$\begin{array}{c}L=S_h^{**}\left(\frac{S_h}{S_h^{**}}-\ln \frac{S_h}{S_h^{**}}\right)+I_h^{**}\left(\frac{I_h}{I_h^{**}}-\ln \frac{I_h}{I_h^{**}}\right)+\frac{g_1{R}_h^{**}}{\gamma }\left(\frac{R_h}{R_h^{**}}-\ln \frac{R_h}{R_h^{**}}\right)\\ +S_v^{**}\left(\frac{S_v}{S_v^{**}}-\ln \frac{S_v}{S_v^{**}}\right)+I_v^{**}\left(\frac{I_v}{I^{**}}-\ln \frac{I_v}{I_v^{**}}\right)\\ +R_v^{**}\left(\frac{R_v}{R_h^{**}}-\ln \frac{R_v}{R^{**}}\right)+V_v^{**}\left(\frac{Vv}{V_v^{**}}-\ln \frac{V_v}{V_v^{**}}\right).\end{array}$ (23)

where;

$g_1=\left(\gamma +{\mu }_h+{\delta }_h\right)\mathrm{,}{g}_2=\left({\sigma }_2+{\mu }_h\right)\mathrm{,}$

$g_3=\left(\alpha +{\mu }_v+{\delta }_v\right)\mathrm{,}{g}_4=\left({\sigma }_v+{\mu }_v\right)\mathrm{,}$

when the above Lyapunov function is differentiated with respect to time, we obtain the equation;

$\begin{array}{c}\frac{\text{d}L}{\text{d}t}=\left(1-\frac{S_h^{**}}{S_h}\right)\frac{\text{d}{S}_h}{\text{d}t}+\left(1-\frac{I_h^{**}}{I_h}\right)\frac{\text{d}{I}_h}{\text{d}t}+\frac{g_1}{\gamma }\left(1-\frac{R^{**}}{R_h}\right)\frac{\text{d}{R}_h}{\text{d}t}\\ +\left(1-\frac{S_v^{**}}{S_v}\right)\frac{\text{d}{S}_v}{\text{d}t}+\left(1-\frac{I_v^{**}}{I_v}\right)\frac{\text{d}{I}_v}{\text{d}t}\\ +\left(1-\frac{R_v^{**}}{R_v}\right)\frac{\text{d}{R}_v}{\text{d}t}+\left(1-\frac{V_v^{**}}{V_v}\right)\frac{\text{d}{V}_v}{\text{d}t}.\end{array}$ (24)

Therefore, this implies that;

$\begin{array}{l}\frac{\text{d}L}{\text{d}t}=\left(1-\frac{S_h^{**}}{S_h}\right)\left[{\Lambda }_h+{\sigma }_h{R}_h^{**}+\beta {\beta }_m^{**}{S}_h^{**}+{\mu }_h{S}_h^{**}-{\Lambda }_h-\sigma R_h-\beta {\beta }_m{S}_h-{\mu }_h{S}_h\right]\\ +\left(1-\frac{I_h^{**}}{I_h}\right)\left[\beta {\beta }_m{S}_h-g_1{I}_h\right]+\frac{g_1}{\gamma }\left(1-\frac{R_h^{**}}{R_h}\right)\left[\gamma I_h-g_2{R}_h\right]\\ +\left(1-\frac{S_v^{**}}{S_v}\right)[\left(1-u_3\right){\Lambda }_v+\lambda \beta {\beta }_m^{**}{S}_v^{**}+{\mu }_v{S}_v^{**}+{\sigma }_v{R}_v^{**}+\tau V_v^{**}\\ -\left(1-u_3\right){\Lambda }_v-\lambda {\beta }_m{S}_v-{\mu }_v{S}_v-{\sigma }_v{R}_v-\tau V_v]\\ +\left(1-\frac{I_v^{**}}{I_v}\right)\left[\lambda {\beta }_m{S}_v+b\lambda {\beta }_m{V}_v-g_3{I}_v\right]+\frac{g_3}{\alpha }\left(1-\frac{R_v^{**}}{R_v}\right)\left[\alpha I_v-g_4{R}_v\right]\\ +\left(1-\frac{V_v^{**}}{V_v}\right)\left[u_3{\Lambda }_v+b\lambda {\beta }_m^{**}{V}_v^{**}+\left(\tau +{\mu }_v\right)V_v^{**}-u_3{\Lambda }_v-b\lambda {\beta }_m{V}_v-\left(\tau +{\mu }_v\right)V_v\right]\end{array}$ (25)

By further simplification;

Arithmetic mean value exceeds the geometric mean value [17] [18]. This follows that;

$\begin{array}{l}2-\frac{S_h^{**}}{S_h}-\frac{S_h}{S_h^{**}}\le 0\\ 1-\frac{R_h}{R_h^{**}}\le 0\\ 1-\frac{R_h^{**}}{R_h}-\frac{g_1{g}_2{S}_h}{\gamma S_h^{**}}\left(1-\frac{R_h^{**}}{R_h}\right)\le 0\\ 1-\frac{{\beta }_m}{{\beta }_h^{**}}-\frac{S_h^{**}}{S_h}-\frac{S_h{\beta }_m{I}_h^{**}}{S_h^{**}{\beta }_m^{**}}\le 0\\ 1-\frac{I_h}{I_h^{**}}-\frac{\gamma I_h}{I_h^{**}}\left(1-\frac{R_h^{**}}{R_h}\right)\le 0\end{array}$

$\begin{array}{l}2-\frac{S_v^{**}}{S_v}-\frac{S_v}{S_v^{**}}\le 0\\ 1-\frac{S_v^{**}}{S_v}-\frac{R_v}{R_v^{**}}+\frac{R_v{S}_v^{**}}{R_v^{**}{S}_v}\le 0\\ 1-\frac{S_v^{**}}{S_v}-\frac{V_v}{V_v^{**}}-\frac{V_v{S}_v^{**}}{V_v^{**}{S}_v}\le 0\\ 1-\frac{S_v^{**}}{S_v}+\frac{{\beta }_m}{{\beta }_h^{**}}-\frac{S_v{\beta }_m{I}_v^{**}}{S_v^{**}{\beta }_m^{**}{I}_v}\le 0\end{array}$

$\begin{array}{l}1-\frac{I_v}{I_v^{**}}-g_3\alpha \frac{I_v}{I_v^{**}}-g_3\alpha \frac{I_v{R}_v^{**}}{I_v^{**}{R}_v}\le 0\\ 1-\frac{R_v}{R_v^{**}}\le 0\\ 2-\frac{V_v^{**}}{V_v}-\frac{V_v}{V_v^{**}}\le 0\end{array}$

$1-\frac{V_v^{**}}{V_v}+\frac{{\beta }_m}{{\beta }_h^{**}}-\frac{V_v{\beta }_m{I}_v^{**}}{V_v^{**}{\beta }_m^{**}{I}_v}$ (26)

Since all model parameters are non-negative, then the derivative of the Lyapunov function is less than zero $\left(\frac{\text{d}L}{\text{d}t}\le 0\right)$, if the $R_{hv}$ of the system in Equation (3) is greater than one. $\left(R_{hv}>1\right)$. By LaSalle’s Invariant Principle [16], as t approaches infinity, all the solution of the system approach the EE point if $R_{hv}>1$.

Backward Bifurcation and Multiple Equilibrium

In this section, we discuss the phenomenon of Backward Bifurcation. The bifurcating EE point exists only if $\left(R_{hv}<1\right)$. The model has exhibited this property and backward bifurcation exists. When backward bifurcation occurs, the range of $R_{hv}$ is between $R_{hv}^{*}<R_{hv}<1$. There exits at least one EE. Usually, at least one is stable and the DFE is not globally stable when $\left(R_{hv}<1\right)$. In this scenario, the infection would exist even when $R_{hv}<1$.

4. Sensitivity Analysis

Basically, the essence of sensitivity analysis is to determine the contribution of each parameter to the reproductive rate. This is help to identify the parameters with high impact on $R_{hv}$. The basic reproductive rate is usually analysed to find out whether or not treatment of the infectives, mortality and vaccination could help in the control or eradication of the disease in the population [19] [20].

Definition 6. The normalised forward sensitivity index of a variable, y, which depends differentially on a parameter, x, defined as:

${\Upsilon }_x^y=\frac{\partial y}{\partial x}\times \frac{x}{y}\mathrm{.}$ (27)

Sensitivity Indices of $R_{hv}$

In epidemiological models, the value of $R_{hv}$ determines the ability of the infection to spread within the population. We will determine the reduction in infection due to the diseases by computing the sensitivity indices of $R_{hv}$, with respect to the parameter values in the model. The sensitivity indices serve as determinants of the contribution of each parameter in the dynamics of the diseases. Considering all the parameters of the system in model (3), we derive the sensitivity of $R_{hv}$ to each of the parameters in the model.

The sensitivity indices of $R_{hv}$ with respect to each of the parameters of the system in model (3), are given in Table 2.

The detailed sensitivity analysis of $\left(R_{hv}\right)$ showed that increasing $\alpha$ would decrease $R_{hv}$. Moreover, decreasing $\alpha$ would increase $\left(R_{hv}\right)$. Also, an increase in the values of ${\Lambda }_h$, ${\Lambda }_v$, $\beta$ and $\lambda$ would cause an increase in $R_{hv}$ and a decrease in the values of ${\Lambda }_h$, ${\Lambda }_v$, $\beta$ and $\lambda$ would cause a decrease in $R_{hv}$.

Table 2. Sensitivity indices of parameters to $R_{hv}$.

5. Optimal Control of the Anthrax Model

The essence optimal control in disease modelling to to determine the best control strategy in fighting the spread of the infection and the costs associated with this strategy. This analysis is carried out to establish the impact of various intervention schemes. Optimal control problem is established by introducing controls into the Anthrax model (3) and the formation of an objective functional which seeks to minimise: $\left(u_1,u_2\mathrm{,}{u}_3\mathrm{,}{u}_4\right)$.

Preventive control measure on susceptible humans, $\left(u_1\right)$ : This is to reduce the acquisition of the disease. Treatment control on infected humans, $\left(u_2\right)$ : This is to minimise infections. Vaccination as a control measure on susceptible animals, $\left(u_3\right)$ : By using antimicrobial drugs. Treatment control on infected animals, $\left(u_4\right)$ : This is to reduce infections.

The objective functional in achieving this purpose is as follows:

$J=\underset{u_1,u_2,u_3,u_4}{\min }{\displaystyle {\int }_0^{t_f}}\left(A_1{I}_v+A_2{I}_h+A_3{u}_1^2+A_4{u}_2^2+A_{5}m{u}_3^2+A_6{u}_4^2\right)\text{d}t.$ (28)

subject to the system in (3). Where; $\left(A_1\mathrm{,}{A}_2\mathrm{,}{A}_3\mathrm{,}{A}_4\mathrm{,}{A}_5\mathrm{,}{A}_5\right)$ are generally referred as balancing cost factors (weight constants) to help balance all terms in the integral thereby avoiding dominance of one over the other. $A_1{I}_h,A_2{I}_v$ are the costs associated with $I_h$ and $I_v$. $A_3{u}_1^2$ is the cost associated with $S_h$. $A_{5}m{u}_3^2$ is the cost associated with vaccination of $S_v$ and m is the number of $I_v$. $A_4{u}_2^2$ is the costs of treatment of $I_h$. $A_6{u}_4^2$ is the cost of treatment of $I_v$. $t_f$ is intervention period. $A_1{I}_h\mathrm{,}{A}_2{I}_v$ represents a linear function for cost of infection and $A_3{u}_1^2\mathrm{,}{A}_4{u}_2^2\mathrm{,}{A}_{5}m{u}_3^2\mathrm{,}{A}_6{u}_4^2$ denotes a quadratic function for cost of controls.

Control efforts of model is by linear combination $u_i^2\left(t\right)$, $\left(i=\mathrm{1,2}\right)$. A quadratic in nature by assumption that cost is non-linear naturally [21] [22] [23]. The purpose is to reduce infection and cost of treatment.

This is finding the optimal functions; $\left(u_1^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}_2^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}_3^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}_4^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\right)$ such that;

$J\left(u_1^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}_2^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}_3^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}_4^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\right)=\underset{u_1,u_2,u_3,u_4\in \cup }{\min }J\left(u_1\mathrm{,}{u}_2\mathrm{,}{u}_3\mathrm{,}{u}_4\right)$ (29)

where;

$\cup =\left\{u:u\mathrm{,0}\le u_i\left(t\right)\le 1,t\in \left[\mathrm{0,}{t}_f\right]\mathrm{,}i=\mathrm{1,2,3,4}\right\}$ and is called the control set.

Pontryagin’s Maximum Principle

This principle gives the necessary conditions that an optimal must satisfy. It changes the system in (3) and (29) into minimisation problem point-wise Hamiltonian, (H) with respect to $\left(u_1,u_2\mathrm{,}{u}_3\mathrm{,}{u}_4\right)$.

$\begin{array}{c}H=A_1{I}_v+A_2{I}_h+A_3{u}_1^2+A_4{u}_2^2+A_{5}m{u}_3^2+A_6{u}_6^2\\ +{\lambda }_1\left\{{\Lambda }_h+{\sigma }_h{R}_h-\left(1-u_1\right)\beta \left(I_v+I_h\right)S_h-{\mu }_h{S}_h\right\}\\ +{\lambda }_2\left\{\left(1-u_1\right){\beta }_m^{*}\beta S_h-(u_2+\gamma )I_h-({\mu }_h+{\delta }_h)I_h\right\}\\ +{\lambda }_3\left\{\left(u_2+\gamma \right)I_h-\left({\sigma }_h+{\mu }_h\right)R_h\right\}\\ +{\lambda }_4\left\{\left(1-u_3\right){\Lambda }_v-\left(1-u_1\right)\lambda \left(I_v+I_h\right)S_v-\left({\mu }_v{S}_v+{\sigma }_v{R}_v+\tau V_v\right)\right\}\end{array}$

$\begin{array}{l}+{\lambda }_5\left\{\left(1-u_1\right){\beta }_m^{*}\lambda S_v+\left(1-u_1\right)b{\beta }_m^{*}\lambda V_v-(u_4+\alpha )I_v-({\mu }_v+{\delta }_v)I_v\right\}\\ +{\lambda }_6\left\{\left(u_4+\alpha \right)I_v-\left({\sigma }_v+{\mu }_v\right)R_v\right\}\\ +{\lambda }_7\left\{u_3{\Lambda }_v-\left(\tau +{\mu }_v\right)V_v-\left(1-u_1\right)b\lambda \left(I_v+I_h\right)V_v\right\}\end{array}$ (30)

where ${\lambda }_1\mathrm{,}{\lambda }_2\mathrm{,}{\lambda }_3\mathrm{,}{\lambda }_4\mathrm{,}{\lambda }_5\mathrm{,}{\lambda }_6$ and ${\lambda }_7$ are the adjoint variables.

The solutions of adjoint are system;

$\begin{array}{l}\frac{\text{d}{\lambda }_1}{\text{d}t}=\left(1-u_1\right)\left(I_v+I_h\right)\beta \left({\lambda }_1-{\lambda }_2\right)+{\mu }_h{\lambda }_1\hfill \\ \frac{\text{d}{\lambda }_2}{\text{d}t}=-A_2+\left(1-u_1\right)\beta S_h\left({\lambda }_1-{\lambda }_2\right)+\left(u_2+\gamma \right)\left({\lambda }_2-{\lambda }_3\right)+{\lambda }_2\left({\mu }_h+{\delta }_h\right)\hfill \\ +\left(1-u_1\right)\lambda S_v\left({\lambda }_4-{\lambda }_5\right)+(1-u_1)b\lambda V_v\left({\lambda }_7-{\lambda }_5\right)\hfill \\ \frac{\text{d}{\lambda }_3}{\text{d}t}=-{\sigma }_h{\lambda }_1+\left({\sigma }_h+{\mu }_h\right){\lambda }_3\hfill \\ \frac{\text{d}{\lambda }_4}{\text{d}t}=\left(1-u_1\right)\lambda \left(I_v+I_h\right)\left({\lambda }_4-{\lambda }_5\right)+{\mu }_v{\lambda }_4\hfill \\ \frac{\text{d}{\lambda }_5}{\text{d}t}=-A_1+\left(1-u_1\right)\beta S_h\left({\lambda }_1-{\lambda }_2\right)+\left(1-u_1\right)\lambda S_v\left({\lambda }_4-{\lambda }_5\right)\hfill \\ +b\lambda \left({\lambda }_7-{\lambda }_5\right)V_v+\left({\lambda }_5-{\lambda }_6\right)\left(u_4+\alpha \right)+{\lambda }_5\left({\mu }_v+{\delta }_v\right)\hfill \\ \frac{\text{d}{\lambda }_6}{\text{d}t}=-{\sigma }_v{\lambda }_4+\left({\sigma }_v+{\mu }_v\right){\lambda }_6\hfill \\ \frac{\text{d}{\lambda }_7}{\text{d}t}=-\tau {\lambda }_4+\left(1-u_1\right)b\lambda \left(I_v+I_h\right)\left({\lambda }_7-{\lambda }_5\right)+\left(\tau +{\mu }_v\right){\lambda }_7\hfill \end{array}\}$ (31)

which satisfies the transversality condition;

${\lambda }_1\left(tf\right)={\lambda }_2\left(tf\right)={\lambda }_3\left(tf\right)={\lambda }_4\left(tf\right)={\lambda }_5\left(tf\right)={\lambda }_6\left(tf\right)={\lambda }_7\left(tf\right)=0.$ (32)

By combining the Pontryagin’s Maximum Principle and the existence of optimal control [24] [25].

Theorem 7. The optimal control vector $\left(u_1^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}_2^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}_3^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}_4^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\right)$ that maximizes the objective function (J) over $\cup$, given by;

$\begin{array}{l}{u}_1^{*}\left(t\right)=\max \{0,\min (1,\frac{\beta \left({\lambda }_2-{\lambda }_1\right)\left(I_v+I_h\right)S_h^{*}}{2A_3}\hfill \\ +\frac{\lambda \left({\lambda }_5-{\lambda }_4\right)\left(I_v+I_h\right)S_v^{*}+b\lambda \left({\lambda }_5-{\lambda }_7\right)\left(I_v+I_h\right)V_v^{*}}{2A_3})\}\hfill \\ u_2^{*}\left(t\right)=\max \left\{0,\min \left(1,\frac{\left({\lambda }_2-{\lambda }_3\right)I_h^{*}}{2A_4}\right)\right\}\hfill \\ u_3^{*}\left(t\right)=\max \left\{0,\min \left(1,\frac{{\Lambda }_v\left({\lambda }_4-{\lambda }_7\right)}{2m{A}_5}\right)\right\}\hfill \\ u_4^{*}\left(t\right)=\max \left\{0,\min \left(1,\frac{\alpha \left({\lambda }_5-{\lambda }_6\right)I_v^{*}}{2A_6}\right)\right\}\hfill \end{array}\}$ (33)

where ${\lambda }_1\mathrm{,}{\lambda }_2\mathrm{,}{\lambda }_3\mathrm{,}{\lambda }_4\mathrm{,}{\lambda }_5\mathrm{,}{\lambda }_6$ and ${\lambda }_7$ are the solutions of Equation (31) and Equation (32).

Proof. Optimal control exists as a result the convexity of the integral of J w.r.t $u_1\mathrm{,}{u}_2\mathrm{,}{u}_3$ and $u_4$, Lipschitz property of the state system with respect to the state variables and a priori boundedness of the state solutions [9]. The system in (31) was obtained by differentiating the Hamiltonian function evaluated at optimal control. By equating the derivatives of Hamiltonian with respect to the controls to zero, we obtained the following;

$\begin{array}{l}{u}_1={\tilde{u}}_1:=\{(\frac{\beta \left({\lambda }_2-{\lambda }_1\right)\left(I_v+I_h\right)S_h^{*}}{2A_3}\hfill \\ +\frac{\lambda \left({\lambda }_5-{\lambda }_4\right)\left(I_v+I_h\right)S_v^{*}+b\lambda \left({\lambda }_5-{\lambda }_7\right)\left(I_v+I_h\right)V_v^{*}}{2A_3})\}\hfill \\ u_2={\tilde{u}}_2:=\left\{\left(\frac{\left({\lambda }_2-{\lambda }_3\right)I_h^{*}}{2A_4}\right)\right\}\hfill \\ u_3={\tilde{u}}_3:=\left\{\left(\frac{{\Lambda }_v\left({\lambda }_4-{\lambda }_7\right)}{2m{A}_5}\right)\right\}\text{and}\hfill \\ u_4={\tilde{u}}_4:=\left\{\left(\frac{\left({\lambda }_5-{\lambda }_6\right)I_v^{*}}{2A_6}\right)\right\}\hfill \end{array}\}$

By standard control arguments which involve bounds on the controls, it can be concluded that;

$\begin{array}{ccc}{u}_1^{*}=\{\begin{array}{ll}0,\hfill & \text{if}{\tilde{u}}_1\le 0\hfill \\ {\tilde{u}}_1\hfill & \text{if}0<{\tilde{u}}_1<1\hfill \\ 1\hfill & \text{if}{\tilde{u}}_1\ge 1\hfill \end{array}& ,& u_2^{*}=\{\begin{array}{ll}0,\hfill & \text{if}{\tilde{u}}_2\le 0\hfill \\ {\tilde{u}}_2\hfill & \text{if}0\prec <{\tilde{u}}_2<1\hfill \\ 1\hfill & \text{if}{\tilde{u}}_2\ge 1\hfill \end{array}\\ u_3^{*}=\{\begin{array}{ll}0,\hfill & \text{if}{\tilde{u}}_3\le 0\hfill \\ {\tilde{u}}_3\hfill & \text{if}0<{\tilde{u}}_3<1\hfill \\ 1\hfill & \text{if}{\tilde{u}}_3\ge 1\hfill \end{array}& ,& u_4^{*}=\{\begin{array}{ll}0,\hfill & \text{if}{\tilde{u}}_4\le 0\hfill \\ {\tilde{u}}_4\hfill & \text{if}0<{\tilde{u}}_4<1\hfill \\ 1\hfill & \text{if}{\tilde{u}}_4\ge 1\hfill \end{array}\end{array}\}$ (34)

The system in (34) leads to the system in (33) in Theorem (7). Uniqueness of optimal of the system is guaranteed by imposing a condition on time interval [26]. □

6. Numerical Results

The numerical simulations of the effects of control strategies on disease dynamics is shown. This is done by solving the optimal system consisting of Equation (3), co-state Equation (31), transversality conditions (32) and characterisation (34). Optimal system were solved by applying an iterative scheme. A fourth order Range-Kutta scheme was applied to solve the state equations with a guess of controls over time. We also apply the current iterations solutions of state systems to solve adjoint equations by backward fourth order Range-Kutta scheme. Finally, an update of controls by using a convex combination of previous controls and value from (33). This process was repeated and iterations stops if values of unknowns at the previous iterations are very close to the ones at the present iterations [27].

Most effective strategies were presented as follows: Combination of $u_1$ and $u_2$. Combination of $u_2$, $u_3$ and $u_4$. Combination of $u_2$ and $u_4$. Combination of $u_2$ and $u_3$. Combination of $u_1$ and $u_3$. Combination of $u_1$, $u_2$ and $u_3$. Combination of $u_3$ and $u_4$. Combination of $u_1$, $u_3$ and $u_4$. Combination of $u_1$ and $u_4$. The plots of the three most effective strategies were selected and presented as follows;

Strategy 1: Optimal treatment of infectious animals and treatment of humans.

Using $u_4$ and $u_2$, we optimised the objective functional, (J) by setting $u_1$ and $u_3$ to zero. Due to these control strategies applied, from Figure 3, it can be seen that the number of $I_v$ and $I_h$ have reduced substantially. This implies spread of anthrax can be curbed through effective treatment of infectious animals and the treatment of infectious humans. This strategy can best be reached by treatment all infectious animals and humans in the system. Figure 4 shows an exponential reduction in the number of recovered animals and humans the system.

Figure 3. Simulation of model indicating the effects of optimal strategies: $u_4$ and $u_2$.

Figure 4. Simulation of model indicating the effects of optimal strategies: $u_4$ and $u_2$.

Figure 5. Simulation of anthrax model indicating the effects of optimal strategies: $u_3$ and $u_1$.

Figure 6. Simulation of Anthrax model indicating the effects of optimal strategies: $u_3$ and $u_1$.

Figure 7. Simulation of model showing the effects of optimal strategies: $u_3$ and $u_4$.

Figure 8. Simulation of model indicating the effects of optimal strategies: $u_3$ and $u_4$.

Strategy 2: Optimal vaccination of animals and prevention of humans.

Using $u_3$ and $u_1$, we optimise J by setting the controls $u_2$ and $u_4$ to zero. As a result of this control strategies used, it can be seen that the number of $I_v$ and $I_h$ have reduced considerably as shown in Figure 5. This means that anthrax spread can be tackled through $u_3$ and $u_1$. This can be done though education of farmers on vaccination of animals against anthrax. Also, educating the public on the dangers associated with consumption of infected meat and products from animals infected with anthrax. Moreover, there are some reductions in the populations of $R_h$ and $R_v$ as shown in Figure 6.

Strategy 3: Optimal vaccination of animals and treatment of infectious animals.

Using $u_3$ and $u_4$, we optimise the objective functional by setting $u_1$ and $u_2$ to zero. As a result of this, it can be observed that there have been a reduction in the population of $I_v$ and $I_h$ as indicated in Figure 7. This means anthrax spread can handled by effective vaccination of susceptible animals and treatment of infectious animals. In order to achieve this, proper vaccination of animals and the treatment of all infectious animals should be conducted. Figure 8 shows a considerable reduction in both human and animal population.

7. Cost Effectiveness Analysis

In this section, a cost-benefit analysis is conducted to determine the costs associated with prevention of susceptible humans, treatment of infectious humans, vaccination of susceptible animals and treatment of infectious animals. In this paper, we consider the Infection Averted Ratio, (IAR) approach.

Intervention Averted Ratio (IAR)

The intervention averted ratio is defined as the ratio of the number of infection averted to the number of recovered. Number of infection averted is the diffeences between infectious population without optimal control and infectious population with optimal control. The IAR for each strategy was determined using the model’s parameter values. Generally, the strategy with the highest ratio is taken as the best and most effective strategy to be considered.

$\text{IAR}=\frac{\text{Number of infection averted}}{\text{Number of recovered}}$ (35)

The values in Table 3 shows the intervention averted ratios of the various strategies. Strategy 2 has the highest ratio and hence, it is the most effective strategy. The second most highest ratio is strategy 3. Hence, anthrax spread can be tackled through $u_1$ and $u_3$. This can be done though education of farmers on vaccination of animals against anthrax. Also, educating the public on the dangers associated with consumption of infected meat and products from animals infected with anthrax. The diagram in Figure 9 shows IAR plots on the effects of control strategies.

Figure 9. IAR plot showing control effects on strategies.

Table 3. Total infection averted ratio (IAR).

8. Conclusions

The model’s qualitative analysis exhibited an existence of multiple endemic equilibria. Its biological implications are that, proper anthrax disease can best be done if $R_{hv}<1$. Moreover, sensitivity analysis of $R_{hv}$ showed that, an increase in animal recovery rate would decrease $R_{hv}$. However, a decrease in animal recovery rate would increase $R_{hv}$. Additionally, an increase in animal and human transmission rates would increase $R_{hv}$ and a decrease in animal and human transmission rates decreases $R_{hv}$.

Analysis of optimal control and cost effectiveness of our model showed that the best and most effective strategy is vaccination of animals and prevention of susceptible humans in the system. Prevention control on humans and vaccination of animals should be considered as priority when fighting anthrax infections. There should be proper campaign on anthrax prevention and more animals should be vaccinated against the disease.

Data Availability Statement

Data supporting this model analysis are from previously published works and reported studies and all these have been cited in this paper. Some parameter values were assumed and others taken from published articles. These articles have also been cited accordingly at relevant places within the text as references.

Funding

This research is not funded by any institution.

Conflicts of Interest

Authors declare that there is no conflict of interest regarding the publication of this paper.

References

[1]	Osman, S. and Makinde, O.D. (2018) A Mathematical Model for Co-Infection of Listeriosis and Anthrax Diseases. International Journal of Mathematics and Mathematical Sciences, 2018, Article ID: 1725671. https://doi.org/10.1155/2018/1725671
[2]	Samad, M.A. (2013) Public Health Threat Caused by Zoonotic Diseases in Bangladesh. Bangladesh Journal of Veterinary Medicine, 9, 95-120. https://doi.org/10.3329/bjvm.v9i2.13451
[3]	Eustace, K.A., Osman, S. and Wainaina, M. (2018) Mathematical Modelling and Analysis of the Dynamics of Cholera. Global Journal of Pure and Applied Mathematics, 14, 1259-1275.
[4]	Muia, D.W., Osman, S. and Wainaina, M. (2018) Modelling and Analysis of Trypanosomiasis Transmission Mechanism. Global Journal of Pure and Applied Mathematics, 14, 1311-1331.
[5]	Osman, S., Musyoki, E.A. and Ndungu, R.M. (2019) A Mathematical Model for the Transmission of Measles with Passive Immunity. International Journal of Research in Mathematical and Statistical Sciences, 6, 1-8.
[6]	Karunditu, J.W., Kimathi, G. and Osman, S. (2019) Mathematical Modeling of Typhoid Fever Disease Incorporating Unprotected Humans in the Spread Dynamics. Journal of Advances in Mathematics and Computer Science, 32, 1-11. https://doi.org/10.9734/jamcs/2019/v32i330144
[7]	Gumel, A.B. and Moghadas, S.M. (2003) A Qualitative Study of a Vaccination Model with Non-Linear Incidence. Applied Mathematics and Computation, 143, 409-419. https://doi.org/10.1016/S0096-3003(02)00372-7
[8]	Kar, T.K. and Batabyal, A. (2011) Stability Analysis and Optimal Control of an Sir Epidemic Model with Vaccination. Biosystems, 104, 127-135. https://doi.org/10.1016/j.biosystems.2011.02.001
[9]	Makinde, O.D. and Okosun, K.O. (2011) Impact of Chemo-Therapy on Optimal Control of Malaria Disease with Infected Immigrants. BioSystems, 104, 32-41. https://doi.org/10.1016/j.biosystems.2010.12.010
[10]	Hattaf, K., Rachik, M., Saadi, S., Tabit, Y. and Yousfi, N. (2009) Optimal Control of Tuberculosis with Exogenous Reinfection. Applied Mathematical Sciences, 3, 231-240.
[11]	Lashari, A.A., Aly, S., Hattaf, K., Zaman, G., Jung, I.H. and Li, X.-Z. (2012) Presentation of Malaria Epidemics Using Multiple Optimal Controls. Journal of Applied Mathematics, 2012, Article ID: 946504. https://doi.org/10.1155/2012/946504
[12]	Rogers, D.J. (1988) A General Model for the African Trypanosomiases. Parasitology, 97, 193-212. https://doi.org/10.1017/S0031182000066853
[13]	Jones, P. and Martin, M. (2003) A Review of the Literature on the Occurrence and Survival of Pathogens of Animals and Humans in Green Compost. WRAP Standards Report. The Waste and Resources Action Programme, Oxon.
[14]	Okosun, K.O., Mukamuri, M. and Makinde, D.O. (2016) Global Stability Analysis and Control of Leptospirosis. Open Mathematics, 14, 567-585. https://doi.org/10.1515/math-2016-0053
[15]	Van den Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Mathematical Biosciences, 180, 29-48. https://doi.org/10.1016/S0025-5564(02)00108-6
[16]	LaSalle, J.P. (1976) The Stability of Dynamical Systems (Regional Conference Series in Applied Mathematics). SIAM, Philadelphia.
[17]	Kanyaa, J.K., Osman, S. and Wainaina, M. (2018) Mathematical Modelling of Substance Abuse by Commercial Drivers. Global Journal of Pure and Applied Mathematics, 14, 1149-1165.
[18]	Makinde, O.D. (2007) Adomian Decomposition Approach to a Sir Epidemic Model with Constant Vaccination Strategy. Applied Mathematics and Computation, 184, 842-848. https://doi.org/10.1016/j.amc.2006.06.074
[19]	Osman, S., Makinde, O.D. and Theuri, D.M. (2018) Stability Analysis and Modelling of Listeriosis Dynamics in Human and Animal Populations. Global Journal of Pure and Applied Mathematics, 14, 115-137. https://doi.org/10.1155/2018/1725671
[20]	Eyaran, W.E., Osman, S. and Wainaina, M. (2019) Modelling and Analysis of Seir with Delay Differential Equation. Global Journal of Pure and Applied Mathematics, 15, 365-382.
[21]	Lashari, A.A. and Zaman, G. (2012) Optimal Control of a Vector Borne Disease with Horizontal Transmission. Nonlinear Analysis: Real World Applications, 13, 203-212. https://doi.org/10.1016/j.nonrwa.2011.07.026
[22]	Okosun, K.O. and Makinde, O.D. (2014) Optimal Control Analysis of Hepatitis c Virus with Acute and Chronic Stages in the Presence of Treatment and Infected Immigrants. International Journal of Biomathematics, 7, Article ID: 1450019. https://doi.org/10.1142/S1793524514500193
[23]	Okosun, K.O. and Makinde, O.D. (2014) A Co-Infection Model of Malaria and Cholera Diseases with Optimal Control. Mathematical Biosciences, 258, 19-32. https://doi.org/10.1016/j.mbs.2014.09.008
[24]	Pontryagin, L.S. (1987) Mathematical Theory of Optimal Processes. CRC Press, Boca Raton.
[25]	Fleming, W.H. and Rishel, R.W. (2012) Deterministic and Stochastic Optimal Control, Volume 1. Springer Science & Business Media, Berlin.
[26]	Okosun, K.O., Ouifki, R. and Marcus, N. (2011) Optimal Control Analysis of a Malaria Disease Transmission Model That Includes Treatment and Vaccination with Waning Immunity. Biosystems, 106, 136-145. https://doi.org/10.1016/j.biosystems.2011.07.006
[27]	Joshi, H.R., Lenhart, S., Li, M.Y. and Wang, L.C. (2006) Optimal Control Methods Applied to Disease Models. Contemporary Mathematics, 410, 187-208. https://doi.org/10.1090/conm/410/07728