Mathematical Modeling of Malaria Transmission Dynamics: Case of Burundi

Egide Ndamuzi¹, Paterne Gahungu²

¹Ecole Normale Supérieure de Bujumbura, Bujumbura, Burundi.
²University of Burundi, Bujumbura, Burundi.

DOI: 10.4236/jamp.2021.910156                

Abstract

Burundi, a country in East Africa with a temperate climate, has experienced in recent years a worrying growth of the Malaria epidemic. In this paper, a deterministic model of the transmission dynamics of malaria parasite in mosquito and human populations was formulated. The mathematical model was developed based on the SEIR model. An epidemiological threshold, R₀, called the basic reproduction number was calculated. The disease-free equilibrium point was locally asymptotically stable if R₀ < 1 and unstable if R₀ > 1. Using a Lyapunov function, we proved that this disease-free equilibrium point was globally asymptotically stable whenever the basic reproduction number is less than unity. The existence and uniqueness of endemic equilibrium were examined. With the Lyapunov function, we proved also that the endemic equilibrium is globally asymptotically stable if R₀ > 1. Finally, the system of equations was solved numerically according to Burundi’s data on malaria. The result from our model shows that, in order to reduce the spread of Malaria in Burundi, the number of mosquito bites on human per unit of time (σ), the vector population of mosquitoes (N_v), the probability of being infected for a human bitten by an infectious mosquito per unit of time (b) and the probability of being infected for a mosquito per unit of time (c) must be reduced by applying optimal control measures.

Keywords

Compartmental Model, Basic Reproduction Number, Local and Global Asymptotic Stability, Disease Free-Equilibrium, Endemic Equilibrium, Lyapunov Function

1. Introduction

Malaria is an infectious disease caused by a parasite of the genus plasmodium spread by the bites of certain species of female mosquitoes, the genus anopheles. Of the one hundred and twenty-three species of the genus plasmodium listed, only four are specifically human [1] : Plasmodium falciparum responsible for a large majority of deaths, and three others that cause “mild” forms of malaria that are not usually fatal: Plasmodium vivax, Plasmodium ovale, and Plasmodium malariae. Malaria is one of the most widespread vector-borne disease. Every year we find several hundred millions of people sick with malaria and between 1 and 3 million of deaths per year, the majority of them are children under 5 years and pregnant women. Eighty percent of all cases are recorded in sub-Saharan Africa. The disease reappears in areas where control efforts were effective and is emerging in areas considered free from the disease despite decades of eradication and control efforts of this epidemic [2] [3] [4] . In fact, it seems important to us, in this context, to make a scientific contribution in the quest for a good understanding of the transmission dynamics of the latter in order to optimize interventions for effective control and eradication of the disease. Globally, it is estimated that nearly forty percent of the population live in areas where malaria is endemic as reported in World Health Organization (WHO) fact sheet (2009). Burundi, one the East African region member, has experienced malaria as a major public health problem and is one of the main national health priority. According to data from the National Health Information System 2017, malaria is the leading cause of morbidity and mortality in Burundi with an incidence rate of 815‰. It constitutes 45.4% of the reasons for general consultations recorded in health facilities in 2017 with a rate of 50.5% in children under 5 years old. It is responsible for 50.8% of hospital deaths [5] . Among the plasmodial species responsible for malaria in humans, three exist in Burundi: Plasmodium falciparum, responsible for severe forms (81.6%), Plasmodium malariae (12.5%) and Plasmodium Ovale (5.8%) [5] .

Better knowledge of the disease, planning for the future and considering appropriate control measures are provided by mathematical models of the dynamics of malaria transmission. Mathematical modeling of malaria originated from the works of Ross [6] and MacDonald who did some modification to the Ross’s model and included superinfection [7] [8] . This has allowed the scientific community to refine these models until today. In order to model the disease, compartmental models of malaria and differential equations are constructed [3] [9] [10] .

The remainder of this paper is structured as follows: we present the formulation of the malaria transmission model in Section (2). In section (3), we develop the existence and positivity of solutions. The basic reproduction number is computed in section (4). In section (5), the stability analysis of disease equilibrium points is developed. Simulation and interpretation of the results are given in section (6). A conclusion is drawn in section (7).

2. Model Formulation

A mathematical model for the transmission dynamics of malaria is formulated. Using compartmental approach, our model is shown in Figure 1. The human population is divided into the SEIR compartmental model which consists of four classes: susceptible $S_h$, latent (exposed) $L_h$, infectious $I_h$ and recovered $R_h$. As the continuous exposure is necessary to ensure immunity [3] , we assume that the immunity is temporary. Those who have recovered have immunity against the disease for a certain period. Acquired immunity exists but the mechanisms are not yet fully understood [11] . A human is susceptible to give birth at ${\alpha }_h$ rate and loses immunity at $\gamma$ rate. The size of this class decreases as a result of natural mortality at ${\mu }_h$ rate as well as when an infectious mosquito bites a susceptible human with a probability b of being infected for a human bitten by a mosquito infectious per unit of time.

On the side of mosquito population, two compartments have been established: susceptible $S_v$ and infectious $I_v$. There is no recovered class for the vector as mosquitoes never recover [9] [10] or latent class (for simplicity of the model). A vector is susceptible to give birth at ${\alpha }_v$ rate. Susceptible mosquitoes that bite infectious human take gametocytes in the blood meal [7] and after fertilization, sporozoites are produced and migrated to the salivary glands ready to infect any susceptible host with probability c of being infected for a mosquito biting an infectious human per unit of time. The vector is then considered as infectious.

2.1. The Assumptions of the Model

1) All new borns are susceptible to the disease in these two populations.

2) Both susceptible, Latent (Exposed), Infectious and Recovered humans die at the same natural rate ${\mu }_h$ except for the class of infectious humans which has an additional death rate $\eta$ caused by the disease.

