Stability and Optimal Control of Tuberculosis Spread with an Imperfect Vaccine in the Case of Co-Infection with HIV

Leontine Nkague Nkamba^1,2, Thomas Timothee Manga³, Noboru Sakamoto^2,4

¹Department of Mathematics, Higher Teacher Training College, University of Yaoundé I, Yaoundé, Cameroon.
²Deustotech Laboratory Chair of Computational Mathematics, University of Deusto, Bilbao, Spain.
³AIDEPY Association des Ingénieurs Diplomés de l’Ecole Polytechnique de Yaoundé, Yaoundé, Cameroon.
⁴Faculty of Science and Engineering, Department of Mechatronics, Nanzan University, Nagoya, Japan.

DOI: 10.4236/ojmsi.2019.72005                

Abstract

This paper focuses on the study and control of a non-linear mathematical epidemic model ( SS_vihVELI ) based on a system of ordinary differential equation modeling the spread of tuberculosis infectious with HIV/AIDS coinfection. Existence of both disease free equilibrium and endemic equilibrium is discussed. Reproduction number R₀ is determined. Using Lyapunov-Lasalle methods, we analyze the stability of epidemic system around the equilibriums (disease free and endemic equilibrium). The global asymptotic stability of the disease free equilibrium whenever R_{vac < 1} is proved, where R₀ is the reproduction number. We prove also that when R₀ is less than one, tuberculosis can be eradicated. Numerical simulations are conducted to approve analytic results. To achieve control of the disease, seeking to reduce the infectious group by the minimum vaccine coverage, a control problem is formulated. The Pontryagin’s maximum principle is used to characterize the optimal control. The optimality system is derived and solved numerically using the Runge Kutta fourth procedure.

Keywords

Tuberculosis, Basic Reproduction Number, Global Stability, Prevalence, CD4 Cells, Immune Deficiency, Optimal Control, Optimality, HIV

1. Introduction

Tuberculosis is one of the top 10 causes of death in the world [1]. In Cameroon, tuberculosis and HIV infection account for 25% of the morbid burden. TB is caused by bacteria (Mycobacterium tuberculosis) that most often affects the lungs. Tuberculosis is curable and preventable. TB is spread from person to person through the air. When people with lung TB cough, sneeze or spit, they propel the TB germs into the air. A person needs to inhale only a few of these germs to become infected. About one-third of the world population has latent TB, which means people have been infected by TB bacteria but are not (yet) ill with the disease and cannot transmit the disease. People infected with TB bacteria have a 10% lifetime risk of falling ill with TB. However, persons with compromised immune systems, such as people living with HIV, malnutrition or diabetes, or people who use tobacco, have a much higher risk of falling ill. When a person develops active TB disease, the symptoms (such as cough, fever, night sweats, or weight loss) may be mild for many months. This can lead to delays in seeking care, and results in transmission of the bacteria to others. People with active TB can infect 10 other people through close contact over the course of a year. The population of Cameroon is estimated in 2017 to 24 million inhabitants. In the same year, 24,905 cases of tuberculosis were reported. The medical coverage is estimated at 52%. The death rate due to the tuberculose is approximatively 0.29% [2]. In 2018, the national budget allocated to the control of the tuberculosis is 10 million American dollars, and distributed as follows: 14% of internal source, 27% of external source and 59% of deficit.

The theme “End the global TB” has been decided by the world health organization and it covers the period 2016-2035 and the overall goal is to end the global TB epidemic. As part of the necessary multidisciplinary research approach, mathematical models have been extensively used to provide a framework for understanding tuberculosis transmission dynamics and control strategies of the infection spread in the host population [3] [4] [5] [6] [7]. Arthur et al. [6] study social and cultural factors in the successful control of tuberculosis. Optimal control theory has been applied to some biological models (VIH/AIDS treatment, Cancer, Fish harvestime) [8] - [14] and transmission model.

Some mathematical models have been used to control tuberculosis. In [15] , optimal control theory is applied to investigate optimal strategies for controlling the spread of tuberculosis using treatment of infected individuals with TB as the system control variables. In [16] , the tuberculosis control is formulated and solved as an optimal control theory problem, indicating how a control term on the chemoprophylaxis should be introduced in the population to reduce the number of individuals with active TB. The papers [17] consider the optimal control of tuberculosis through education, diagnosis campaign and chemoprophylaxis of latently infected. In [6] , it talks about the Impact of an effective multidrug-resistant tuberculosis control programme in the setting of an immature HIV epidemic.

This paper deals with the stability analysis of an $SV{S}_{vih}\mathrm{,}ELI$ transmission model and uses optimal control technique to find and evaluate the impact of a mass vaccination schedule in the spread of TB/VIH coinfection. Individuals are classified as one of susceptible (S) V vaccinated, AIDS patients ( $S_{vih}$ ), earlier latent (E) late latent (L) infectious or tuberculous (I), but allow that susceptible individuals may be given an imperfect vaccine that reduces their susceptibility to the disease, the V-compartment of vaccinated individuals is considered [18] as a susceptible compartment. AIDS patients are considered as another type of susceptible with hight level of susceptibility due to the fact that their immune system are compromised.

Since V and $S_{vih}$ are considered as the susceptible compartments, thus we are dealing with a differential susceptibility system with bilinear mass action as in Hyman and Li [19]. However, we include one-way flow between these two compartments, to denote the prevalence of AIDS in community. For analysis and stability issues;We use Lyapunov-LaSalle methods, we fully resolve the global dynamics of the model for the full parameter space. We demonstrate that the model exhibits threshold behavior with a globally stable disease-free equilibrium if the basic reproduction number is less than unity and a globally stable endemic equilibrium if the basic reproduction number is greater than unity. Our goal is to minimise number of TB infectious persons and the overall cost of the vaccine during a fixed period.

