A Mathematical Modelling of the Effect of Treatment in the Control of Malaria in a Population with Infected Immigrants

Olaniyi S. Maliki¹, Ngwu Romanus¹, Bruno O. Onyemegbulem²

¹Department of Mathematics, Michael Okpara University of Agriculture, Umudike, Nigeria.
²African Center of Excellence in Phytomedicine, Research and Development University of Jos, Jos, Nigeria.

DOI: 10.4236/am.2018.911081                

Abstract

In this work, we developed a compartmental bio-mathematical model to study the effect of treatment in the control of malaria in a population with infected immigrants. In particular, the vector-host population model consists of eleven variables, for which graphical profiles were provided to depict their individual variations with time. This was possible with the help of MathCAD software which implements the Runge-Kutta numerical algorithm to solve numerically the eleven differential equations representing the vector-host malaria population model. We computed the basic reproduction ratio R₀ following the next generation matrix. This procedure converts a system of ordinary differential equations of a model of infectious disease dynamics to an operator that translates from one generation of infectious individuals to the next. We obtained R₀ = , i.e., the square root of the product of the basic reproduction ratios for the mosquito and human populations respectively. R_0m explains the number of humans that one mosquito can infect through contact during the life time it survives as infectious. R_0h on the other hand describes the number of mosquitoes that are infected through contacts with the infectious human during infectious period. Sensitivity analysis was performed for the parameters of the model to help us know which parameters in particular have high impact on the disease transmission, in other words on the basic reproduction ratio R₀.

Keywords

Malaria Control, Infected Immigrants, Basic Reproduction Ratio, Differential Equations, MathCAD Simulation

1. Introduction

Malaria is a highly prevalent infectious disease especially in the tropical and subtropical areas. Figure 1 below is a map obtained from WHO Malaria Report 2010 [1] , depicting the countries where malaria was endemic in 2009 (shaded region).

In addition to being widespread, malaria is also a deadly disease. This is because statistics has shown that for Africa in particular, annually 145,000 million to 150,000 million infections are reported, among which, 800 to 850 cases result in deaths as shown in Table 1. Most of the deaths are either children under five or pregnant women. Typical symptoms of malaria infections start with headache, followed by periodic bouts of fevers and chills, and sometimes even coma. The period of cyclical fevers lasts several days, during which time a high probability of dying has been observed for children, since their immune systems are weak. Such fever can also lead to abortions in pregnant women.

1.1. Brief Analysis of Malaria Data

It is interesting to do a quick statistical analysis of the data in Table 1, for the malaria cases in Africa as provided by WHO (Figure 2). We perform a nonlinear regression analysis for both the reported cases (C) and deaths (D) against time (T). The result follows from SPSS.

Figure 1. Malaria endemic countries 2009.

Table 1. Estimates of malaria cases and deaths in Africa by WHO, 2000-2009.

Figure 2. Quadratic regression model for malaria cases 2000-2009.

Model Summary and Parameter Estimates. Dependent Variable: C (Numbers of cases)

The independent variable is T.

Observation: It is quite clear from the WHO data, for the number of malaria cases reported over the 10 year period that the incidence of malaria infection follows a parabolic curve, rising sharply initially, to reach a maximum and then declining sharply thereafter (Figure 3). The equation of the parabola is given by:

$C=165283.3+7609.85T-647.73T^2$ with goodness of fit $R^2=0.981$ .

Model Summary and Parameter Estimates. Dependent Variable: D (Number of deaths)

The independent variable is T.

Observation: The number of malaria related deaths over the 10 year period as depicted in the above graph, follows a parabolic curve, rising from a high value initially, then reaching a maximum and then declining sharply thereafter. The equation of the parabola is given by:

$D=882.883+12.07T-2.89T^2$ with goodness of fit $R^2=0.992$ .

Figure 3. Quadratic regression model for deaths caused by malaria 2000-2009.

1.2. Life Cycle of Malaria Parasites

Malaria is a vector-borne disease [2] . Malaria parasites are transferred between humans through mosquitoes. The malaria parasite life cycle is divided into two parts, one is within host (human) body and the other is within vector (mosquito) body.

Human infection starts from a blood meal of an infectious female mosquito. The parasites existing in the infectious mosquito’s saliva, called sporozoites at this stage, enter the bloodstream of the human through mosquito bites and migrate to the liver. Within minutes after entering in the human body, sporozoites infect hepatocytes, and multiply asexually and asymptomatically in liver cells for a period of 5 - 30 days [3] . This period is called the exo-erythrocytic stage. At the end of this stage, thousands of merozoites (schizonts) emerge inside an infected liver cell. These merozoites rupture their host cells undetectably by wrapping themselves in the membrane of infected liver cells. Then, merozoites escape into the bloodstream and get ready to infect red blood cells. Once entering the bloodstream, free merozoites undergo the so-called erythrocytic stage, in which merozoites invade red blood cells to develop ring forms before experiencing asexual or sexual maturation. Within the red blood cells, a proportion of parasites keep multiplying asexually and periodically break out of infected old red blood cells to invade fresh red blood cells. Such amplification cycles may cause the symptom of waves of fever. The remaining parasites follow sexual maturation and produce male (micro-) and female (macro-) gametocytes which may be taken up by bites of female mosquitoes. Finally, when it has developed into an infectious form, it spreads the disease to a new mosquito that bites the infectious human.

1.3. Malaria Control and Treatments

According to the transmission procedure of malaria, there are three conditions for the prevalence of the disease:

1) High density of Anopheles mosquitoes,

2) High density of human population,

3) Large rate of transmission of parasites between human beings and mosquitoes.