3) Both susceptible and Infectious mosquitoes die at the same natural rate ${\mu }_v$.

4) Neither human nor mosquito total population is constant.

5) Recovered humans are still able to transmit the disease but at a lower rate, therefore negligible in our case.

6) Human recover from infection (without immunity) and move right back to

Figure 1. Flow diagram for the transmission dynamics of malaria.

susceptible state or gain temporary immunity before losing it and returning to the susceptible class.

The differential equations describing the transmission dynamics of disease are given in the following section according to the compartmental model in Figure 1 and the Balance Law.

2.2. Establishment of the Model Equations

The equations modeling this dynamic are:

$\{\begin{array}{l}\frac{\text{d}{S}_h}{\text{d}t}={\alpha }_h{N}_h-\frac{\sigma b{S}_h{I}_v}{N_h}+\gamma R_h-{\mu }_h{S}_h,\\ \frac{\text{d}{L}_h}{\text{d}t}=\frac{\sigma b{S}_h{I}_v}{N_h}-\left({\mu }_h+\delta \right)L_h,\\ \frac{\text{d}{I}_h}{\text{d}t}=\delta L_h-\left(\eta +{\mu }_h+\lambda \right)I_h,\\ \frac{\text{d}{R}_h}{\text{d}t}=\lambda I_h-\left({\mu }_h+\gamma \right)R_h,\\ \frac{\text{d}{S}_v}{\text{d}t}={\alpha }_v{N}_v-\frac{\sigma c{S}_v{I}_h}{N_h}-{\mu }_v{S}_v,\\ \frac{\text{d}{I}_v}{\text{d}t}=\frac{\sigma c{S}_v{I}_h}{N_h}-{\mu }_v{I}_v.\end{array}$ (1)

2.3. Description of Parameters and Variables

1) The considered parameters are:

${\alpha }_h$ : birth rate for humans;

${\alpha }_v$ : birth rate for mosquitoes;

$\sigma$ : number of bites of a mosquito to a human per unit of time;

b: the probability of being infected for a human bitten by an infectious mosquito per unit of time;

c: the probability of being infected for a mosquito biting an infectious human per unit of time;

$\delta$ : rate of progression from latent to infectious state for humans;

$\lambda$ : rate of progression from infectious to recovered state for humans;

$\gamma$ : rate of progression from recovered to susceptible state for humans;

${\mu }_h$ : human (natural) mortality rate;

$\eta$ : human malaria-induced death rate;

${\mu }_v$ : mosquito mortality rate.

2) The considered variables are:

$N_h$ : total population for humans (host);

$N_v$ : total population for mosquitoes (vector);

$S_h$ : class (number) of susceptible humans;

$L_h$ : class (number) of latent humans;

$I_h$ : class (number) of infectious humans;

$R_h$ : class (number) of recovered humans;

$S_v$ : class (number) of susceptible mosquitoes;

$I_v$ : class (number) of infectious mosquitoes.

Remark. All parameters of the model are assumed to be non-negative.

3. Existence and Positivity of Solutions

We show that the system is epidemiologically and mathematically well-posed in the feasible region $\Omega$ given by

$\Omega =\left\{\left(S_h\mathrm{,}{L}_h\mathrm{,}{I}_h\mathrm{,}{R}_h\mathrm{,}{S}_v\mathrm{,}{I}_v\right)\in {ℝ}_{+}^6\mathrm{<mo>\; :\; </mo>0}\le {\mu }_h\le {\alpha }_h\mathrm{,\; <mi>\; </mi>0}\le {\mu }_v\le {\alpha }_v\right\}$ (2)

Theorem 1. The system (1) always admits positive solutions for all positive initial conditions and the biological region $\Omega \in {ℝ}_{+}^6$ is positively invariant and globally attractive for the same system.

Proof. The system (1) can be written in the form $\stackrel{\dot{}}{x}\left(t\right)=f\left(x\left(t\right)\right)$ where

$x=\left(x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4\mathrm{,}{x}_5\mathrm{,}{x}_6\right)=\left(S_h\mathrm{,}{L}_h\mathrm{,}{I}_h\mathrm{,}{R}_h\mathrm{,}{S}_v\mathrm{,}{I}_v\right)$ (3)

and

$f\left(x\right)=\left(f_1\left(x\right)\mathrm{,}{f}_2\left(x\right)\mathrm{,}{f}_3\left(x\right)\mathrm{,}{f}_4\left(x\right)\mathrm{,}{f}_5\left(x\right)\mathrm{,}{f}_6\left(x\right)\right)$ (4)

Then, the system (1) becomes

$\{\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}={\stackrel{\dot{}}{x}}_1={\alpha }_h\left(x_1+x_2+x_3+x_4\right)-\frac{\sigma b{x}_1{x}_6}{x_1+x_2+x_3+x_4}+\gamma x_4-{\mu }_h{x}_1,\\ \frac{\text{d}{x}_2}{\text{d}t}={\stackrel{\dot{}}{x}}_2=\frac{\sigma b{x}_1{x}_6}{x_1+x_2+x_3+x_4}-A{x}_2,\\ \frac{\text{d}{x}_3}{\text{d}t}={\stackrel{\dot{}}{x}}_3=\delta x_2-B{x}_3,\\ \frac{\text{d}{x}_4}{\text{d}t}={\stackrel{\dot{}}{x}}_4=\lambda x_3-C{x}_4,\\ \frac{\text{d}{x}_5}{\text{d}t}={\stackrel{\dot{}}{x}}_5={\alpha }_v\left(x_5+x_6\right)-\frac{\sigma c{x}_3{x}_5}{x_1+x_2+x_3+x_4}-{\mu }_v{x}_5,\\ \frac{\text{d}{x}_6}{\text{d}t}={\stackrel{\dot{}}{x}}_6=\frac{\sigma c{x}_3{x}_5}{x_1+x_2+x_3+x_4}-{\mu }_v{x}_6.\end{array}$ (5)