For numerical simulation, the programs used in this paper are designed so that no knowledge of MATLAB is required. For the control problem, there is a user-friendly interface that will guide you through. We have two different MATLAB programs, plotTB.m and codeTB.m. Here, .m is the extension given to all files intended for use in MATLAB.The file codeTB.m is the Runge-Kutta based, forward-backward sweep solver. It takes as input the values of the various parameters in the problem and outputs the solution to the optimality system. The ?le plotTB.m is the user friendly interface. It will ask you to enter the values of the parameters one by one, compile codeTB.m with these values, and plot the resulting solutions. All the files must be in the directory that MATLAB treats as the home directory. This is usually the Work directory. This paper is organized as follows.

In next section, model is described, in Section 3 we investigate stability analysis for the ( $S{S}_{vih}VELI$ ) epidemic model in this section We followed the methods of Nkamba, Leontine Nkague et al. 2019 [20] A control system for the optimality and its existence,and the optimal control are derived in Section 4. In Section 5, utilizing the representation of the optimal control, we describe a numerical solution of the optimality system consisting of the original state system, the adjoint system, and their boundary conditions. In Section 6, we describe, in detail, a real application of our optimal control theory. Finally, we conclude in Section 7.

2. Model Description

2.1. Basic Model of $S{S}_{vih}VELI$,Preliminaries

When first infected with TB bacteria, a person typically goes through a latent, asymptomatic and non-infectious period during which the body’s immune system fights the TB bacteria. There are two distinct stages of the latent TB infection. During the first two years, the risk of developing active disease is much higher, whereas during the later stage, the progression to active disease is much slower.

Compartmental modeling is used among epidemiologists to simulate disease dynamics. These models treat each disease state as a different compartment that contains a homogeneous population of individuals. Using a compartmental approach, the total host population can be partitioned into seven compartments: susceptible individuals (S), susceptible infected with AIDS ( $S_{vih}$ ), vaccinated individuals (V) early latent (E) late latent (L) individuals, individuals with active TB disease (I) and recovered individuals (R). $S\left(t\right)\mathrm{,}{S}_{vih}\left(t\right)\mathrm{,}V\left(t\right)\mathrm{,}E\left(t\right)\mathrm{,}L\left(t\right)$ and $I\left(t\right)$ denote the density of populations in the four corresponding compartments at time t. Only individuals in compartment I are infectious, and new infections result on the one hand from contacts between a susceptible and an infectious individual, with an incidence rate $\beta S\left(t\right)I\left(t\right)$; from contacts between a HIV patient and Tuberculous, with an incidence $\left(1+\sigma \right)\beta S_{vih}\left(t\right)I\left(t\right)$,and on the other hand from contacts between a vaccinated and infectious individual, with an incidence ${\theta }_1\beta V\left(t\right)I\left(t\right)$,due to the fact that the vaccine does not confer a total immunity, but Vaccination reduces the risk of infection by a factor ${\theta }_1\in \left[\mathrm{0,1}\right]$ and the efficacy of the vaccine is $1-{\theta }_1$. AIDS increases the risk of infection by a factor, the immune deficiency rate $\sigma \in \left[\mathrm{0,1}\right]$ $\sigma =1-T_{CD4}$ $T_{CD4}$ is the rate of CD4 cells. Let us pose ${\theta }_2=\left(1+\sigma \right)$ the cost induced by the immune deficiency status in the transmission of tuberculosis. The per capita death rates for susceptibles, HIV patients, early latents, late latents and infectious individuals are ${\mu }_S$, ${\mu }_{vih}$, ${\mu }_E$, ${\mu }_L$ and ${\mu }_I$ respectively. Once infected, individuals progress through the early latent stage with an average rate $\omega$. A fraction $p;0<p\le 1$; of these individuals progress directly to the active TB stage, and the remaining $1-p$ fraction progresses to the late latent stage. Once there, the rate of progression to active disease is at a lower rate $\nu$. The recruitment makes respectively into the susceptible class, the vaccinated class, the VIH/patient class with the constant rate ${\pi }_1$ ${\pi }_2$ and ${\pi }_3$. $\alpha$ is the vaccination coverage rate. ${\pi }_2$ is the recruitment of vaccinated a few day after birth, so we suppose that immunity is passed during the birth. ${\pi }_3$ is VIH/AIDS vertical transmission recruitment, it’s means some peoples born with VIH/AIDS infection.

The dynamical transfer among the seven compartments is depicted in the transfer diagram (Figure 1).

2.2. Description of Variables and Parameters

All parameters described in Table 1 are assumed to be positive.

2.3. Compartmental Diagram and Differential Equations of $S{S}_{vih}VELT$ Model

Our model consists of the following system of ordinary differential equations:

Figure 1. This is a flow chart for our model. The seven boxes represent the seven groups of individuals. The arrows show the movement between groups, and into and out of the population.

Table 1. Description of variables and parameters.