Obviously, not too much can be done in respect to (2). So, (1) and (3) are naturally targeted. That is, either controlling the population of Anopheles female mosquitoes at a lower level, or avoiding biting by mosquitoes can reduce the chance of malaria becoming endemic. In the middle of the last century, people in Africa have already knew how to remove or poison the breeding grounds of mosquitoes or the aquatic habitats of the larva stages, such as by filling or applying oil to places with standing water, to control the population of mosquitoes [4] . Later, pesticide was widely employed to eliminate mosquitoes. On the other hand, mosquito nets, bedclothes and mosquito-repellent incense (indoor residual spraying) also help to keep mosquitoes far away from people and minimize the biting rate, greatly reducing the chance of infection and transmission of malaria. There are some effective drugs for malaria patients currently. For example, Chloroquine, Quinuine, Primaquine and combinations of some other drugs like sulfadoxine and pyrimethamine (SP) are effective medicines for treating infections caused by the five major parasites. Although malaria is an entirely preventable or curable disease thanks to these effective medicines, there are still millions of people suffering from this disease, who are too poor to afford full treatments. Moreover, insufficient treatments due to poor economic conditions, may result in drug resistance and lead to emergence of new (drug resistant) strains of malaria parasites. For instance, the first case of resistance to Chloroquine was documented in 1957. Chloroquine, Quinine and Sulfadoxine-pyrimethamine resistance cases have been reported in almost all disease endemic areas [5] .

1.4. Control of Mosquito-Borne Infections

In order to control mosquito-borne infections one can adopt the following measures;

• Reduce vector population: Make environment less mosquito-friendly by draining stagnant water.

• Use insecticides; not without problems: for example some mosquitoes become insecticide resistant.

• Prevent mosquitoes biting people. Insecticide-laced bed nets, although this is ineffective against mosquitoes that mainly bite during the day (e.g. A. aegypti).

• Vaccines and drug treatments. Not always available, there are problems with drugs and drug resistance.

1.5. The Ross-Macdonald Malaria Model

The first and simplest model of malaria was developed by [6] Ross and later extended by Macdonald [7] . This so-called Ross-Macdonald model is the best-known and most widely used model. Despite its simple structure as shown below, it enables us to interpret and compare a broad range of epidemiological models.

1.6. Remark

In the Ross-Macdonald model of malaria transmission, the flow of human from a susceptible class to an infected class and through recovery from infection, the reverse is shown in the upper part of the Figure 4. The flow of mosquitoes from susceptible class to an infected class, and finally to an infectious class is shown further down. The human and mosquito population are linked through the transmission process.

1.7. Statement of the Problem

The development of the means intended to reduce the spread of malaria infections and eradication necessitates decisive measures to curb the malaria epidemic. In particular, sustained minimization of the number of humans with incidence of malaria as a result of adequate control, can be attained by developing a suitable mathematical model which can enable us to understand better the dynamics and control of the vector-host endemic.

In developing the model, the human population is compartmentalized into seven classes including the susceptible, infected, exposed, treated, non-treated, recovered, and protected classes. For the mosquito population, we have four classes, namely; class of mosquito larva, susceptible mosquitoes, infected mosquitoes and exposed mosquitoes. We assume free interaction between the vector and host populations. The mathematical analysis of the compartmental models leads us to eleven coupled systems of nonlinear ordinary differential equations.

2. Construction of the Compartmental Model

In this section we develop a compartmental bio-mathematical model (Figure 5)

Figure 4. The ross-macdonald malaria model.

Figure 5. Compartmental model for human-mosquito interaction.

to study the effect of treatment in the control of malaria in a population with infected immigrants.

From the above compartmental model we obtain the following equations for the dynamics of the human-mosquito interaction.

2.1. Human Population

$\frac{\text{d}{S}_H}{\text{d}t}=\left(1-q\right){\Lambda }_H+{\alpha }_{2}A+\rho R_H-\frac{{\beta }_H{I}_M{S}_H}{N_H}-{\alpha }_1{S}_H-{\delta }_H{S}_H$ (1)

$\frac{\text{d}{E}_H}{\text{d}t}=\frac{{\beta }_H{I}_M{S}_H}{N_H}-g{E}_H-{\delta }_H{E}_H$ (2)

$\frac{\text{d}{I}_H}{\text{d}t}=g{E}_H+q{\Lambda }_H-k_1{I}_H-k_2{I}_H-{\delta }_H{I}_H$ (3)

$\frac{\text{d}{I}_{HN}}{\text{d}t}=k_2{I}_H-\left({\omega }_H+{\delta }_H\right)I_{HN}$ (4)

$\frac{\text{d}{T}_H}{\text{d}t}=k_1{I}_H-\gamma T_H-{\delta }_H{T}_H$ (5)

$\frac{\text{d}{R}_H}{\text{d}t}=\gamma T_H-\left(\mu +\rho +{\delta }_H\right)R_H$ (6)

$\frac{\text{d}A}{\text{d}t}={\alpha }_1{S}_H+\mu R_H-{\alpha }_{2}A-{\delta }_{H}A$ (7)

2.2. Mosquitoe Population

$\frac{\text{d}{L}_M}{\text{d}t}={\Lambda }_M-m{L}_M-{\delta }_M{L}_M$ (8)

$\frac{\text{d}{S}_M}{\text{d}t}=m{L}_M-\frac{{\beta }_M{I}_H{S}_M}{N_H}-{\delta }_M{S}_M$ (9)

$\frac{\text{d}{E}_M}{\text{d}t}=\frac{{\beta }_M{I}_H{S}_M}{N_H}-\varphi E_M-{\delta }_M{E}_M$ (10)

$\frac{\text{d}{I}_M}{\text{d}t}=\varphi E_M-{\delta }_M{I}_M$ (11)

2.3. Remark

The state variables and parameters are defined in Table 2 and Table 3 respectively.

Table 2. State variables of the basic malaria model.

Table 3. Parameters of the basic malaria model.

2.4. Invariant Region