where $A={\mu }_h+\delta ,\mathrm{<mi>\; </mi>}B=\eta +{\mu }_h+\lambda$ and $C={\mu }_h+\gamma$.

As the second member of the system (1) is locally Lipschitzian, then it admits a unique solution. It is clear that $\forall j=\mathrm{1,}\cdots \mathrm{,6}$, $f_j\left(x\right)\ge 0$ if $x\in {\left[\mathrm{0,}+\infty \right]}^6$, hence the positivity of the solution.

Added member to member, the first four equations of the system (1) and the last two equations of the same system, we notice that the total population of human and vector satisfies the following equations:

$\{\begin{array}{l}{\stackrel{\dot{}}{N}}_h={\alpha }_h{N}_h-{\mu }_h{N}_h-\eta I_h,\\ {\stackrel{\dot{}}{N}}_v={\alpha }_v{N}_v-{\mu }_v{N}_v\end{array}$ (6)

From the Equation (6), we have ${\stackrel{\dot{}}{N}}_h\le {\alpha }_h{N}_h-{\mu }_h{N}_h\Rightarrow {\stackrel{\dot{}}{N}}_h\le 0$ if ${\alpha }_h\le {\mu }_h$. So $N_h\left(t\right)\le N_h\left(0\right)\exp \left[\left({\alpha }_h-{\mu }_h\right)t\right]$. In the same way, we have ${\stackrel{\dot{}}{N}}_v\le 0$ if ${\alpha }_v\le {\mu }_v$ and $N_v\left(t\right)\le N_v\left(0\right)\exp \left[\left({\alpha }_v-{\mu }_v\right)t\right]$. We conclude that the region $\Omega$ is positively invariant.

Moreover, if ${\alpha }_h<{\mu }_h$ and ${\alpha }_v<{\mu }_v$, we see that the solution of the system (1) enter in $\Omega$ in finite time, i.e. $N_h\left(t\right)$ asymptotically approaches ${\alpha }_h$ and $N_v\left(t\right)$ approaches ${\alpha }_v$ asymptotically.

Therefore, all solutions in ${ℝ}_{+}^6$ eventually enter $\Omega$ i.e. the biological region $\Omega$ is globally attractive for the same system. So the problem is then mathematically and epidemiologically well posed. Hence, every solution of the model (1) with initial conditions in $\Omega$ remains in $\Omega$ for all $t>0$. o

4. Basic Reproduction Number

In order to assess the stability of disease equilibria, the computation of the basic reproduction number $R_0$ is needed. The next generation matrix method for calculating $R_0$ as established by Watmough and Van Den Driessche [12] will be used for our case. We know that

$N_h=S_h+L_h+I_h+R_h\text{with}\{\begin{array}{l}{L}_h\mathrm{,}{I}_h\mathrm{<mo>\; :\; </mo>}\text{infected states}\mathrm{,}\\ S_h\mathrm{,}{R}_h\mathrm{<mo>\; :\; </mo>}\text{uninfected states}\mathrm{.}\end{array}$ (7)

and

$N_v=S_v+I_v\text{with}\{\begin{array}{l}{I}_v\mathrm{<mo>\; :\; </mo>}\text{infected states}\mathrm{,}\\ S_v\mathrm{<mo>\; :\; </mo>}\text{uninfected states}\mathrm{.}\end{array}$ (8)

In the infection free state, we have $L_h=I_h=R_h=0\Rightarrow N_h=S_h$ and $I_v=0\Rightarrow N_v=S_v$. So for $\left(L_h\mathrm{,}{I}_h\mathrm{,}{I}_v\right)$ and using the expressions for $N_h$ and $N_v$ that we have just found above, we have the following linear system:

$\{\begin{array}{l}\frac{\text{d}{L}_h}{\text{d}t}=\sigma b{I}_v-\left({\mu }_h+\delta \right)L_h,\\ \frac{\text{d}{I}_h}{\text{d}t}=\delta L_h-\left(\eta +{\mu }_h+\lambda \right)I_h,\\ \frac{\text{d}{I}_v}{\text{d}t}=\frac{\sigma c{N}_v{I}_h}{N_h}-{\mu }_v{I}_v.\end{array}$ (9)

Let $X={\left(L_h\mathrm{,}{I}_h\mathrm{,}{I}_v\right)}^{\prime }$ where the prime denotes the transpose. The infection subsystem linearized is then in the form:

$\stackrel{\dot{}}{X}=\left(T+\Sigma \right)X$ (10)

where T is the transmission matrix and $\Sigma$ is the transition matrix such that:

$T=\left(\begin{array}{ccc}0& 0& \sigma b\\ 0& 0& 0\\ 0& \frac{\sigma c{N}_v}{N_h}& 0\end{array}\right)\text{and}\Sigma =\left(\begin{array}{ccc}-\left({\mu }_h+\delta \right)& 0& 0\\ \delta & -\left(\eta +{\mu }_h+\lambda \right)& 0\\ 0& 0& -{\mu }_v\end{array}\right)$ (11)

So, the next generation matrix $K=-T{\Sigma }^{-1}$ where ${\Sigma }^{-1}=\frac{1}{\left|\Sigma \right|}\cdot {\Sigma }^{Adj.}$, ${\Sigma }^{Adj\mathrm{.}}={\left[Com\left(\Sigma \right)\right]}^{\prime }$ and $Com\left(\Sigma \right)$ is the matrix of cofactors of $\Sigma$. Thus,

$K=\left(\begin{array}{ccc}0& 0& \frac{\sigma b}{{\mu }_v}\\ 0& 0& 0\\ \frac{\delta \sigma c{N}_v}{N_h\left({\mu }_h+\delta \right)\left(\eta +{\mu }_h+\lambda \right)}& \frac{\sigma c{N}_v}{N_h\left(\eta +{\mu }_h+\lambda \right)}& 0\end{array}\right)$ (12)

By computation the dominant eigenvalue of K which is the reproduction number basic, we have

$R_0=\rho \left(K\right)=\sigma \mathrm{<mi>\; </mi>}\sqrt{\frac{\delta }{{\mu }_h+\delta }\frac{bc{N}_v}{{\mu }_v{N}_h\left(\eta +{\mu }_h+\lambda \right)}}$ (13)

As $det\left(K\right)=0$, we notice that the K found is the next generation matrix with large domain denoted $K_L$. So, from this $K_L$, we deduce $K^{\mathrm{<mo>\; *\; </mo>}}$ which is the next generation matrix (i.e. $K^{\mathrm{<mo>\; *\; </mo>}}$ is the reduction of $K_L$ ), knowing that $\rho \left(K_L\right)=\rho \left(K^{\mathrm{<mo>\; *\; </mo>}}\right)$. We have

$K^{\mathrm{<mo>\; *\; </mo>}}=E^{\prime }{K}_{L}E=-E^{\prime }T\Sigma E\mathrm{<mi>\; </mi>}\text{with}\mathrm{<mi>\; </mi>}E=\left(\begin{array}{cc}1& 0\\ 0& 0\\ 0& 1\end{array}\right)$

Thus,

$K^{\mathrm{<mo>\; *\; </mo>}}=\left(\begin{array}{cc}0& \frac{\sigma b}{{\mu }_v}\\ \frac{\delta \sigma c{N}_v}{N_h\left({\mu }_h+\delta \right)\left(\eta +{\mu }_h+\lambda \right)}& 0\end{array}\right)$ (14)

Since $K^{\mathrm{<mo>\; *\; </mo>}}$ is a $2\times 2$ matrix, we have

$\begin{array}{c}{R}_0=\rho \left(K^{\mathrm{<mo>\; *\; </mo>}}\right)=\frac{1}{2}\left[tr\left(K^{\mathrm{<mo>\; *\; </mo>}}\right)+\sqrt{tr{\left(K^{\mathrm{<mo>\; *\; </mo>}}\right)}^2-4det\left(K^{\mathrm{<mo>\; *\; </mo>}}\right)}\right]\\ =\sigma \mathrm{<mi>\; </mi>}\sqrt{\frac{\delta }{{\mu }_h+\delta }\frac{bc{N}_v}{{\mu }_v{N}_h\left(\eta +{\mu }_h+\lambda \right)}}\end{array}$ (15)

5. Stability Analysis of Disease Equilibria

There are two types of disease equilibria, which are obtained by solving the system (1), for ${\stackrel{\dot{}}{S}}_h={\stackrel{\dot{}}{L}}_h={\stackrel{\dot{}}{I}}_h={\stackrel{\dot{}}{R}}_h={\stackrel{\dot{}}{S}}_v={\stackrel{\dot{}}{I}}_v=0$. In particular the disease-free equilibrium point (DFE), $E_0$, which is obtained by adding the condition $L_h=I_h=R_h=I_v=0$. With these conditions, the resolution of the system (1) then gives

$E_0=\left(\frac{{\alpha }_h{N}_h}{{\mu }_h}\mathrm{<mo>\; ;\; </mo>0\; <mo>\; ;\; </mo>0\; <mo>\; ;\; </mo>0\; <mo>\; ;\; </mo>}\frac{{\alpha }_v{N}_v}{{\mu }_v}\mathrm{<mo>\; ;\; </mo>0}\right)$ (16)

And the endemic equilibrium “EE” point, $E_1$, is obtained with the condition $S_h\ne \mathrm{0\; <mo>\; ;\; </mo>}{L}_h\ne \mathrm{0\; <mo>\; ;\; </mo>}{I}_h\ne \mathrm{0\; <mo>\; ;\; </mo>}{R}_h\ne \mathrm{0\; <mo>\; ;\; </mo>}{S}_v\ne \mathrm{0\; <mo>\; ;\; </mo>}{I}_v\ne 0$. Implicitily, we have

$E_1=\left(S_h^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}{L}_h^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}{I}_h^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}{R}_h^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}{S}_v^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}{I}_v^{\mathrm{<mo>\; *\; </mo>}}\right)$ (17)