$\{\begin{array}{l}\stackrel{\dot{}}{S}={\pi }_1-\beta SI-\left({\mu }_S+\epsilon +\alpha \right)S,\\ \stackrel{\dot{}}{V}={\pi }_2+\alpha S-{\theta }_1\beta VI-{\mu }_{V}V\\ {\stackrel{\dot{}}{S}}_{vih}={\pi }_3+\epsilon S-{\mu }_{vih}{S}_{vih}-{\theta }_2\beta S_{vih}I,\\ \stackrel{\dot{}}{E}=\beta SI+{\theta }_2\beta S_{vih}I+{\theta }_1\beta VI-\left({\mu }_E+\omega \right)E,\\ \stackrel{\dot{}}{L}=\left(1-p\right)\omega E-\left({\mu }_L+\nu \right)L,\\ \stackrel{\dot{}}{I}=p\omega E+\nu L-\left({\mu }_I+\tau \right)I,\\ \stackrel{\dot{}}{R}=\tau I-{\mu }_R\\ N=S+S_{vih}+V+E+L+I+R\end{array}$ (1)

with initial conditions $\left(S\left(0\right)\mathrm{,}{S}_{vih}\left(0\right)\mathrm{,}V\left(0\right)E\left(0\right)\mathrm{,}L\left(0\right)\mathrm{,}I\left(0\right)R\left(0\right)\right)\in {ℝ}_{+}^5$. We have also ${\theta }_2=1+\sigma$ and $\sigma =1-T_{cd4}$.

3. Stability Analysis

3.1. A Compact Positively Invariant Absorbing Set

In order that the model be well-posed, it is necessary that the state variables $S\left(t\right)$, $S_{vih}\left(t\right)$, $V\left(t\right)$ $E\left(t\right)$, $R\left(t\right)$ $I\left(t\right)$ and $I\left(t\right)$ remain nonnegative for all $t\ge 0$. That is, the nonnegative orthant ${ℝ}_{+}^5$ must be positively invariant.

Let

$\Gamma =\left\{\left(S\mathrm{,}{S}_{vih}\mathrm{,}V\mathrm{,}E\mathrm{,}L\mathrm{,}I\right)\in {ℝ}_{+}^5\mathrm{<mo>\; :\; </mo>0}\le S\le S_0\mathrm{,0}\le S_{vih}\le S_{vih0}\mathrm{,0}\le V\le V_0\mathrm{,0}\le N\le \frac{\Lambda }{\mu }\right\}$ (2)

where $\mu =\min \left\{{\mu }_S,{\mu }_{vih}{\mu }_V,{\mu }_E,{\mu }_L,{\mu }_I\right\}$ and $\Lambda ={\pi }_1+{\pi }_2+{\pi }_3$ with $\left(S_0{V}_0{S}_{vih0}0000\right)$ the disease free equilibrium.

Lemma 1. The compact set $\Gamma$ is a positively invariant and attracting.

3.2. Existence of an Disease Free Equilibrium (DFE)

It is easy to check that model 1 always has the disease-free equilibrium

$P_0=\left(S_0,S_{vih0},V_0,0,0,0,0\right)$

where

$S_0=\frac{{\pi }_1}{{\mu }_S+\epsilon +\alpha },S_{vih0}=\frac{\epsilon S_0+{\pi }_3}{{\mu }_{vih}}\text{and}{V}_0=\frac{\alpha S_0+{\pi }_2}{{\mu }_V}$ (3)

In order to assume that vaccinated people don’t produced infected more than susceptible people we should have $0\le {\theta }_1\le \frac{S_0}{V_0}$.

Additionally, an endemic equilibrium $P^{*}=\left(S^{*},S_{vih}^{*},V^{*}{E}^{*},L^{*},I^{*}\right)$ may also exist.

To consider the existence and uniqueness of endemic equilibrium $P^{*}=\left(S^{*},S_{vih}^{*},V^{*}{E}^{*},L^{*},I^{*}\right)$,we firstly study the basic reproductive number induced by vaccine $R_{vac}$ of model.

3.3. Basic Reproduction Ratio

Using the method of James Watmouth and all the next generation matrix [21] [22] , the basic reproduction number ${\mathcal{R}}_0$ is giving by

${\mathcal{R}}_0=\frac{\beta \omega \left(S_0+{\theta }_1{V}_0+{\theta }_2{S}_{vih0}\right)\left[\nu \left(1-p\right)+p\left({\mu }_L+\nu \right)\right]}{\left({\mu }_I+\tau \right)\left({\mu }_E+\omega \right)\left({\mu }_L+\nu \right)}\mathrm{.}$ (4)

We will see in the Section 3.4 theorem 1 that, when ${\mathcal{R}}_0$ is less than unity, infection can disappear in the population. Numerical simulations will confirm our results.

3.4. Stability of Equilibriums

Stability of Disease Free Equilibrium (DFE)

In this section, we show that the disease-free equilibrium $P_0$ is globally asymptotically stable with respect to if $R_0\le 1$; and $P_0$ is unstable if $R_0>1$ :

Theorem 1. If ${\mathcal{R}}_0\le 1$,then the disease-free equilibrium is globally asymptotically stable.

Proof. Consider a Lyapunov function,

$\mathcal{V}=\mathcal{V}\left(S\mathrm{,}{S}_{vih}\mathrm{,}V\mathrm{,}E\mathrm{,}L\mathrm{,}I\right)=p\omega \left({\mu }_L+\nu \right)E+\nu \left({\mu }_E+\omega \right)L+\left({\mu }_E+\omega \right)\left({\mu }_L+\nu \right)I\mathrm{.}$ (5)

Direct calculation leads to

$\begin{array}{c}\stackrel{\dot{}}{\mathcal{V}}=p\omega \left({\mu }_L+\nu \right)\stackrel{\dot{}}{E}+\nu \left({\mu }_E+\omega \right)\stackrel{\dot{}}{L}+\left({\mu }_E+\omega \right)\left({\mu }_L+\nu \right)\stackrel{\dot{}}{I}\\ =p\omega \left({\mu }_L+\nu \right)\left(\beta SI+{\theta }_1\beta VI+{\theta }_2\beta S_{vih}\right)I+\nu \left(1-p\right)\left({\mu }_E+\omega \right)E\\ +\left({\mu }_E+\omega \right)\left({\mu }_L+\nu \right)\left({\mu }_I+\tau \right)I\end{array}$ (6)

at equilibrium we have the relation

$\beta SI+{\theta }_2\beta S_{vih}I+{\theta }_1\beta VI=\left({\mu }_E+\omega \right)E$

then

$\begin{array}{l}\stackrel{\dot{}}{\mathcal{V}}=\left({\mu }_E+\omega \right)\left({\mu }_L+\nu \right)\left({\mu }_I+\tau \right)\\ \times \left[\frac{\left(\beta S+{\theta }_2\beta S_{vih}+{\theta }_1\beta V\right)\left[p\omega \left({\mu }_L+\nu \right)+\nu \omega \left(1-p\right)\right]}{\left({\mu }_E+\omega \right)\left({\mu }_L+\nu \right)\left({\mu }_I+\tau \right)}-1\right]I\end{array}$ (7)

Because $S\mathrm{,}V\mathrm{,}{S}_{vih}\in \Gamma$ we have:

$\begin{array}{l}\stackrel{\dot{}}{\mathcal{V}}\le \left({\mu }_E+\omega \right)\left({\mu }_L+\nu \right)\left({\mu }_I+\tau \right)\\ \times \left[\frac{\left(\beta S_0+{\theta }_2\beta S_{vih0}+{\theta }_1\beta V_0\right)\left[p\omega \left({\mu }_L+\nu \right)+\nu \omega \left(1-p\right)\right]}{\left({\mu }_E+\omega \right)\left({\mu }_L+\nu \right)\left({\mu }_I+\tau \right)}-1\right]I\\ \stackrel{\dot{}}{\mathcal{V}}\le \left({\mu }_E+\omega \right)\left({\mu }_L+\nu \right)\left({\mu }_I+\tau \right)\left[{\mathcal{R}}_{vac}-1\right]I\mathrm{.}\end{array}$ (8)

Furthermore

$\stackrel{\dot{}}{\mathcal{V}}=0\iff I=0\Rightarrow \stackrel{\dot{}}{I}=0\Rightarrow L=E=0$

Therefore, the largest compact invariant set in $G=\left\{\left(S;S_{vih};V;E;L;I\right)\in \Gamma :\stackrel{\dot{}}{\mathcal{V}}=0\right\}$; when ${\mathcal{R}}_{vac}\le 1$,is the singleton $\left\{P_0\right\}$ : LaSalle’s Invariance Principle implies that all solutions in $\Gamma$ converges to $P_0$. This establishes the theorem.

Theorem 1 completely determines the global dynamics of (1) in $\Gamma$ when $R_0\le 1$. It establishes the basic reproduction number $R_0$ as a sharp threshold parameter. Namely, if $R_0<1$; all solutions in the feasible region converge to the disease-free equilibrium $P_0$; and the TB will die out from the population irrespective of the initial conditions. If $R_0>1$; $P_0$ is unstable, it could exist and endemic equilibrium and the system is uniformly persistent, TB epidemic will always become endemic.

3.5. Numerical Simulation of DFE and Endemic Equilibrium

In matlab code, the parameters are named as followed:

pi1 = π₁ pi2 = π₂ pi3 = π₁; beta = β varepsilon = ε; muvih = μ_vih; sigma = σ; theta1 = θ₁; theta2 = θ₂ p; muS = μ_S muV = μ_V muE = μ_E omega = ω muL = μ_L nu = ν muI = μ_I tau = τ A; alpha = α S0 = S₀ V0 = S₀ Svih0 = S₀ E0 = E₀ L0 = L₀ I0 = I₀.

Let us take the following set of parameters

pi1 = 15; pi2 = 0.1; pi3 = 0.1; beta = 0.085; varepsilon = 0.2; muvih = 0.05; sigma = 0.01; theta1 = 0.001; theta2 = 1.01; p = 0.2; muS = 0.01; muV = 0.01; muE = 0.02; omega = 0.0645; muL = 0.02; nu = 0.00375; muI = 0.3; tau = 0.5; alpha = 0.2; S0 = 80; V0 = 50; Svih0 = 10; E0 = 20; L0 = 15; I0 = 25.

When $\beta =0.01$ and $R_0=0.57$ $R_0$ is less than unity, the trajectory of TB patients (I) reach axis axe by 10 years and remains at that position with time (Figure 2) black line.

We remark also that when $\beta =0.085$ $R_0=4.58$ greater than unity, the TB patients remains in the community, so it could exist an endemic equilibrium stable in Figure 2 of the blue line.

4. Optimal Control of System 1

In this section, an optimal control is formulated and it examined to study properties of optimal control strategies.

Let $u_1\left(t\right)=\alpha \left(t\right)$ (vaccine coverage), the first control, be the percentage of susceptible individuals being vaccinated per unit of time. As vaccination of the entire susceptible population is impossible, we bound the control with $0\le u\left(t\right)\le 0.9$. In cash we seek to minimise the infectious group with the minimum possible of vaccine coverage. We consider an optimal control problem

Figure 2. Dynamic and stability of system 1.

to minimize the objective functional

$\underset{u}{\min }{\displaystyle {\int }_1^T}I+Au{\left(t\right)}^2\text{d}x$ (9)

$\begin{array}{l}\text{subject}\text{to}\\ \stackrel{\dot{}}{S}\left(t\right)={\pi }_1-\beta S\left(t\right)I\left(t\right)-\left({\mu }_S+u\left(t\right)+\epsilon \right)S,S\left(0\right)=S_0\\ \stackrel{\dot{}}{V}\left(t\right)={\pi }_2+u_1\left(t\right)S\left(t\right)-{\theta }_1\beta V\left(t\right)I\left(t\right)-{\mu }_{V}V\left(t\right),V\left(0\right)=S_0\\ {\stackrel{\dot{}}{S}}_{vih}\left(t\right)=\left({\pi }_3+\epsilon \right)S-{\mu }_{vih}{S}_{vih}\left(t\right)-{\theta }_2\beta S_{vih}\left(t\right)I\left(t\right),S_{vih}\left(0\right)=S_{vih0}\end{array}$

$\begin{array}{l}\stackrel{\dot{}}{E}\left(t\right)=\beta S\left(t\right)I\left(t\right)+{\theta }_2\beta S_{vih}\left(t\right)I\left(t\right)+{\theta }_1\beta V\left(t\right)I\left(t\right)-\left({\mu }_E+\omega \right)E\left(t\right),E\left(0\right)=E_0\\ \stackrel{\dot{}}{L}\left(t\right)=\left(1-p\right)\omega E\left(t\right)-\left({\mu }_L+\nu \right)L\left(t\right),L\left(0\right)=L_0\\ \stackrel{\dot{}}{I}\left(t\right)=p\omega E\left(t\right)+\nu L\left(t\right)-\left({\mu }_I+\tau \right)I\left(t\right),I\left(0\right)=I_0\\ 0\le u\left(t\right)\le 0.8\end{array}$ (10)

Let us pose

$\begin{array}{l}\lambda \left(t\right)=\left({\lambda }_i\left(t\right)\right)\text{with}i=1,2,\cdots ,6\\ X\left(t\right)=\left(X_1,X_2,X_3,\cdots ,X_6\right)\end{array}$ (11)

Hamiltonian of our control problem is

$H\left(tX\left(t\right)u\left(t\right)\lambda \left(t\right)\right)=I+A{u}_1{\left(t\right)}^2+{\displaystyle \underset{i=1}{\overset{6}{\sum }}}{\lambda }_i\left(t\right){\stackrel{\dot{}}{X}}_i(t)$

$\begin{array}{l}H\left(tX\left(t\right)u\left(t\right)\lambda \left(t\right)\right)\\ =I+Au{\left(t\right)}^2+{\lambda }_1\left[{\pi }_1-\beta S\left(t\right)I\left(t\right)-\left({\mu }_S+u\left(t\right)+\epsilon \right)S\right]\\ +{\lambda }_2\left[{\pi }_2+u\left(t\right)S\left(t\right)-{\theta }_1\beta V\left(t\right)I\left(t\right)-{\mu }_{V}V\left(t\right)\right]\\ +{\lambda }_3\left[\left({\pi }_3+\epsilon \right)S-{\mu }_{vih}{S}_{vih}\left(t\right)-{\theta }_2\beta S_{vih}\left(t\right)I\left(t\right)\right]\end{array}$

$\begin{array}{l}+{\lambda }_4\left[\beta S\left(t\right)I\left(t\right)+{\theta }_2\beta S_{vih}\left(t\right)I\left(t\right)+{\theta }_1\beta V\left(t\right)I\left(t\right)-\left({\mu }_E+\omega \right)E\left(t\right)\right]\\ +{\lambda }_5\left[\left(1-p\right)\omega E\left(t\right)-\left({\mu }_L+\nu \right)L\left(t\right)\right]\\ +{\lambda }_6\left[p\omega E\left(t\right)+\nu L\left(t\right)-\left({\mu }_I+\tau \right)I\left(t\right)\right]\end{array}$ (12)

The adjoint equations and transversality conditions can be obtained by using Pontryagin’s Maximum Principle such that

$\begin{array}{l}-\frac{\partial H}{\partial S}={{\lambda }^{\prime }}_1=\left({\lambda }_1-{\lambda }_4\right)\beta I^{*}+{\lambda }_1\left({\mu }_S+u+\epsilon \right)-{\lambda }_{2}u-{\lambda }_3\epsilon ,{\lambda }_1\left(t_f\right)=0\\ -\frac{\partial H}{\partial V}={{\lambda }^{\prime }}_2=\left({\lambda }_2-{\lambda }_4\right){\theta }_1\beta I^{*}+{\lambda }_2{\mu }_V,{\lambda }_2\left(t_f\right)=0\\ -\frac{\partial H}{\partial S_{vih}}={{\lambda }^{\prime }}_3=\left({\lambda }_3-{\lambda }_4\right){\theta }_2\beta I^{*}+{\lambda }_3{\mu }_{vih},{\lambda }_3\left(t_f\right)=0\\ -\frac{\partial H}{\partial E}={{\lambda }^{\prime }}_4={\lambda }_4\left({\mu }_E+\omega \right),{\lambda }_4\left(t_f\right)=0\end{array}$

$\begin{array}{l}-\frac{\partial H}{\partial L}={{\lambda }^{\prime }}_5={\lambda }_5\left({\mu }_L+\nu \right),{\lambda }_5\left(t_f\right)=0\\ -\frac{\partial H}{\partial I}={{\lambda }^{\prime }}_6=-1+{\lambda }_6\left({\mu }_I+\tau \right)+\beta \left({\lambda }_1-{\lambda }_4\right)S+\beta {\theta }_1\left({\lambda }_2-{\lambda }_4\right)V\\ +\beta {\theta }_2\left({\lambda }_3-{\lambda }_4\right)S_{vih},{\lambda }_6\left(t_f\right)=0\end{array}$ (13)

Taking into account the bounds on control $0\le u\le 0.8$ Optimal control $u^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)$,is derived using the following optimality conditions:

$\{\begin{array}{l}u=0\text{if}\frac{\partial H}{\partial u}\ge 0\\ 0\le u\le 0.8\text{if}\frac{\partial H}{\partial u}=0\\ u=0.8\text{if}\frac{\partial H}{\partial u}\le 0\end{array}$ (14)

From relations 14 we have:

$\frac{\partial H}{\partial u}=2A{u}^{*}+S^{*}\left({\lambda }_2-{\lambda }_1\right)=0\text{when}0\le u^{*}\le 0.8$ (15)

$\frac{\partial H}{\partial u}=0\iff 2A{u}^{*}+S^{*}\left({\lambda }_2-{\lambda }_1\right)=0$ (16)

From relations 16 we have:

$u^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)=\frac{S^{\mathrm{<mo>\; *\; </mo>}}\left({\lambda }_1-{\lambda }_2\right)}{2A}\text{where}0\le u^{\mathrm{<mo>\; *\; </mo>}}\le 0.8$ (17)

Taking account the optimality conditions 14 induce by the bounds conditions of control u we have

$u^{*}\left(t\right)=\min \left(0.8,\max \left(0,\frac{S^{*}\left({\lambda }_1-{\lambda }_2\right)}{2A}\right)\right)$ (18)

5. Numerical Simulations of the Spread of Tuberculosis/VIH Coinfection and Mass Vaccination Schedule