The total population sizes $N_H$ and $N_M$ can be determined by $N_H=S_H+E_H+I_H+I_{HN}+T_H+A+R_H$ and $N_M=L_M+S_M+E_M+I_M$ . Thus

$\frac{\text{d}{N}_H}{\text{d}t}={\Lambda }_H-{\delta }_H{N}_H-{\omega }_H{I}_{HN}$ (12)

Without loss of generality, we can write

$\frac{\text{d}{N}_H}{\text{d}t}\le {\Lambda }_H-{\delta }_H{N}_H,\frac{\text{d}{N}_M}{\text{d}t}\le {\Lambda }_M-{\delta }_M{N}_M$ (13)

2.5. Lemma

The model system has solution which are contained in the feasible $\Omega ={\Omega }_H\times {\Omega }_M$ .

Proof: let $\Omega =\left\{S_H,E_H,I_H,I_{HN},T_H,A,R_H,L_M,S_M,E_M,I_M\right\}\in {ℝ}_{+}^{11}$ be any solution of the system with non-negative initial conditions. From Equation (13)

$\frac{\text{d}{N}_H}{\text{d}t}\le {\Lambda }_H-{\delta }_H{N}_H$ (14)

Adopting Birhoff and Rotta [8] theorem on differential inequality, we have

$0\le N_H\le \frac{{\Lambda }_H}{{\delta }_H},{\Lambda }_H-{\delta }_H{N}_H\ge C{\text{e}}^{-{\delta }_{H}t}$ (15)

where C is a constant.

Therefore, all feasible solutions of the human population only of the model system are in the region.

${\Omega }_H=\left\{\left(S_H,E_H,I_H,I_{HN},T_H,A,R_H\right)\in {ℝ}_{+}^7:N_H\le \frac{{\Lambda }_H}{{\delta }_H}\right\}$

Similarly the feasible set for model of the mosquitoes population only are in the region

${\Omega }_M=\left\{\left(L_M,S_M,E_M,I_M\right)\in {ℝ}_{+}^4:N_M\le \frac{{\Lambda }_M}{{\delta }_M}\right\}$

Therefore the feasible set for the model system is given by

$\begin{array}{l}\Omega =\{\left(S_H,E_H,I_H,I_{HN},T_H,A,R_H,L_M,S_M,E_M,I_M\right)\in {ℝ}_{+}^{11}:N_H\le \frac{{\Lambda }_H}{{\delta }_H}=N_H^{*},\\ N_M\le \frac{{\Lambda }_M}{{\delta }_M}=N_M^{*}\}\end{array}$ (16)

2.6. Mathematical Analysis of the Model

The nonlinear system (1)-(11) will be qualitatively analyzed so as to find the conditions for existence and stability of disease free equilibrium points. Analysis of the model allows us to determine the impact of treatment on the transmission of malaria infection in a population. Also on finding the reproductive number $R_0$ , one can determine if the disease become endemic in a population or not [9] . However, one can see that adding the human equation of the model, with the case that there is no disease -induced death. From Equation (13)

$\frac{\text{d}{N}_H}{\text{d}t}={\Lambda }_H-{\delta }_H{N}_H$ , hence $N_H\left(t\right)\to \frac{{\Lambda }_H}{{\delta }_H}$ as $t\to \infty$ .

Thus $\frac{{\Lambda }_H}{{\delta }_H}$ is the upper bound of $N_H\left(t\right)$ provided that $N_H\left(0\right)\le \frac{{\Lambda }_H}{{\delta }_H}$ . Similarly,

$\frac{\text{d}{N}_M}{\text{d}t}={\Lambda }_M-{\delta }_M{N}_M\Rightarrow N_M\left(t\right)\to \frac{{\Lambda }_M}{{\delta }_M}$ as $t\to \infty$ .

Thus $\frac{{\Lambda }_M}{{\delta }_M}$ is the upper bound of $N_M\left(t\right)$ provided that $N_M\left(0\right)\le \frac{{\Lambda }_M}{{\delta }_M}$ . Hence the invariant region is

$\begin{array}{l}\Omega =\{\left(S_H,E_H,I_H,I_{HN},T_H,A,R_H,L_M,S_M,E_M,I_M\right)\in {ℝ}_{+}^{11}:N_H\le \frac{{\Lambda }_H}{{\delta }_H}=N_H^{*},\\ N_M\le \frac{{\Lambda }_M}{{\delta }_M}=N_M^{*}\}\end{array}$

is positively invariant. Hence no solution path leaves through and boundary of $\Omega$ . Since path cannot leave $\Omega$ , solution remains non-negative for non negative initial conditions. This means that the solution exists for all positive time t. Therefore the model (1)-(11) is mathematically and epidemiological well-posed [10] .

For convenience and to simplify the analysis of our model, we rewrite the model system (1)-(11) in terms of the proportions of individual in each class. Let

$\begin{array}{l}{s}_h=\frac{S_H}{N_H},e_h=\frac{E_H}{N_H},i_h=\frac{I_H}{N_H},t_h=\frac{T_H}{N_H},r_h=\frac{R_H}{N_H},i_{hn}=\frac{I_{HN}}{N_H},\\ z=\frac{A}{N_H},l_m=\frac{L_M}{N_H},s_m=\frac{S_M}{N_H},e_m=\frac{E_M}{N_H},i_m=\frac{I_M}{N_H}.\end{array}$

Let $\pi =\frac{N_M}{N_H}$ be the female mosquito?human ratio, that is, the number of female mosquito per human host. The ratio $\pi =\frac{N_M}{N_H}$ is constant because a

mosquito takes a fixed number of blood meals per unit independent of the population density of the host [11] . Also let

${\Lambda }_H={\Lambda }_h,{\Lambda }_M={\Lambda }_m,{\beta }_H={\beta }_h,{\delta }_H={\delta }_h,{\beta }_M={\beta }_m,{\delta }_M={\delta }_m,{\omega }_H={\omega }_h.$