5.1. Stability of DFE

Theorem 2. The disease-free equilibrium, $E_0$, (DFE) of Equation (1), given by Equation (16), is locally asymptotically stable if $R_0<1$, and unstable if $R_0>1$.

Proof. Let $J\left(M\right)$ be the Jacobian of Equation (1) at point $M=\left(S_h\mathrm{,}{L}_h\mathrm{,}{I}_h\mathrm{,}{R}_h\mathrm{,}{S}_v\mathrm{,}{I}_v\right)$

$J\left(M\right)=\left(\begin{array}{cccccc}-\frac{\sigma b{I}_v}{N_h}-{\mu }_h& 0& 0& \gamma & 0& -\frac{\sigma b{S}_h}{N_h}\\ \frac{\sigma b{I}_v}{N_h}& -\left({\mu }_h+\delta \right)& 0& 0& 0& \frac{\sigma b{S}_h}{N_h}\\ 0& \delta & -\left(\eta +{\mu }_h+\lambda \right)& 0& 0& 0\\ 0& 0& \lambda & -\left({\mu }_h+\gamma \right)& 0& 0\\ 0& 0& -\frac{\sigma c{S}_v}{N_h}& 0& -\frac{\sigma c{I}_h}{N_h}-{\mu }_v& 0\\ 0& 0& \frac{\sigma c{S}_v}{N_h}& 0& \frac{\sigma c{I}_h}{N_h}& -{\mu }_v\end{array}\right)$ (18)

The Jacobian matrix in Equation (1) at point $E_0$ then becomes:

$J\left(E_0\right)=\left(\begin{array}{cccccc}-{\mu }_h& 0& 0& \gamma & 0& -\frac{\sigma b{\alpha }_h}{{\mu }_h}\\ 0& -\left({\mu }_h+\delta \right)& 0& 0& 0& \frac{\sigma b{\alpha }_h}{{\mu }_h}\\ 0& \delta & -\left(\eta +{\mu }_h+\lambda \right)& 0& 0& 0\\ 0& 0& \lambda & -\left({\mu }_h+\gamma \right)& 0& 0\\ 0& 0& -\frac{\sigma c{\alpha }_v{N}_v}{{\mu }_v{N}_h}& 0& -{\mu }_v& 0\\ 0& 0& \frac{\sigma c{\alpha }_v{N}_v}{{\mu }_v{N}_h}& 0& 0& -{\mu }_v\end{array}\right)$ (19)

Consider $P\left({\lambda }^{\mathrm{<mo>\; *\; </mo>}}\right)$ the characteristic polynomial of $J\left(E_0\right)$, we then have:

$\left|J\left(E_0\right)-{\lambda }^{\mathrm{<mo>\; *\; </mo>}}{\mathbb{I}}_6\right|=P(\; \lambda \; *\; )$

$P\left({\lambda }^{\mathrm{<mo>\; *\; </mo>}}\right)=\left|\begin{array}{cccccc}-{\mu }_h-{\lambda }^{\mathrm{<mo>\; *\; </mo>}}& 0& 0& \gamma & 0& -\frac{\sigma b{\alpha }_h}{{\mu }_h}\\ 0& -\left({\mu }_h+\delta \right)-{\lambda }^{\mathrm{<mo>\; *\; </mo>}}& 0& 0& 0& \frac{\sigma b{\alpha }_h}{{\mu }_h}\\ 0& \delta & -\left(\eta +{\mu }_h+\lambda \right)-{\lambda }^{\mathrm{<mo>\; *\; </mo>}}& 0& 0& 0\\ 0& 0& \lambda & -\left({\mu }_h+\gamma \right)-{\lambda }^{\mathrm{<mo>\; *\; </mo>}}& 0& 0\\ 0& 0& -\frac{\sigma c{\alpha }_v{N}_v}{{\mu }_v{N}_h}& 0& -{\mu }_v-{\lambda }^{\mathrm{<mo>\; *\; </mo>}}& 0\\ 0& 0& \frac{\sigma c{\alpha }_v{N}_v}{{\mu }_v{N}_h}& 0& 0& -{\mu }_v-{\lambda }^{\mathrm{<mo>\; *\; </mo>}}\end{array}\right|$ (20)

By expanding, we have:

$\begin{array}{c}P\left({\lambda }^{\ast }\right)={\lambda }^{\ast }{}^6+\left(A+B+C+{\mu }_h+2{\mu }_v\right){\lambda }^{\ast }{}^5\\ +[{\mu }_v\left(2{\mu }_h+{\mu }_v\right)+\left({\mu }_h+2{\mu }_v\right)\left(A+B+C\right)+AB+AC+BC]{\lambda }^{\ast }{}^4\\ +[{\mu }_h{\mu }_v^2+{\mu }_v\left(2{\mu }_h+{\mu }_v\right)\left(A+B+C\right)\\ +\left({\mu }_h+2{\mu }_v\right)\left(AB+AC+BC\right)-D]{\lambda }^{\ast }{}^3\end{array}$

$\begin{array}{l}+[{\mu }_h{\mu }_v^2\left(A+B+C\right)+{\mu }_v\left(2{\mu }_h+{\mu }_v\right)\left(AB+AC+BC\right)\\ +\left({\mu }_h+2{\mu }_v\right)ABC-\left(C+{\mu }_h+{\mu }_v\right)D]{\lambda }^{\ast }{}^2\\ +[{\mu }_h{\mu }_v^2\left(AB+AC+BC\right)+{\mu }_v\left(2{\mu }_h+{\mu }_v\right)ABC\\ -\left({\mu }_h{\mu }_v+C{\mu }_h+C{\mu }_v\right)D]{\lambda }^{\ast }+ABC{\mu }_h{\mu }_v^2-CD{\mu }_h{\mu }_v\end{array}$ (21)

where $A={\mu }_h+\delta \mathrm{<mi>\; </mi>},\mathrm{<mi>\; </mi>}B={\mu }_h+\eta +\lambda ,\mathrm{<mi>\; </mi>}C={\mu }_h+\gamma$ and $D=\frac{bc{\sigma }^2\delta {\alpha }_h{\alpha }_v{N}_v}{{\mu }_h{\mu }_v{N}_h}$.

With Descartes’ rule of signs [13] , if $R_0<1$, all the coefficients of $P\left({\lambda }^{\mathrm{<mo>\; *\; </mo>}}\right)$ are strictly positive then $P\left({\lambda }^{\mathrm{<mo>\; *\; </mo>}}\right)$ does not have a positive root. So, $J\left(E_0\right)$ has all its eigenvalues with strictly negative real part, hence $E_0$ is locally asymptotically stable by the Poincarré-Lyapunov theorem [13] . Otherwise, with the same rule, if $R_0>1$, there is a variation of coefficients of $P\left({\lambda }^{\mathrm{<mo>\; *\; </mo>}}\right)$, then $P\left({\lambda }^{\mathrm{<mo>\; *\; </mo>}}\right)$ has at least one positive root. So, $J\left(E_0\right)$ has at least one eigenvalue with a strictly positive real part, hence then $E_0$ is unstable by the Poincarré-Lyapunov theorem. o

Theorem 3. The disease-free equilibrium point (DFE) of the model, denoted $E_0$, is globally asymptotically stable if $R_0<1$ and unstable if $R_0>1$.