5.1. The Numerical Algorithm

The numerical algorithm presented below is a classical Rung -Kutta 4 method.

We discretize the interval [t0,T = tf] at the points $ti=t0+ih$ ( $i=0,1,\cdots ,n$ ), where h is the time step such that $tn=T=M$, $h2=h/2$ and $j=M+2-i$.

Next, we define the state and adjoint variables respectivily $S\left(t\right)$, $S_{vih}\left(t\right)$, $V\left(t\right)$ $E\left(t\right)$, $R\left(t\right)$ $I\left(t\right)$ and ${\lambda }_1\left(t\right)$, ${\lambda }_2\left(t\right)$, ${\lambda }_3\left(t\right)$, ${\lambda }_4\left(t\right)$, ${\lambda }_5\left(t\right)$, ${\lambda }_6(t))$

The control u in terms of nodal points $S\left(i\right)$, $S_{vih}\left(i\right)$, $V\left(i\right)$, $E\left(i\right)$, $R\left(i\right)$, $I\left(i\right)$ and ${\lambda }_1\left(j\right)$, ${\lambda }_2\left(j\right)$, ${\lambda }_3\left(j\right)$, ${\lambda }_4\left(j\right)$, ${\lambda }_5\left(j\right)$, ${\lambda }_6\left(j\right)$ $u\left(i\right)$.