The simplified model now becomes modified human and mosquito population models.

2.7. Modified Human Population

$\frac{\text{d}{s}_h}{\text{d}t}=\left(1-q\right){\Lambda }_h+{\alpha }_{2}z+\rho r_h-{\beta }_h{i}_m{s}_h-{\alpha }_1{s}_h-{\delta }_h{s}_h$ (17)

$\frac{\text{d}{e}_h}{\text{d}t}={\beta }_h{i}_m{s}_h-g{e}_h-{\delta }_h{e}_h$ (18)

$\frac{\text{d}{i}_h}{\text{d}t}=g{e}_h+q{\Lambda }_h-k_1{i}_h-k_2{i}_h-{\delta }_h{i}_h$ (19)

$\frac{\text{d}{i}_{hn}}{\text{d}t}=k_2{i}_h-\left({\omega }_h+{\delta }_h\right)i_{hn}$ (20)

$\frac{\text{d}{t}_h}{\text{d}t}=k_1{i}_h-\gamma t_h-{\delta }_h{t}_h$ (21)

$\frac{\text{d}{r}_h}{\text{d}t}=\gamma t_h-\left(\mu +\rho +{\delta }_h\right)r_h$ (22)

$\frac{\text{d}z}{\text{d}t}={\alpha }_1{s}_h+\mu r_h-{\alpha }_{2}z-{\delta }_{h}z$ (23)

2.8. Modified Mosquitoes Population

$\frac{\text{d}{l}_m}{\text{d}t}={\Lambda }_m-m{l}_m-{\delta }_m{l}_m$ (24)

$\frac{\text{d}{s}_m}{\text{d}t}=m{l}_m-{\beta }_m{i}_h{s}_m-{\delta }_m{s}_m$ (25)

$\frac{\text{d}{e}_m}{\text{d}t}={\beta }_m{i}_h{s}_m-\varphi e_m-{\delta }_m{e}_m$ (26)

$\frac{\text{d}{i}_m}{\text{d}t}=\varphi e_m-{\delta }_m{i}_m$ (27)

2.9. Positivity of Solutions

It is necessary to prove that all solutions of system (17)-(27) with positive initial data will remain positive for all times $t>0$ . This will be established by the following theorem.

2.10. Theorem

Let the initial data be

$\begin{array}{l}\{s_h\left(0\right)\ge 0,i_h\left(0\right)\ge 0,i_{hn}\left(0\right)\ge 0,t_h\left(0\right)\ge 0,z\left(0\right)\ge 0,r_h\left(0\right)\ge 0,\\ e_h\left(0\right)\ge 0,s_m\left(0\right)\ge 0,l_m\left(0\right)\ge 0,e_m\left(0\right)\ge 0,i_m\left(0\right)\ge 0\}\in \Omega \end{array}$

Then the solution set $\left(s_h,e_h,i_h,i_{hn},t_h,z,r_h,l_m,s_m,e_m,i_m\right)\left(t\right)$ of the model system (4) is positive for all $t>0$ .

Proof: From first equation of (17)

$\frac{\text{d}{s}_h}{\text{d}t}=\left(1-q\right){\Lambda }_h+{\alpha }_{2}z+\rho r_h-{\beta }_h{i}_m{s}_h-{\alpha }_1{s}_h-{\delta }_h{s}_h\ge -\left({\beta }_h{i}_m+{\alpha }_1+{\delta }_h\right)s_h$

$\Rightarrow {\displaystyle \int \frac{1}{s_h}\text{d}\left(s_h\right)}\ge -{\displaystyle \int \left({\beta }_h{i}_m+{\alpha }_1+{\delta }_h\right)\text{d}t}$

$\therefore s_h\left(t\right)\ge s_h\left(0\right){\text{e}}^{-\left({\beta }_h{i}_m+{\alpha }_1+{\delta }_h\right)t}\ge 0$

Following the above procedure, from equations (18)-(23), we obtain respectively the positivity conditions;

$\begin{array}{l}{e}_h\left(t\right)\ge e_h\left(0\right){\text{e}}^{-\left(g+{\delta }_h\right)t}\ge 0,i_h\left(t\right)\ge i_h\left(0\right){\text{e}}^{-\left(k_1+k_2+{\delta }_h\right)t}\ge 0,\\ i_{hn}\left(t\right)\ge i_{hn}\left(0\right){\text{e}}^{-\left({\omega }_h+{\delta }_h\right)t}\ge 0,t_h\left(t\right)\ge t_h\left(0\right){\text{e}}^{-\left(\gamma +{\delta }_h\right)t}\ge 0,\\ r_h\left(t\right)\ge r_h\left(0\right){\text{e}}^{-\left(\mu +\rho +{\delta }_h\right)t}\ge 0,z\left(t\right)\ge z\left(0\right){\text{e}}^{-\left({\delta }_h+{\alpha }_2\right)}{}^t\ge 0.\end{array}$

Similarly for the modified mosquito population, equations (20)-(27) gives the positivity conditions;

$\begin{array}{l}{l}_m\left(t\right)\ge l_m\left(0\right){\text{e}}^{-\left(m+{\delta }_m\right)t}\ge 0,s_m\left(t\right)\ge s_m\left(0\right){\text{e}}^{-\left({\beta }_m{i}_h+{\delta }_m\right)t}\ge 0,\\ e_m\left(t\right)\ge e_m\left(0\right){\text{e}}^{-\left(\varphi +{\delta }_m\right)t}\ge 0,i_m\left(t\right)\ge i_m\left(0\right){\text{e}}^{-{\delta }_{m}t}\ge 0.\end{array}$

2.11. Existence and Stability of Steady-State Solutions

Let $E^0=\left(s_h^0,e_h^0,i_h^0,i_{hn}^0,t_h^0,z^0,r_h^0,l_m^0,s_m^0,e_m^0,i_m^0\right)$ be the steady-state of the system (17)-(27) which can be calculated by setting the right hand side of the model (17)-(27) to zero, giving us the following;