Proof. Consider then the following Lyapunov candidate function for $E_0$ :

$V=\frac{1}{2}{\left(S_h-\frac{{\alpha }_h{N}_h}{{\mu }_h}\right)}^2+\frac{1}{2}{L}_h^2+\frac{1}{2}{I}_h^2+\frac{1}{2}{R}_h^2+\frac{1}{2}{\left(S_v-\frac{{\alpha }_v{N}_v}{{\mu }_v}\right)}^2+\frac{1}{2}{I}_v^2$ (22)

Its first derivative gives us

$\begin{array}{c}\stackrel{\dot{}}{V}=2{\alpha }_h{N}_h{S}_h+\gamma S_h{R}_h+\frac{\sigma b{\alpha }_h{S}_h{I}_v}{{\mu }_h}+\frac{\sigma b{S}_h{L}_h{I}_v}{N_h}+\delta L_h{I}_h+\lambda I_h{R}_h\\ +2{\alpha }_v{N}_v{S}_v+\frac{\sigma c{\alpha }_v{N}_v{S}_v{I}_h}{{\mu }_v{N}_h}+\frac{\sigma c{S}_v{I}_v{I}_h}{N_h}-(\frac{\sigma b{S}_h^2{I}_v}{N_h}+{\mu }_h{S}_h^2+\frac{{\alpha }_h^2{N}_h^2}{{\mu }_h}\\ +\frac{\gamma {\alpha }_h{N}_h{R}_h}{{\mu }_h})-[\left({\mu }_h+\delta \right)L_h^2+\left(\eta +\lambda +{\mu }_h\right)I_h^2+\left(\gamma +{\mu }_h\right)R_h^2+\frac{\sigma c{S}_v^2{I}_h}{N_h}\\ +{\mu }_v{S}_v^2+\frac{{\alpha }_v^2{N}_v^2}{{\mu }_v}+{\mu }_v{I}_v^2]\end{array}$ (23)

We notice that if $R_0<1$, the Lyapunov candidate function V thus defined is in the strict sense in the whole domain of $\Omega$. Hence, the largest compact invariant set in $\left\{\left(S_h,L_h,I_h,R_h,S_v,I_v\right)\in \Omega :\stackrel{\dot{}}{V}<0\right\}$ is the DFE, $E_0$. Therefore, by Lyapunov-Lasalle’s Invariance Principle, the disease free equilibrium is globally asymptotically stable in $\Omega$ if $R_0<1$ and unstable if $R_0>1$. o

5.2. Existence and Uniqueness of EE

Proposition 4. Endemic Equilibrium “EE” exists if and only if $R_0>1$.

Proof. Recall that the endemic equilibrium point, $E_1$, found previously was determined implicitly due to the complexity of this system. With the notions of the equilibrium point, the system (1) becomes:

$\{\begin{array}{l}0={\alpha }_h{N}_h-\frac{\sigma b{S}_h^{*}{I}_v^{*}}{N_h}+\gamma R_h^{*}-{\mu }_h{S}_h^{*},\\ 0=\frac{\sigma b{S}_h^{*}{I}_v^{*}}{N_h}-\left({\mu }_h+\delta \right)L_h^{*},\\ 0=\delta L_h^{*}-\left(\eta +{\mu }_h+\lambda \right)I_h^{*},\\ 0=\lambda I_h^{*}-\left({\mu }_h+\gamma \right)R_h^{*},\\ 0={\alpha }_v{N}_v-\frac{\sigma c{S}_v^{*}{I}_h^{*}}{N_h}-{\mu }_v{S}_v^{*},\\ 0=\frac{\sigma c{S}_v^{*}{I}_h^{*}}{N_h}-{\mu }_v{I}_v^{*}.\end{array}$ (24)

From this system, we deduce that

$\{\begin{array}{l}{S}_h^{*}=\frac{N_h\left({\alpha }_h{N}_h+\gamma R_h^{*}\right)}{\sigma b{I}_v^{*}+{\mu }_h{N}_h},\\ L_h^{*}=\frac{\sigma b{S}_h^{*}{I}_v^{*}}{N_h\left({\mu }_h+\delta \right)},\\ I_h^{*}=\frac{\delta }{\eta +{\mu }_h+\lambda }{L}_h^{*},\\ R_h^{*}=\frac{\lambda }{{\mu }_h+\gamma }{I}_h^{*},\\ S_v^{*}=\frac{{\alpha }_v{N}_v{N}_h}{\sigma c{I}_h^{*}+{\mu }_v{N}_h},\\ I_v^{*}=\frac{\sigma c{S}_v^{*}{I}_h^{*}}{N_h{\mu }_v}\end{array}$ (25)

By making substitutions, we notice that the EE satisfies the following polynomial function:

$P\left(I_h^{*}\right)=I_h^{*}\left({\psi }_1{I}_h^{*}+{\psi }_2\right)=0$ (26)

where

$\begin{array}{l}{\psi }_1={\sigma }^{2}bc\delta \gamma \lambda {\alpha }_v{N}_v\mathrm{<mi>\; </mi>}-ABC\left({\sigma }^{2}bc{\alpha }_v{N}_v+\sigma c{\mu }_h{\mu }_v{N}_h\right)\mathrm{,}\\ {\psi }_2=N_{h}C\left({\sigma }^{2}bc\delta {\alpha }_v{\alpha }_h{N}_v\mathrm{<mi>\; </mi>}-AB{\mu }_h{\mu }_v^2{N}_h\right)\mathrm{.}\end{array}$ (27)

From Equation (26), we have: $I_h^{*}=0$ which is the solution of DFE, or ${\psi }_1{I}_h^{*}+{\psi }_2=0$, hence the uniqueness of endemic equilibrium whenever $R_0>1$, otherwise, no endemic equilibrium. o

5.3. Stability of EE

Theorem 5. The endemic equilibrium point (EE), noted $E_1$, is globally asymptotically stable if $R_0>1$.

Proof. The proof of this theorem is based on the Lyapunov-Lasalle’s Invariance Principle. Consider then the following Lyapunov candidate function for $E_1$ :

$\begin{array}{c}V=\frac{1}{2}{\left(S_h-S_h^{*}\right)}^2+\frac{1}{2}{\left(L_h-L_h^{*}\right)}^2+\frac{1}{2}{\left(I_h-I_h^{*}\right)}^2+\frac{1}{2}{\left(R_h-R_h^{*}\right)}^2\\ +\frac{1}{2}{\left(S_v-S_v^{*}\right)}^2+\frac{1}{2}{\left(I_v-I_v^{*}\right)}^2\end{array}$ (28)

By deriving this function once, we have:

$\begin{array}{c}\stackrel{\dot{}}{V}=\left(S_h-S_h^{*}\right){\stackrel{\dot{}}{S}}_h+\left(L_h-L_h^{*}\right){\stackrel{\dot{}}{L}}_h+\left(I_h-I_h^{*}\right){\stackrel{\dot{}}{I}}_h\\ +\left(R_h-R_h^{*}\right){\stackrel{\dot{}}{R}}_h+\left(S_v-S_v^{*}\right){\stackrel{\dot{}}{S}}_v+\left(I_v-I_v^{*}\right){\stackrel{\dot{}}{I}}_v\end{array}$

$\begin{array}{l}=[S_h\left({\alpha }_h{N}_h+\gamma R_h\right)+S_h^{*}\left(\frac{\sigma b{S}_h{I}_v}{N_h}+{\mu }_h{I}_h\right)+\frac{\sigma b{S}_h{I}_v{L}_h}{N_h}+\left({\mu }_h+\delta \right)L_h{L}_h^{*}\\ +\delta L_h{I}_h+\left(\eta +{\mu }_h+\lambda \right)I_h{I}_h^{*}+\lambda I_h{R}_h+\left({\mu }_h+\gamma \right)R_h{R}_h^{*}+{\alpha }_v{N}_v{S}_v+\frac{\sigma c{S}_v{I}_h{S}_v^{*}}{N_h}\\ +{\mu }_v{S}_v{S}_v^{*}+\frac{\sigma c{S}_v{I}_h{I}_v}{N_h}+{\mu }_v{I}_v{I}_v^{*}]-[S_h\left(\frac{\sigma b{S}_h{I}_v}{N_h}+{\mu }_h{s}_h\right)+S_h^{*}\left({\alpha }_h{N}_h+\gamma R_h\right)\\ +\left({\mu }_h+\delta \right)L_h^2+\frac{\sigma b{S}_h{I}_v{L}_h^{*}}{N_h}+\left(\eta +{\mu }_h+\lambda \right)I_h^2+\delta L_h{I}_h^{*}+\left({\mu }_h+\gamma \right)R_h^2+\lambda I_h{R}_h^{*}\\ +\frac{\sigma c{S}_v^2{I}_h}{N_h}+{\mu }_v\left(S_v^2+I_v^2\right)+{\alpha }_v{N}_v{S}_v^{*}+\frac{\sigma c{S}_v{I}_h{I}_v^{*}}{N_h}]\end{array}$ (29)

Consequently, the largest compact invariant set in $\left\{\left(S_h\mathrm{,}{L}_h\mathrm{,}{I}_h\mathrm{,}{R}_h\mathrm{,}{S}_v\mathrm{,}{I}_v\right)\in \Omega \mathrm{<mo>\; :\; </mo>}\stackrel{\dot{}}{V}\le 0\right\}$ is the EE, $E_1$. Therefore, by Lyapunov-Lasalle’s Invariance Principle, it follows that every solution in $\Omega$ approaches $E_1$ for $R_0>1$ as $t\to \infty$. So the EE, $E_1$, is globally asymptotically stable. The epidemiological implication of this result is that malaria will persist and invade population whenever $R_0>1$. o

6. Simulation and Interpretation

In this section, we will analyze the model (1) numerically. In particular, we will be interested in the numerical analysis of the infected human population with malaria. To do this, we collect a reasonable set of data for the parameters of the model. Most of these are obtained from models developed in the literature and other parameters will be assumed.

In what follows, we will use the below parameters:

${\alpha }_h$ : The human birth rate is 38.377 per year per 1000 inhabitants [14] .

${\alpha }_v$ : The vector birth rate is 130 adult female mosquitoes per year per 1000 female mosquitoes [15] .

$\sigma$ : 20.425 bites per man per might (average bites per man per day inside and outside the home) [Assumed].

b: The probability of being infected for a human bitten by an infectious mosquito per unit of time is 0.0181 [5] .

c: The probability of being infected for a mosquito biting an infectious human per unit of time is 0.11 per day [16] .

$\delta$ : Rate of progression from latent to infectious state for humans is 0.1 [Assumed].

$\lambda$ : Rate of progression from infectious to recovered state for humans is 0.0022 per day [17] .

$\gamma$ : Rate of progression from recovered to susceptible state for humans 5.5 × 10⁻⁴ [17] .

${\mu }_h$ : Human (natural) mortality rate is 7.766 per year per 1000 inhabitants [14] .

$\eta$ : Human malaria-induced death rate is 32% [18] .

${\mu }_v$ : Mosquito mortality rate is 0.003 [15] .

The data of confirmed cases used, are obtained from the Ministry of Public Health and the Fight against HIV/AIDS for the period 2012 to 2018 [19] . As given in the literature, we initialize parameters and use minimize function with noise in python for obtaining parameter values corresponding to the data in study. Initial conditions of states are $S_h=498901$ ; $L_h=374176$ ; $I_h=823$ ; $R_h=560$ ; $S_v=3582$ ; $I_v=66$ for the following Figure 2.

Using the parameters derived from the fitted model above and taking $N_v=1000$ [Assumed]; $N_h=11890781$ [13] , the basic reproduction number as mentioned in section (4) is 0.8092009797722537.

According to this value and to the section (5), we conclude that the average of an individual infected produces less than one new infected individual during their infectious period, hence the disease-free equilibrium, $E_0$, (DFE) of Equation (1), given by Equation (16), is locally asymptotically stable.

When we assume $N_v=10000$ [assumed], using the same parameters and let fix $N_h=11890781$ [14] then, the value of basic reproduction is 2.5589181809201627. That implies that each infected individual produces an average of more than one; the infection and disease persist and invade the human population. It implies that the infection is greater than recovery rate. Therefore the strategy is to

Figure 2. Infected human population.

take all possible measures to avoid/reduce transmission and improve recovery.

7. Conclusion

In this paper, a malaria model was developed and analyzed to study the stability of both disease-free and endemic equilibrium point. Using matrix generation approach, the basic reproduction number $R_0$ was computed. Therefore, the disease-free equilibrium of the model obtained is both locally and globally stable for $R_0<1$. It is also shown that the endemic equilibrium solution of the model is globally asymptotically stable if $R_0>1$. The results from our model show that, in order to reduce the spread of disease, the number of mosquito bites to a human per unit of time ( $\sigma$ ), the vector population of mosquitoes ( $N_v$ ), the probability of being infected for a human bitten by an infectious mosquito by unit of time (b) and the probability of being infected for a mosquito biting a human infectious per unit of time (c) must be decreased significantly. Hence there is a need to use mosquito nets and insecticides or other strategies that can reduce these parameters.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Faniran, T.S., Falade, A.O. and Ogunsanwo, J. (2020) Sensitivity Analysis of an Untreated Liver-Stage Malaria. Journal of Computer Science & Computational Mathematics, 10, 61-68. https://www.jcscm.net/fp/179.pdf
[2]	Bakare, E.A. and Nwozo, C.R. (2015) Mathematical Analysis of the Dynamics of Malaria Disease Transmission Model. International Journal of Pure and Applied Mathematics, 99, 411-437. https://doi.org/10.12732/ijpam.v99i4.3
[3]	Tumwiine, J., Mugisha, J.Y.T. and Luboobi, L.S. (2008) Threshold and Stability Results for a Malaria Model in a Population with Protective Intervention among High-Risk Groups. Mathematical Modelling and Analysis, 13, 443-460. https://doi.org/10.3846/1392-6292.2008.13.443-460
[4]	Singh, K., Wester, W.C., Gordon, M. and Trenholme, G.M. (2003) Problems in the Therapy for Imported Malaria in the United States. Archives of Internal Medicine, 163, 2027-2030. https://doi.org/10.1001/archinte.163.17.2027
[5]	Ministry of Public Health and the Fight against HIV/AIDS (2018) Accelerate the Reduction of the Incidence of Malaria in Burundi, National Strategic Plan for the Fight against Malaria 2018-2023. Ministry of Public Health and the Fight against HIV/AIDS, Burundi.
[6]	Ross, R. (1910) The Prevention of Malaria. John Murray, London.
[7]	Labadin, J., Kon, C.M.L. and Juan, S.F.S. (2009) Deterministic Malaria Transmission Model with Acquired Immunity. World Congress on Engineering and Computer Science, Vol. 2, San Francisco, 20-22 October 2009, 779-784.
[8]	Georges, M. (1957) The Epidemiology and Control of Malaria. Oxford University Press, Oxford.
[9]	Chiyaka, C., Garira, W. and Dube, S. (2007) Transmission Model of Endemic Human Malaria in a Partially Immune Population. Mathematical and Computer Modelling, 46, 806-822. https://doi.org/10.1016/j.mcm.2006.12.010
[10]	Chitnis, N., Cushing, J.M. and Hyman, J.M. (2006) Bifurcation Analysis of a Mathematical Model for Malaria Transmission. SIAM Journal on Applied Mathematics, 67, 24-45. https://doi.org/10.1137/050638941
[11]	Li, J. (2008) A Malaria Model with Partial Immunity in Humans. Mathematical Biosciences and Engineering, 5, 789-801. https://doi.org/10.3934/mbe.2008.5.789
[12]	Van Den Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Treshold Endemic Equilibria for Compartmental of Disease Transmission. Mathematical Biosciences, 180, 29-48. https://doi.org/10.1016/S0025-5564(02)00108-6
[13]	Diouf, M.L. (2016) Analyse de modèles épidémiologiques à plusieurs classes d’infectés: stabilité et observabilité. Mathématiques [math]. Université Gaston BERGER de Saint-Louis, Sénégal.
[14]	https://perspective.usherbrooke.ca
[15]	Zongo, P. (2009) Mathematical Modeling of the Transmission Dynamics of Malaria Doctoral Thesis Ouagadougou University, Burkina Faso.
[16]	Olaniyi, S., Okosun, K., Adesanya, S. and Lebelo, R. (2020) Modelling Malaria Dynamics with Partial Immunity and Protected Travellers: Optimal Control and Cost-Effectiveness Analysis. Journal of Biological Dynamics, 14, 90-115. https://doi.org/10.1080/17513758.2020.1722265
[17]	Chitnis, N., Hyman, J.M. and Cushing, J.M. (2008) Determining Important Parameters in the Spread of Malaria through the Sensitivity Analysis of a Mathematical Model. Bulletin of Mathematical Biology, 70, Article No. 1272. https://doi.org/10.1007/s11538-008-9299-0
[18]	https://reliefweb.int
[19]	Ministry of Public Health and the Fight against HIV/AIDS (2018) Statistical Directory of Data from Health Centers and Hospitals, 2012-2018. Ministry of Public Health and the Fight against HIV/AIDS, Burundi.