5.2. Parameters Description and Values

Most of parameters values are from Cameroon, like natural rate of mortality and rate of birth according to the World Health organisation report 2017. Other parameters are extracted to the data collected at the Hospital Jamot Center, where is housed the screening center of tuberculosis. This center receive approximatively one thousand new cases of tuberculosis each year. Data have been collected during one year (31 March 2016 to 31 March 2017). Those data concerned new cases of pulmonary and extrapulmonary TB patients; number of TB patients tested HIV; Number of TB patients co-infected by VIH; The distribution of TB patients by 11 age and sex. Helped by theses data we found out that, about three new cases are detected by day, the mortality date is 0.1; the percentage of TB patients tested for HIV is 90% and about 35% of them are HIV positive.The mortality rate due to the infection is 0.1% Recovery rate and rate of apparition of clinical symptoms are coming from WHO. See Table 2 for the description of parameters and their based line or range value.

5.3. Optimal Strategies of Mass Vaccination

Some key parameters like the effective contact rate $\beta$ medical coverage rate of VIH/Patients $\sigma$,HIV vertical transmission rate ${\pi }_3$ and HIV prevalence $\epsilon$ have a great impact in the spread of TB infectious. we are going to simulate five different scenarios and observed the mass vaccination optimal strategy induced respectively by the low effective contact rate, the hight effective contact rate, the absence of HIV medical coverage, the hight HIV medical coverage and the combination (hight HIV medical coverage, low HIV vertical transmission rate, and the low HIV prevalence rate) which assure the eradication of TB infectious.

Table 2. Description and values of parameters.

For simulations we have built two matlab file:

The first one is tuberculosesida.m the main file and the second one is val_t uberculosesida.m in order to plot some figures.

In matlab code, the parameters are named as followed.

pi1 = π₁; pi2 = π₂ pi3 = π₁; beta = β varepsilon = ε; muvih = μ_vih; sigma = σ; theta1 = θ₁; theta2 = θ₂ p; muS = μ_S; muV = μ_V; muE = μ_E; omega = ω; muL = μ_L; nu = ν; muI = μ_I; tau = τ A; alpha = α; S0 = S₀; V0 = S₀; Svih0 = S₀; E0 = E₀; L0 = L₀; I0 = I₀.

5.3.1. Optimal Control Induced by the Low Effective Contact Rate $\beta =0.08$

Set of parameters values

pi1 = 50; pi2 = 1; pi3 = 0.1; beta = 0.08; varepsilon = 0.00001; muvih = 0.05; sigma = 0.0001; theta1 = 0.3; theta2 = 1.0001; p = 0.6; muS = 0.01; muV = 0.01; muE = 0.02; omega = 0.0605; muL = 0.02; nu = 0.00375; muI = 0.12; tau = 0.1; A = 1; alpha = 0.2; S0 = 80;V0 = 50; Svih0 = 10; E0 = 8; L0 = 7; I0 = 5.

With the precedents values of parameters we obtain a Reproduction number ( $R_0$ ) bigger than unity. In this case Figure 3, the system will reach definitively and endemic equilibrium. Controlling the system is not necessary until the trajectory of TB patients approach the endemic equilibrium and then optimal control appears in step 110 and goes until step 200 in order to have a better endemic equilibrium.

5.3.2. Optimal Control Induced by A High Effective Contact Rate $\beta =0.6$

Set of parameters values