$\left(1-q\right){\Lambda }_h+{\alpha }_{2}z+\rho r_h-{\beta }_h{i}_m{s}_h-{\alpha }_1{s}_h-{\delta }_h{s}_h=0$ (28)

${\beta }_h{i}_m{s}_h-g{e}_h-{\delta }_h{e}_h=0$ (29)

$g{e}_h+q{\Lambda }_h-k_1{i}_h-k_2{i}_h-{\delta }_h{i}_h=0$ (30)

$k_2{i}_h-\left({\omega }_h+{\delta }_h\right)i_{hn}=0$ (31)

$k_1{i}_h-\gamma t_h-{\delta }_h{t}_h=0$ (32)

$\gamma t_h-\left(\mu +\rho +{\delta }_h\right)r_h=0$ (33)

${\alpha }_1{s}_h+\mu r_h-{\alpha }_{2}z-{\delta }_{h}z=0$ (34)

${\Lambda }_m-m{l}_m-{\delta }_m{l}_m=0$ (35)

$m{l}_m-{\beta }_m{i}_h{s}_m-{\delta }_m{s}_m=0$ (36)

${\beta }_m{i}_h{s}_m-\varphi e_m-{\delta }_m{e}_m=0$ (37)

$\varphi e_m-{\delta }_m{i}_m=0$ (38)

2.12. Disease-Free Equilibrium Point

Disease-free equilibrium points (DFE) are steady-state solutions where there is no disease (malaria). The disease free equilibrium of the normalized model (17)-(27) is obtained by setting

$\frac{\text{d}{s}_h}{\text{d}t}=\frac{\text{d}{\stackrel{\dot{}}{e}}_h}{\text{d}t}=\frac{\text{d}{i}_h}{\text{d}t}=\frac{\text{d}{i}_{hn}}{\text{d}t}=\frac{\text{d}{t}_h}{\text{d}t}=\frac{\text{d}{r}_h}{\text{d}t}=\frac{\text{d}z}{\text{d}t}=\frac{\text{d}{l}_m}{\text{d}t}=\frac{\text{d}{s}_m}{\text{d}t}=\frac{\text{d}{e}_m}{\text{d}t}=\frac{\text{d}{i}_m}{\text{d}t}=0$

At disease free equilibrium we have,

$\begin{array}{l}{s}_h=\frac{{\Lambda }_h}{{\delta }_h},s_m=\frac{m{\Lambda }_m}{{\delta }_m\left(m+{\delta }_m\right)},\\ e_h=i_h=i_{hn}=t_h=r_h=l_m=e_m=i_m=z=q=0.\end{array}$

Therefore the disease free equilibrium (DFE) denoted by $E^0$ of the system (28)-(38) is given by

$\begin{array}{c}{E}^0=\left(s_h^0,e_h^0,i_h^0,i_{hn}^0,t_h^0,z^0,r_h^0,l_m^0,s_m^0,e_m^0,i_m^0\right)\\ =\left(\frac{{\Lambda }_h}{{\delta }_h},0,0,0,0,0,0,0,\frac{m{\Lambda }_m}{{\delta }_m\left(m+{\delta }_m\right)},0,0\right)\end{array}$

that represents the state in which there is no infection in the society and is known as the disease-free equilibrium point (DFE). This implies that at the disease-free equilibrium, the susceptible human population is equal to the total human population and the susceptible mosquito population is equal to the total mosquito population.

2.13. Local Stability of DFE

The disease free equilibrium of the model (17)-(27) was given by

$\begin{array}{c}{E}^0=\left(s_h^0,e_h^0,i_h^0,i_{hn}^0,t_h^0,z^0,r_h^0,l_m^0,s_m^0,e_m^0,i_m^0\right)\\ =\left(\frac{{\Lambda }_h}{{\delta }_h},0,0,0,0,0,0,0,\frac{m{\Lambda }_m}{{\delta }_m\left(m+{\delta }_m\right)},0,0\right)\end{array}$

2.14. Basic Reproduction Ratio

R₀ is often found through the study and computation of the eigenvalues of the Jacobian at the disease- or infectious-free equilibrium Diekmann [12] follow a different approach which is the next generation matrix method. This procedure converts a system of ordinary differential equations of a model of infectious disease dynamics to an operator (or matrix) that translate from one generation of infectious individuals to the next. The basic reproductive number is then defined as the spectral radius (dominant eigenvalue) of this operator. Van den Driessche and Watmough [9] describe such a method in detail for general deterministic compartmental models.

The dynamics of the model is specified by the IVP;

$\frac{\text{d}{x}_i}{\text{d}t}=f_i\left(x\right),x\left(0\right)\in {ℝ}_{+}^n$ (39)

We define ${\Theta }_0$ as the set of all disease-free states as

${\Theta }_0=\left\{x\in {ℝ}_{+}^n:x_i=0,1\le i\le m\right\}$ (40)

Next we recast the IVP (4.39) in the form;

$\frac{\text{d}{x}_i}{\text{d}t}=F_i\left(x\right)-V_i\left(x\right)$ (41)

where $F_i\left(x\right)$ is the rate of new infections entering compartment i, and

$V_i=V_i^{-}\left(x\right)-V_i^{+}\left(x\right)$ (42)

where $V_i^{+}\left(x\right)$ is the rate of transfer into compartment i by any other means, and $V_i^{-}\left(x\right)$ is the rate of transfer out of compartment i. Given a disease-free equilibrium point $x_{DFE}$ of (39), with $x_{DFE}$ and $f\left(x\right)$ satisfying certain important assumptions [12] , then we define the square matrices F and V of dimension $m\times m$ as follows;

$F_{ij}={\frac{\partial F_i\left(x\right)}{\partial x_j}|}_{x_{DFE}},V_{ij}={\frac{\partial V_i\left(x\right)}{\partial x_j}|}_{x_{DFE}}\text{for}1\le i,j\le m$ (43)

It then follows that $F{V}^{-1}$ is the next generation matrix and the basic reproduction ratio $R_0$ is the spectral radius of $F{V}^{-1}$ ,

$\Rightarrow R_0=\rho \left(F{V}^{-1}\right)$ (44)

Rewriting the system (41) starting with the infected compartments for both populations; $e_h,i_h,e_m,i_m,i_{hn},t_h$ and then followed by uninfected classes; $s_h,z,r_h,l_m,s_m$ also from the two populations, gives;

$\frac{\text{d}{e}_h}{\text{d}t}={\beta }_h{i}_m{s}_h-g{e}_h-{\delta }_h{e}_h$

$\frac{\text{d}{i}_h}{\text{d}t}=g{e}_h+q{\Lambda }_h-k_1{i}_h-k_2{i}_h-{\delta }_h{i}_h$

$\frac{\text{d}{e}_m}{\text{d}t}={\beta }_m{i}_h{s}_m-\varphi e_m-{\delta }_m{e}_m$

$\frac{\text{d}{i}_m}{\text{d}t}=\varphi e_m-{\delta }_m{i}_m$

$\frac{\text{d}{i}_{hn}}{\text{d}t}=k_2{i}_h-\left({\omega }_h+{\delta }_h\right)i_{hn}$

$\frac{\text{d}{t}_h}{\text{d}t}=k_1{i}_h-\gamma t_h-{\delta }_h{t}_h$

$\frac{\text{d}{s}_h}{\text{d}t}=\left(1-q\right){\Lambda }_h+{\alpha }_{2}z+\rho r_h-{\beta }_h{i}_m{s}_h-{\alpha }_1{s}_h-{\delta }_h{s}_h$

$\frac{\text{d}z}{\text{d}t}={\alpha }_1{s}_h+\mu r_h-{\alpha }_{2}z-{\delta }_{h}z$

$\frac{\text{d}{r}_h}{\text{d}t}=\gamma t_h-\left(\mu +\rho +{\delta }_h\right)r_h$

$\frac{\text{d}{l}_m}{\text{d}t}={\Lambda }_m-m{l}_m-{\delta }_m{l}_m$

$\frac{\text{d}{s}_m}{\text{d}t}=m{l}_m-{\beta }_m{i}_h{s}_m-{\delta }_m{s}_m$

The method of next generation matrix has been used to show the rate of appearance of new infection in compartments; $e_h$ and $e_m$ , from the system (12);

$F=\left(\begin{array}{c}{\beta }_h{i}_m{s}_h\\ 0\\ {\beta }_m{i}_h{s}_m\\ 0\\ 0\\ 0\end{array}\right),V=\left(\begin{array}{c}\left(g+{\delta }_h\right)e_h\\ -g{e}_h-q{\Lambda }_h+\left(k_1+k_2+{\delta }_h\right)i_h\\ \left(\varphi +{\delta }_m\right)e_m\\ -\varphi e_m+{\delta }_m{i}_m\\ -k_2{i}_h+\left({\omega }_h+{\delta }_h\right)i_{hn}\\ -k_1{i}_h+\left(\gamma +{\delta }_h\right)t_h\end{array}\right)$

By linearization approach, the associated matrix at disease free equilibrium is obtained as

$F=\left(\begin{array}{cccccc}0& 0& 0& \frac{{\Lambda }_h{\beta }_h}{{\delta }_h}& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& \frac{m{\Lambda }_m{\beta }_m}{{\delta }_m\left(m+{\delta }_m\right)}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)$

$F{V}^{-1}=\left(\begin{array}{cccccc}0& 0& \frac{{\Lambda }_h{\beta }_h\varphi }{{\delta }_h{\delta }_m\left(\varphi +{\delta }_m\right)}& \frac{{\Lambda }_h{\beta }_h}{{\delta }_h{\delta }_m}& 0& 0\\ 0& 0& 0& 0& 0& 0\\ \frac{m{\Lambda }_m{\beta }_{m}g}{{\delta }_m\left(m+{\delta }_m\right)\left(g+{\delta }_h\right)\left(k_1+k_2+{\delta }_h\right)}& \frac{m{\Lambda }_m{\beta }_m}{{\delta }_m\left(m+{\delta }_m\right)k_1+k_2+{\delta }_h}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)$

$\therefore R_0=\sqrt{\left(\frac{{\Lambda }_h{\beta }_h\varphi }{{\delta }_h{\delta }_m\left(\varphi +{\delta }_m\right)}\right)\left(\frac{m{\Lambda }_m{\beta }_{m}g}{{\delta }_m\left(m+{\delta }_m\right)\left(g+{\delta }_h\right)\left(k_1+k_2+{\delta }_h\right)}\right)}$

Here the term $\frac{{\Lambda }_h{\beta }_h\varphi }{{\delta }_h{\delta }_m\left(\varphi +{\delta }_m\right)}$ explains the number of humans that one mosquito infect through contact during the life time it survives as infectious. On the other hand $\frac{m{\Lambda }_m{\beta }_{m}g}{{\delta }_m\left(m+{\delta }_m\right)\left(g+{\delta }_h\right)\left(k_1+k_2+{\delta }_h\right)}$ describes the number of mos-

quitoes that are infected through contacts with the infectious human during infectious period. Hence

$R_0=\sqrt{R_{0m}\times R_{0h}}$

where $R_{0m}=\frac{{\Lambda }_h{\beta }_h\varphi }{{\delta }_h{\delta }_m\left(\varphi +{\delta }_m\right)}$ and $R_{0h}=\frac{m{\Lambda }_m{\beta }_{m}g}{{\delta }_m\left(m+{\delta }_m\right)\left(g+{\delta }_h\right)\left(k_1+k_2+{\delta }_h\right)}$ .