pi1 = 50; pi2 = 1; pi3 = 0.1; beta = 0.6; varepsilon = 0.00001; muvih = 0.05; sigma = 0.0001; theta1 = 0.3; theta2 = 1.0001; p = 0.6; muS = 0.01; muV = 0.01; muE = 0.02; omega = 0.0605; muL = 0.02; nu = 0.00375; muI = 0.12; tau = 0.1; A = 1; alpha = 0.2; S0 = 80; V0 = 50; Svih0 = 10; E0 = 8; L0 = 7; I0 = 5.

Here in Figure 4, the effective contact rate $\beta$ is higher than the precedent case, optimal control is more aggressive and appears earlier in the step 70, and is stopped before step 200; the consequence is that the trajectory of vaccinated

Figure 3. Optimal control induced by low effective contact rate.

Figure 4. Optimal control induced by high effective contact rate.

reached more later abscise axis and TB patients leave abscise axis later also, this means that spread of epidemic is retarded. we could conclude that more earlier system is controlling more later epidemic occurs in the population.

5.3.3. Optimal Control to a Absence of HIV Medical Coverage: $\sigma =1-T_{cd4}=1$ and $\beta =0.4$

Set of parameters values

pi1 = 50; pi2 = 1; pi3 = 0.1; beta = 0.4; varepsilon = 0.00001; muvih = 0.05; sigma = 1; theta1 = 0.3; theta2 = 2; p = 0.6; muS = 0.01; muV = 0.01; muE = 0.02; omega = 0.0605; muL = 0.02; nu = 0.00375; muI = 0.12; tau = 0.1; A = 1; alpha = 0.2; S0 = 80; V0 = 50; Svih0 = 10; E0 = 8; L0 = 7; I0 = 5.

In this case 5 $T_{cd4}$ is nul, that’s suppose HIV patients are not cured. In Figure 5, we can see that optimal control appears quickly at years fourth, in the same time hen TB patients are leaving the abscise axis and epidemic occurs. We remark also that optimal control disappears after few months. We can conclude that:the lack of HIV medical coverage accelerate the propagation of TB infectious, and when the velocity of epidemic is hight optimal control appears quickly and disappears quickly also. Also mass vaccination is useless when we have lack of HIV medical coverage.

5.3.4. Optimal Control Induced by a VIH Medical Coverage: $\sigma =1-T_{cd4}=0.01$

When VIH patients are cured, optimal control appears later after 40 years and stay more than five years Figure 6 here Optimal control is better than the precedent case 5 more than 200 TB patients at equilibrium. In Figure 6, we have above 150 TB patients at equilibrium.

5.3.5. Optimal Control Induced by a Hight VIH Medical Coverage: $\sigma =1-T_{cd4}=0.0001$

Set of parameters values

pi1 = 50; pi2 = 1; pi3 = 0.1; beta = 0.4; varepsilon = 0.00001; muvih = 0.05; sigma = 0.0001; theta1 = 0.3; theta2 = 1.0001; p = 0.6; muS = 0.01; muV = 0.01;

Figure 5. Optimal control induced by the lack of HIV medical coverage.

Figure 6. Optimal control induced by a VIH medical coverage.

muE = 0.02; omega = 0.0605; muL = 0.02; nu = 0.00375; muI = 0.12; tau = 0.1; A = 1; alpha = 0.2; S0 = 80; V0 = 50; Svih0 = 10; E0 = 8; L0 = 7; I0 = 5.

When HIV medical coverage is hight, control optimal appears more later step 80 and stay until step 200. In Figure 7, optimal control is more efficient because at endemic equilibrium we have 100 TB patients lesser than the precedent case where he have 150 TB patients in endemic equilibrium.

5.3.6. Eradication of Desease with Hight VIH Medical Coverage $\sigma =0.0001$,the Low VIH Prevalence Rate $\epsilon =0.00001$ and the Low VIH Vertical Transmission ${\pi }_3=0.1$

Set of parameters values

pi1 = 50; pi2 = 1; pi3 = 0.01; beta = 0.4; varepsilon = 0.00001; muvih = 0.05; sigma = 0.0001; theta1 = 0.3; theta2 = 1.0001; p = 0.6; muS = 0.01; muV = 0.01; muE = 0.02; omega = 0.0605; muL = 0.02; nu = 0.00375; muI = 0.12; tau = 0.1; A = 1; alpha = 0.2; S0 = 80; V0 = 50; Svih0 = 10; E0 = 8; L0 = 7; I0 = 5.

Here we have reproduction number $R_0$ less than unity, and the system 1 reach a disease free equilibrium (Figure 8). The number of TB patients decreases and reach abscise axis at step 50. Nowadays optimal control appears at step 70 and the vaccination rate coverage decreases progressively until 0.5% at step 100. The combination (hight VIH medical coverage $\sigma =0.0001$,the low VIH

Figure 7. Optimal control induced by a a hight VIH medical coverage.

Figure 8. Eradication of TB infection.

prevalence rate $\epsilon =0.00001$ and the low VIH vertical transmission ${\pi }_3=0.1$ ) assures the eradication of TB infectious in community. We can conclude that, well controlling the spread of VIH infectious (hight medical coverage of HIV patients, voluntary HIV screening campaign and HIV awareness campaigns) has positive effects in the propagation of TB infectious. Mass vaccination is not necessary when we have at least a good percentage of peoples who are vaccinated at birth.

All the scenarios are resumed in the following Figure 9.

5.3.7. Mass Vaccination Strategies Scenarios Synthesis

See Figure 9.

6. Discussions and Conclusion

The goals of this paper were to study the overall and asymptotic stability of the system around the point of equilibrium on one hand and, to use optimal control techniques to find mass vaccination strategies for each situation and assess impact on the second hand. We simulated the spread of tuberculosis/HIV coinfection and mass vaccination schedule. The database used was essentially made of the World Health organisation report 2017 and the data collected at the Centre Jamot Hospital where is housed the screening center of tuberculosis. The