3. Sensitivity Analysis of the Model Parameters

In this section, we carry out the sensitivity analysis of the model parameter to help us know the parameters that have high impact on the disease transmission, which is on the reproduction ratio $R_0$ .

We used the normalized forward sensitivity index of a variable to parameter approach used in Okosun [13] .

3.1. Sensitivity Analysis of R₀

We compute the sensitivity of $R_0$ to each of the parameters described in Table 4. Using the formula

${\gamma }_n^m=\frac{\partial m}{\partial n}\times \frac{n}{m}$

where n represents the variables of the model, and m the parameters.

Sensitivity index of $\varphi$ given by $-\frac{1}{2}\left(\frac{\varphi }{\varphi +{\delta }_m}\right)$

Table 4. Sensitivity index of parameters.

Sensitivity index of g given by $-\frac{1}{2}\left(\frac{g}{g+{\delta }_h}\right)$

Sensitivity index of ${\delta }_m$ given by $\frac{1}{2}\left(-2-\frac{{\delta }_m}{\varphi +{\delta }_m}-\frac{{\delta }_m}{m+{\delta }_m}\right)$

Sensitivity index of ${\delta }_h$ given by $\frac{1}{2}\left(-1-\frac{{\delta }_h}{g+{\delta }_h}-\frac{{\delta }_h}{k_1+k_2+{\delta }_h}\right)$

Sensitivity index of $k_1$ given by $-\frac{1}{2}\left(\frac{k_1}{k_1+k_2+{\delta }_h}\right)$

Sensitivity index of $k_2$ given by $-\frac{1}{2}\left(\frac{k_2}{k_1+k_2+{\delta }_h}\right)$

Sensitivity index of m given by $\frac{1}{2}\left(1-\frac{m}{m+{\delta }_m}\right)$

Sensitivity index of ${\Lambda }_m={\beta }_m={\beta }_h={\Lambda }_h=\frac{1}{2}$

Remark: Sensitivity indices of R₀ evaluated at the baseline parameter values are given in the Table 5.

From Table 5, the sensitivity index may be a complex expression, depending on different parameters of the system. But it can also be a constant value. Example, the sensitivity index of ${\beta }_M$ , ${\beta }_H$ = +0.5, means that increasing (or decreasing) ${\beta }_M$ , ${\beta }_H$ by 10% increases (or decreases) R₀ by 5%.

3.2. Math Cad Simulation of the Model

Parameter values:

$\begin{array}{l}q:=0.1,{\Lambda }_H:=0.5,{\alpha }_1:=0.8,{\alpha }_2:=0.6,\rho :=0.02,{\beta }_H:=0.5,{\delta }_H:=0.3,\\ g_1:=0.9k_1:=0.8k_2:=0.5,{\omega }_H:=0.5,\gamma :=0.7,\mu :=0.4,N_H:=100,\\ {\Lambda }_M:=0.2,m_1:=0.3,{\delta }_M:=0.1,{\beta }_M:=0.4,\varphi :=0.12\end{array}$

$D\left(t,Y\right):=\left[\begin{array}{c}\left(1-q\right){\Lambda }_H+{\alpha }_2{Y}_6+\rho Y_5-\frac{{\beta }_H{Y}_{10}{Y}_0}{N_H}-{\alpha }_1{Y}_0-{\delta }_H{Y}_0\\ \frac{{\beta }_H{Y}_{10}{Y}_0}{N_H}-g_1{Y}_1-{\delta }_H{Y}_1\\ g_1{Y}_1+q_1{\Lambda }_H-k_1{Y}_2-k_2{Y}_2-{\delta }_H{Y}_2\\ k_2{Y}_2-\left({\omega }_H+{\delta }_H\right)Y_3\\ k_1{Y}_2-\gamma Y_4-{\delta }_H{Y}_4\\ \gamma Y_4-\left(\mu +\rho +{\delta }_H\right)Y_5\\ {\alpha }_1{Y}_0+\mu Y_5-{\alpha }_2{Y}_6-{\delta }_H{Y}_6\\ {\Lambda }_M-m_1{Y}_7-{\delta }_M{Y}_7\\ m_1{Y}_7-\frac{{\beta }_M{Y}_8{Y}_2}{N_H}-{\delta }_M{Y}_8\\ \frac{{\beta }_M{Y}_8{Y}_2}{N_H}-\varphi Y_9-{\delta }_M{Y}_9\\ \varphi Y_9-{\delta }_M{Y}_{10}\end{array}\right]$

Vector of derivative values at any solution point (t, Y):

Define additional arguments for the ODE solver:

$t0:=0$ : Initial value of independent variable

$t1:=0$ : Initial value of independent variable

$Y0:={\left[\begin{array}{ccccccccccc}50& 15& 25& 2& 2& 4& 2& 5& 3& 2& 1\end{array}\right]}^{\text{T}}$ : Vector of initial function values

$num:=1\times {10}^3$ : Number of solution values on [t0, t1]

$S1:=\text{Rkadapt}\left(Y0,t0,t1,num,D\right)$ : Solution matrix

Human (Table 6)

$t:=S{1}^{\langle 0\rangle }$ : Independent variable values

$S_H:=S{1}^{\langle 1\rangle }$ : First solution function values

$E_H:=S{1}^{\langle 2\rangle }$ : Second solution function values

$I_H:=S{1}^{\langle 3\rangle }$ : Third solution function values

$I_{HN}:=S{1}^{\langle 4\rangle }$ : Fourth solution function values

$T_H:=S{1}^{\langle 5\rangle }$ : Fifth solution function values

$R_H:=S{1}^{\langle 6\rangle }$ : Sixth solution function values

$A_H:=S{1}^{\langle 7\rangle }$ : Seventh solution function values

Table 5. Sensitivity indices of R₀ evaluated at the baseline parameter values.

Table 6. Solution matrix S1 for the system of ODEs

Mosquitoes

$L_M:=S{1}^{\langle 8\rangle }$ : Eighth solution function values

$S_M:=S{1}^{\langle 9\rangle }$ : Ninth solution function values

$E_M:=S{1}^{\langle 10\rangle }$ : Tenth solution function values

$I_M:=S{1}^{\langle 11\rangle }$ : Eleventh solution function values

3.3. Results and Discussion

The susceptible human population $S_H$ against time (Figure 6(a)), clearly shows a rapid exponential decline from the initial value to zero. Similarly, the variation of exposed human population $E_H$ against time (Figure 6(b)) depicts an exponential decline from the initial value to zero. The variation of the infected human population $I_H$ against time (Figure 6(c)), also depicts an exponential decline from the initial value to zero. The graphical profile of the variation of the non treated human population $I_{HN}$ against time (Figure 6(d)), shows a sharp rise from the initial value to reach a maximum, and thereafter exhibits an exponential decline to zero. The variation of treated human population $T_H$ against time (Figure 6(e)), shows a sharp rise from the initial value to reach a maximum, and thereafter declines exponentially to zero. Similarly, the graphical profile of the variation of the removed human population $R_H$ against time (Figure 6(f)), depicts a rise from the initial value to reach a maximum, and thereafter declines exponentially to zero. The graphical profile of the variation of the protected human population $A_H$ against time (Figure 6(g)), shows a sharp rise from the initial value to reach a maximum, and thereafter declines exponentially to a steady state. From the graphical profile of the variation of population of mosquito larva $L_M$ against time (Figure 6(h)), we observe an exponential decline from the initial value to reach a steady state. The variation of the susceptible mosquito population $S_M$ against time (Figure 6(i)), depicts a rise from the initial value to reach a maximum, and thereafter exhibits a sharp decline. In the same manner, the variation of the exposed mosquito population $E_M$ against time (Figure 6(j)), shows a decline from the initial value to reach a steady state. Finally, the variation of the infected mosquito population $I_M$ against time (Figure 6(k)), depicts a rise from the initial value to reach a maximum and then exhibits a decline.

3.4. Conclusion

Despite the availability of drugs, the malaria disease is still endemic in many parts of the world including developed countries. Elimination of malaria requires maintaining the effective reproduction number R₀ less than unity, as well as achieving low levels of susceptibility. In this research work, we developed a compartmental bio-mathematical model to study the effect of treatment in the control of malaria in a population with infected immigrants. We obtained the basic reproduction number, R₀ and studied the stability of the disease-free equilibrium of the model. Sensitivity analysis of R₀ with respect to the model parameters was carried out on the compartmental vector-host malaria model with

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k)

Figure 6. (a) Population of susceptible humans against time; (b) Population of exposed humans against time; (c) Population of infected humans against time; (d) Population of non-treated infected human against time; (e) Population of treated humans against time; (f) Population of recovered humans against time; (g) Population of protected humans against time; (h) Population of mosquitoes larva against time; (i) Population of susceptible mosquitoes against time; (j) Population of exposed mosquitoes against time; (k) Population of infected mosquitoes against time.

eleven compartments. From the literature on modelling of vector-host malaria models, we discovered that many researchers failed to consider protective measures in their models, though some discussed it theoretically. Our major contribution to the existing body of knowledge is incorporating the protective measure in our mathematical model.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	WHO (2010) Estimates of Malaria Cases and Deaths in Africa. 2000-2009. http://www.rollbackmalaria.org/about-malaria/key-facts
[2]	Matson, A. (1957) The History of Malaria in Nandi. East African Medical Journal, 34, 431-441.
[3]	Johansson, P. and Leander, J. (2010) Mathematical Modeling of Malaria-Methods for Simulation of Epidemics. A Report from Chalmers University of Technology Gothenburg.
[4]	Killeen. G.F. and Smith, T.A. (2007) Exploring the Contributions of Bed Nets, Cattle, Insecticides and Excitorepellency to Malaria Control: A Deterministic Model of Mosquito Host-Seeking Behaviour and Mortality. Transactions of the Royal Society of Tropical Medicine and Hygiene, 101, 867-880. https://doi.org/10.1016/j.trstmh.2007.04.022
[5]	Yang, H.M. (2000) Malaria Transmission Model for Different Levels of Acquired Immunity and Temperature-Dependent Parameters (Vectors). Revista de Saudepublica, 34, 223-231. https://doi.org/10.1590/S0034-89102000000300003
[6]	Ross, R. (1910) The Prevention of Malaria. J. Murray, London.
[7]	McDonald, G. (1957) The Epidemiology and Control of Malaria. Oxford University Press, London.
[8]	Birhoff, G. and Rotta, G. (1989) Ordinary Differential Equations. 4th Edition.
[9]	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
[10]	Hethcote, H.W. (2004) The Mathematics of Infectious Diseases. SIAM Review, 42, 599-653. https://doi.org/10.1137/S0036144500371907
[11]	Iyare, B.S.E., Okuonghae, D. and Osagiede, F.E.U. (2014) A Model for the Transmission Dynamics of Malaria with Infective Immigrants and Its Optimal Control Analysis. Journal of the Association of Mathematical Physics, 28, 163-176.
[12]	Diekmann, O., Heesterbeek, J.A.P. and Metz, J.A.J. (1990) On the Definition and the Computation of the Basic Reproduction Ratio R0 in Models for Infectious Diseases in Heterogeneous Populations. Journal of Mathematical Biology, 28, 365-382.
[13]	Okosun, K.O. and Makinde, O.D. (2011) Modeling the Impact of Drug Resistance in Malaria Transmission and Its Optimal Control Analysis. International Journal of the Physical Science, 28, 6479-6487.