Figure 9. Mass vaccination strategy scenarios.

approach of optimal control used in various cases with adapted tools, leads to some lessons learnt. The higher the effective contact rate will be, the earlier we should start the mass vaccination that could economically expensive (Figure 3 and Figure 4). The better the medical coverage will be, the the later we should engage mass vaccination as we reach the equilibriums (Figures 5-7). We eradicate the TB propagation by ensuring a better medical coverage, by reducing the HIV transmission from the mother to the child and by reducing the prevalence of HIV within the population (Figure 8).

In the area of HIV and TB co-infection, to reach the target of eradication of the TB propagation, we need to control the HIV propagation and make an emphasis on the immunization against the TB infection, the early screening and treatment of TB patients.

Acknowledgements

The first author acknowledges with thanks, the Chair of Computational Mathematics at Deusto Tech Laboratory in the University of Deusto, Bilbao (Basque Country, Spain), the Center for Tuberculosis Screening and Treatment Jamot Hospital, Yaounde, Cameroon.

Data Availability Statement

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

Funding

This work was supported by Women for Africa Foundation in collaboration with Diputacion Foral de Bizkaia and Deustotech-Deusto Foundation.

Conflicts of Interest

No potental conflict of interest was reported by the authors.

References

[1]	WHO, et al. (2016) Global Tuberculosis Report.
[2]	WHO. Tuberculosis Country pro_le (Cameroon).
[3]	Rubel, A.J. and Garro, L.C. (1992) Social and Cultural Factors in the Successful Control of Tuberculosis. Public Health Reports, 107, 626.
[4]	Narain, J.P., Raviglione, M.C. and Kochi, A. (1992) HIV-Associated Tuberculosis in Developing Countries: Epidemiology and Strategies for Prevention. Tubercle and Lung Disease, 73, 311-321. https://doi.org/10.1016/0962-8479(92)90033-G
[5]	Brewer, T.F., Heymann, S.J., Colditz, G.A., Wilson, M.E., Auerbach, K., Kane, D. and Fineberg, H.V. (1996) Evaluation of Tuberculosis Control Policies Using Computer Simulation. JAMA, 276, 1898-1903. https://doi.org/10.1001/jama.1996.03540230048034
[6]	Atun, R.A., Lebcir, R., Drobniewski, F. and Coker, R.J. (2005) Impact of an Effective Multidrugresistant Tuberculosis Control Programme in the Setting of an Immature HIV Epidemic: System Dynamics Simulation Model. International Journal of STD & AIDS, 16, 560-570. https://doi.org/10.1258/0956462054679124
[7]	West, R.W. and Thompson, J.R. (1997) Modeling the Impact of HIV on the Spread of Tuberculosis in the United States. Mathematical Biosciences, 143, 35-60. https://doi.org/10.1016/S0025-5564(97)00001-1
[8]	Lenhart, S. and Workman, J.T. (2007) Optimal Control Applied to Biological Models. CRC Press.
[9]	Zaman, G., Kang, Y.H. and Jung, I.H. (2008) Stability Analysis and Optimal Vaccination of an Sir Epidemic Model. BioSystems, 93, 240-249. https://doi.org/10.1016/j.biosystems.2008.05.004
[10]	Morton, R. and Wickwire, K.H. (1974) On the Optimal Control of a Deterministic Epidemic. Advances in Applied Probability, 6, 622-635. https://doi.org/10.2307/1426183
[11]	Kar, T. 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
[12]	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.
[13]	Trélat, E., Zhu, J. and Zuazua, E. (2018) Allee Optimal Control of a System in Ecology, Mathematical Models and Methods in Applied Sciences.
[14]	Ibanez, A. (2017) Optimal Control of the Lotka-Volterra System: Turnpike Property and Numerical Simulations. Journal of Biological Dynamics, 11, 25-41. https://doi.org/10.1080/17513758.2016.1226435
[15]	Agusto, F. and Adekunle, A. (2014) Optimal Control of a Two-Strain Tuberculosis-HIV/AIDS Co-Infection Model. Biosystems, 119, 20-44. https://doi.org/10.1016/j.biosystems.2014.03.006
[16]	Bowong, S. (2010) Optimal Control of the Transmission Dynamics of Tuberculosis. Nonlinear Dynamics, 61, 729-748. https://doi.org/10.1007/s11071-010-9683-9
[17]	Moualeu, D.P., Weiser, M., Ehrig, R. and Deuhard, P. (2015) Optimal Control for a Tuberculosis Model with Undetected Cases in Cameroon. Communications in Nonlinear Science and Numerical Simulation, 20, 986-1003. https://doi.org/10.1016/j.cnsns.2014.06.037
[18]	Nkamba, L., Ntaganda, J., Abboubakar, H., Kamgang, J. and Castelli, L. (2017) Global Stability of a Sveir Epidemic Model: Application to Poliomyelitis Transmission Dynamics. Open Journal of Modeling and Simulation, 5, 98-112. https://doi.org/10.4236/ojmsi.2017.51008
[19]	Hyman, J. and Li, J. (2005) Differential Susceptibility Epidemic Models. Journal of Mathematical Biology, 50, 626-644. https://doi.org/10.1007/s00285-004-0301-7
[20]	Nkamba, L.N., et al. (2019) Mathematical Model to Assess Vaccination and Effective Contact Rate Impact in the Spread of Tuberculosis. Journal of Biological Dynamics, 13, 26-42. https://doi.org/10.1080/17513758.2018.1563218
[21]	Diekmann, O. and Heesterbeek, J.A.P. (2000) Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation. Vol. 5, John Wiley & Sons, Hoboken.
[22]	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