SEIRS Mathematical Model for Malaria with Treatment

In this paper a deterministic mathematical model for the spread of malaria in human and mosquito populations are presented. The model has a set of eight non – linear differential equations with five state variables for human and three for mosquito populations respectively. Susceptible humans can be infected when they are bitten by an infectious mosquito. They then progress through the exposed, infectious, treatment and recovered or immune classes before coming back to the susceptible class. Susceptible mosquitoes can become infected when they bite infectious humans, and once infected they move through exposed and infectious class. However, mosquitoes once infected will never recover from the disease during their lifetime. That is, infected mosquitoes will remain infectious until they die. Formula for the basic reproduction number R0 is established and used to determine whether the disease dies out or persists in the populations. It is shown that the disease – free equilibrium point is locally asymptotically stable using the magnitude of Eigen value and Routh – Hurwitz stability Criterion. Result and detailed discussion of the analysis as well as the simulation study is incorporated in the text of the paper lucidly.


Introduction
Malaria is a disease having huge social, economic and health burden on human beings since ancient times. It is predominantly present in the tropical countries where rains and humid are more. Though the disease has been identified hundreds of years back and many prevention and treatment methods are invented it still remains as a major public health problem.
The WHO estimated that there were 219 millions of cases of malaria world -wide resulting in 0.66 million deaths in the year 2010 [11,13]. However, others have estimated that the number of cases of falciparum malaria alone lies between 350 and 550 millions and deaths rose from 1.0 million in 1990 to 1.24 millions in 2010. Further, these estimations argue that the WHO's statistical figures are much smaller than the actual [8,9,12].
Majority of malaria cases accounting about 65 percent occur in children of below 15 years old [9]. About125 million pregnant women are at risk of the infection in Sub -Saharan Africa and the maternal malaria is leading to an estimated 0.2 million infant deaths each year [7]. There are about 10,000 malaria cases per year in Western Europe and 1300 -1500 in the United States [13]. About 900 people died of the disease in Europe during a decade from 1993 to 2003 [7].
Both global incidences of the disease and resulting mortalities have been observed declined in recent years. According to world health organization WHO deaths attributable to malaria in 2010 were reduced by over one third from 2000 having 985,000 incidents. The reduction in mortality is largely resulted due to the widespread use of insecticide -treated nets and artemisinin -based combination ACT therapies [6].
Mathematical models are particularly helpful as they consider and include the relative effects of various sociological, biological and environmental factors on the spread of a disease. The models have been playing an important role in the development of various controlling techniques for malaria epidemic. Analysis of mathematical models is important because they help in understanding the spreads of malaria so that suitable control techniques can be adopted for controlling the disease.
In this study, treatment compartment has been included in the basic SEIRS epidemic model and extended it to SEITRS model. The inclusion of this compartment and extension of the model is justified since treatment plays a crucial role in controlling the spread of malaria. The existing SEIRS mathematical model has four classes or compartments viz., Susceptible, Exposed, Infected and Recovered. The inclusion of the Treatment class to SEIRS model extended it to SEITRS model. Treatment represented by the letter T is a technique used to control the spread of malaria disease. The present modified model is an extension or modification of the existing mathematical models used to deal with malaria epidemic viz., SIR, SEIR, SEIS, SEIIR [1,3,4,10,15]. Thus here in this study, the SEITRS model is presented and shown that it describes the dynamics and controlling mechanism of malaria disease. The effect of treatment technique in the spread of malaria is also analyzed. Simulation studies of the model with variable values of sensitive parameters of the spread of malaria are performed and the results are incorporated in the text of the paper. Important conclusions are drawn and necessary recommendations have been made.

The SEITRS Model
Malaria disease occurs due to the bites of infected female anopheles mosquitoes. That is, the rate of malaria disease can be decreased by accelerating mortality or death rate as well as decelerating infection rate of the female Anopheles mosquitoes.
Hence, it is needed to invent intervention strategies to increase mortality rate of the female anopheles mosquitoes in addition to the natural mortality and to decrease infection rate of such mosquitoes. However, the authors strongly believe that treatment of humans is much preferable than killing mosquitoes. Killing is offensive but treatment is defensive. Defense i.e., treatment strategy gives better results. Hence, in this study more attention is given to treatment.
Many strategies have been designed and implemented to increase the mortality rate. These intervention strategies can increase the mortality rate of female anopheles mosquitoes and help directly to control the birth and spread of malaria. In the recent times many such intervention strategies are adopted worldwide. The important intervention strategies include (i) indoor residual spraying, IRS and (ii) insecticide treated bed net, ITN.
In the present study treatment of infected human population is considered so that the female anopheles mosquitoes will not get infection and will not pass the disease to susceptible humans. Further, it is shown that the intervention strategy treatment plays a significant role in controlling the spread of malaria disease.

Classification of Human Compartments
The present mathematical model is built with an objective of analyzing the effect of intervention strategies such as treatment on the spread of malaria disease. Thus, the treatment class is considered and added to SEIRS model so as to extend it to SEITRS model. In the SEITRS model whole the human population is divided into five classes or compartments based on the properties they exhibit. The description of the classes and their properties are briefly explained here under.
Susceptible human class : It contains humans those do not have malaria disease but are likely to be bitten by infected female anopheles mosquitoes causing malaria. The people in this class do not take any protective measurements against malaria disease. In general, all the common people are included in this class. Humans of the susceptible class have the potential to remain either in the same class or to enter into exposed class.
Exposed human class : It contains humans that have got infection but are not infectious. That is, these humans exhibit the symptoms and suffer from the disease but are not capable of propagating the disease to female anopheles mosquitoes. People from susceptible class come into exposed class.
Infected human class : It contains humans those are already infected and got malaria disease. These people while suffering from the disease are capable of transferring the disease to female anopheles mosquitoes. People come into infected class from exposed class. People from infected class have got the potential to enter into recovered class directly if the disease is naturally cured or enter into treatment class for treatment.
Treatment human class : It contains humans those are already infected and got malaria disease. The people of this class are given medical care or clinical treatment. People into the treatment class come from only infected class. After completion of the treatment successfully, the people from treatment class will go to recovered class.
Recovered human class : It contains people who recover from the malaria disease either by treatment or by natural reasons and return to normal status of health. People come into the recovered class from both infected and treatment classes. After recovery the recovered people will go again into susceptible class.
Further, it is reasonable to assume that sick people cannot travel far away or they cannot move from place to place. Thus, the immigrations of malaria infected people were not included in this model.

Classification of Mosquito Compartments
Anopheles is a kind of mosquito that has the capacity of transmitting malaria parasite to humans. However, Anopheles male mosquitoes cannot carry plasmodium parasite and hence they do not cause for malaria in humans. Thus, anopheles male mosquitoes are not included in the model. In contraction to males, infected female anopheles mosquitoes through bites can cause malaria in humans. Thus, only female anopheles mosquito population is considered and included in the present malaria model.
Further, it is to be noted that only female anopheles mosquitoes bite humans for blood meals but not male ones. In the present model all the female anopheles mosquitoes are divided into three classes or compartments based on the properties they exhibit. The three compartments of such mosquitoes and the respective properties are described briefly here under.
Susceptible mosquito class : It contains female anopheles uninfected mosquitoes. The mosquitoes of this class are not infected so far and thus their bites will not lead to malaria disease in humans. The population size of these mosquitoes grows with the constant recruitment rate of and diminishes with the natural death rate of α . The susceptible female anopheles mosquitoes get infected on feeding themselves with the blood of malaria infected humans and then they leave susceptible mosquito class and enter infected mosquito class. Susceptible mosquitoes bite the infected humans with the rate of . Further, it is assumed that the susceptible mosquitoes get infected with a rate of when they bite infected humans.

Exposed mosquito class
: It contains exposed female anopheles mosquitoes. Though these mosquitoes are exposed to the disease they are not infectious i.e., these are not capable of transferring the disease to susceptible humans. The exposed mosquitoes will reduce with a natural death rate ofα .
Infected mosquito class : It contains infected female anopheles mosquitoes. The mosquitoes of this class are already infected and their bites can lead to malaria disease in humans. The population size of this class grows with a transfer rate proportional to from the susceptible mosquito class. The infected mosquitoes diminish with the natural death rate of α . The infected female anopheles mosquitoes transfer infection to susceptible humans when the former feed them self with the blood of the latter. Infected mosquitoes bite the susceptible humans with the rate of . Further, it is assumed that the infected mosquitoes transfer infection with a rate of to humans when they bite susceptible humans.

Formulation of the Model
The mathematical model of malaria transmission and dynamics considered in this study consists of SEITRS compartments for human population and SEI for mosquito population. The model is formulated for the spread of malaria in the human and mosquito populations with the total population sizes at time denoted by and respectively.
The human populations are further compartmentalized into epidemiological classes as susceptible , exposed , infected , treated and recovered . The mosquito populations are similarly compartmentalized into epidemiological classes as susceptible , exposed and infected . The vector component in the present model does not include any immune class since infected mosquitoes never recover from the infection. That is, their infective period ends with their death due to their relatively short lifecycle. Hence, the total population sizes of both humans and mosquitoes are the sum of the sizes of the respective classes: = + + + + and = + + . The immune class in the mosquito population is negligible due to the reason already stated and death occurs equally in all classes. However, this model can also be used for diseases that persist in a population for a long period of time with vital dynamics.
The present model is constructed incorporating a set of basic assumptions: (i) Total population sizes of both human and vectors are assumed to be constant (ii) The recovered individuals of human population do not develop permanent immunity (iii) The population sizes of both human and vector compartments are non -negative. This fact is proved and presented in the text in form of theorem 1 (iv) Also, all the parameters involved in the model equations are nonnegative quantities. Description of these parameters is given in Table 1 (v) All newborns of both human and mosquitoes are susceptible to infection (vi) the development of malaria starts when the infectious female mosquito bites susceptible human host or susceptible female mosquito bites infectious human host (vii) individuals move from one class to another as per their status with respect to the evolution of the disease (viii) Humans enter the susceptible class through the constant recruitment rate , leave from the susceptible class through death rate and exposed rate (ix) susceptible humans enter the exposed class through an exposed rate of and leaves this class with an infection rate of and (x) all human individuals whatever their status is subjected to a natural death occurring at a rate of and disease induced death rate of . Hence the flow chart of the model is presented here under with the above assumptions.

Parameters
Interpretations Recruitment rate of susceptible humans Force of infection of humans from susceptible state to exposed state α Natural death rate of humans Rate of progression of humans from exposed to infectious state Rate of progression of humans from infectious to treatment state Disease -Induced death rate for humans Flow rate of humans from infectious to recovered state due to treatment Flow rate of humans from infectious to recovered state due to natural reasons Flow rate from recovered to susceptible class due to loss of immunity for humans Recruitment rate of susceptible mosquitoes Flow rate of mosquitoes from susceptible to exposed state : Natural death rate of mosquitoes Flow rate of mosquitoes from exposed to infectious state

Variables
Interpretations Population size of susceptible human class E Population size of exposed human class I Population size of infectious human class T Population size of treatment human class R Population size of recovered or immune human class S Population size of susceptible mosquitoes class E Population size of exposed mosquitoes class I Population size of infectious mosquitoes class Susceptible humans get infected at a certain probability when they contact infectious mosquitoes. They then progress through the exposed, infected, treated and recovered classes before reentering the susceptible class. Susceptible mosquitoes get infected at a certain probability when they bite infectious humans and then move through the exposed and infectious classes. Both humans and mosquito populations grow naturally following logistic function.
Applying the assumptions, definitions of state variables and parameters as stated above the coupled system of non -linear differential equations describing the dynamics of malaria in the human and mosquito populations are formulated as Further, the initial population sizes are considered as initial conditions of the system of model equations. They are S 0 S % , E 0 E % , I 0 I % , T 0 T % , R 0 R % and S 0 S % , E 0 E % , I 0 I % . Thus, the initial conditions representing population sizes are non -negative. That is, the initial conditions are non -negatives such that S 0 ' 0, 0 ' 0, 0 ' 0 0 ' 0, 0 ' 0, S 0 ' 0, 0 ' 0, I 0 ' 0 . Also, the total population sizes and ( can be determined by

Mathematical Analysis of the Model Equations
The mathematical analysis of the present model described by the system of differential equations (1) is conducted and the results are presented here. The model equations are required be analyzed in a feasible region. Hence, it is mandatory to identify in the first step the region. All state variables and parameters are supposed to be positive since they represent populations and their growths. The invariant region for the model formulated can be found to be Therefore, any solution of the system of ordinary differential equations (1) is feasible for all t > 0if that enters the invariant region Ω = Ω × Ω .

Positivity of the Solution
In order that the model equations (1) are biologically and epidemiologically meaningful and well posed it is appropriate to show that all the state variables are nonnegative. This requirement is stated as a theorem and its proof is provided as follows: Theorem 1: Any solution of the model equations (1) together with the initial conditions is non -negative. Alternatively it can be stated as follows: if S 0 > 0, 0 > 0, I 0 > 0, 0 > 0, R 0 > 0, S 0 > 0, 0 > 0 and I 0 > 0, then the solution region 4S t , t , I t , t , R t , S t , t , I t 7 of the system of equations (1) is always non -negative. Proof: In order to show that every solution of the dynamical system (1) is positive each differential equation of the system is considered separately and treated.
Positivity of infective mosquito population I : Consider the last differential equation of the system of differential equations (1) as dI 9 ⁄ = − . Note that the term is a positive quantity and thus by dropping it the foregoing equation can be expressed as an inequality as dI 9 ⁄ ≥ − . Further, separation of variables reduces it to dI ⁄ ≥ − dt and finally integration yields a solution as ≥ ; % e = > ? @ AB C D E . Since the initial population size of infected mosquitoes is non -negative i.e., % ≥ 0and the exponential function is always positive it can be concluded that ≥ 0. Thus, the solution of t is always non -negative. Positivity of exposed mosquito population E : Consider the differential equation dE dt ⁄ = − − of the system of differential equations (1). Note that the term is a positive quantity and thus by dropping it the foregoing equation can be expressed as an inequality as Since the initial population size of exposed humans is nonnegative i.e., E % ≥ 0 and the exponential function is always positive it can be concluded thatE ≥ 0. Thus, the solution of E t is always non -negative. Positivity of susceptible human population S : Consider the first differential equation 9 9 ⁄ = + − − α S of the system of differential equations (1). Since the term + is a positive quantity thus by dropping it the foregoing equation can be expressed as an inequality as 9 9 ⁄ ≥ − + α S . Further, separation of variables reduces it to 9 S ⁄ ≥ − + α 9 and finally integration yields a solution as S ≥ ;S % e = > M H N @ 5α Q AB C D E . Since the initial population size of exposed humansis non -negative i.e., S % ≥ 0 and the exponential function is always positive it can be concluded that S ≥ 0. Thus, the solution of S t is always non -negative.

Boundedness of the Solution Region
Theorem 2: The non -negative solutions of the model system (1) are bounded. That is, the population sizes of both humans and mosquitoes are bounded.
Proof: Note that if the total population sizes of both human and mosquitoes are bounded then the sizes of each compartment are also bounded. Hence, it is sufficient to prove that the total sizes of both the populations are bounded. That is, the solutions lie in the bounded region.
Boundedness of human population ≤ α ⁄ : The rate of change of total human population size = + E t + + T t + can be obtained as 9 9 ⁄ = 9 9 ⁄ + 9 9 ⁄ + 9 9 ⁄ + 9 9 ⁄ + 9 9 ⁄ . On substituting the model equations (1) and after some algebraic simplifications the fore going differential equation takes the simplified form as 9 9 ⁄ = − − α N . Now, since the term appearing on the right hand side is a positive quantity, without loss of generality it can be expressed as an inequality as 9 9 ⁄ ≤ − α N . Thus, the solution of this differential inequality is found to be ≤ α ⁄ + S % − α ⁄ TF =U Q V giving that ≤ α ⁄ as → ∞. The term % denotes the initial total human population. It can be interpreted that the total human population grows and asymptotically converges to a positive quantity α ⁄ . Thus it is an upper bound of the total human population . Whenever the initial human population starts off either lower below or higher above α ⁄ then it grows or decays over time and finally reaches the upper asymptotic value.
Boundedness of mosquito populationN t = ⁄ : Similarly, the rate of change of total mosquito population size = + E t + is obtained on differentiating as 9 9 ⁄ = 9 9 ⁄ + 9 9 ⁄ + 9 9 ⁄ . On substituting the model equations (1) and after some algebraic simplifications the foregoing differential equation takes the simplified form 9 9 ⁄ = − . Also, solution of this differential equation is found to be = + S % − F =? @ Y ⁄ T ⁄ giving that N t → ⁄ as → ∞. The term N % denotes the initial total mosquito population. It can be interpreted that the total mosquito population grows and asymptotically converges to the positive quantity ⁄ . Thus ⁄ is an upper bound of the total mosquito population N t .

Disease Free Equilibrium
Disease -free equilibrium points are steady state solutions of the model equations. Further, at these points malaria does not present in the human population and similarly plasmodium parasite does not present in the mosquito population.
Clearly equilibrium points are the solutions of the model equations satisfying the conditions 9 9 ⁄ = 9 9 ⁄ = 9 9 ⁄ = 9 9 ⁄ = 9 9 ⁄ = 9 9 ⁄ = 9 9 ⁄ = 9 9 ⁄ = 0. Making uses of these conditions in the system of non -linear differential equations of model (1) it can be obtained as Further, disease -free implies that = 0 and = 0 since they are diseased classes of human and mosquito populations respectively. Using disease -free conditions the above set of equations reduce to

Basic Reproduction Number
The basic reproduction number R % is the number of secondary infections that one infectious individual would create over the duration of infectious period provided that everyone else is susceptible. To find formula for it the next generation operator approach can be used as described by Diekmann et al., (1990). Reproduction number is the threshold for many epidemiology models which determines whether a disease can invade in a population or not.
If % < 1 then each infected individual produces on average less than one new infected individual so it is expected that the disease would die out. On the other hand if % > 1 then each individual produces more than one new infected individual so it is expected that the disease would continue spreading in the population. This means that the threshold quantity for eradicating the disease is to reduce the value of % to be less than one.
To determine the basic reproduction number % using the next generation approach for the present model the following steps are to be followed: The reproduction number % is defined as the largest eigenvalue of the next generation matrix. The formulation of this matrix involves determining two classes viz., infected and non-infected from the model. That is, the basic reproduction number cannot be determined from the structure of the mathematical model alone but depends on the definition of infected and uninfected compartments.
Assuming that there are \ compartments in the model of which the first ] compartments are of infected individuals.^_ ` = ^_ = ` − ^_ 5 ` Where ^_ 5 ` is the rate of transfer of individuals into compartment aby all other means and ^_ = ` is the rate of transfer of individual out of the a Y compartment. It is assumed that each function is continuously differentiable at least twice in each variable. The disease transmission model consists of nonnegative initial conditions together with the following system of equations: `b _ = ℎ _ ` = c _ ` − ^_ ` , a = 1,2,3, … \ where xb is the rate of change of `.
The next is the computation of the square matrices c and ^ of order ] × ], where ] is the number of infected classes, defined by c = hic _ i`j ⁄ `% k and ^ = hi^_ i`j ⁄ `% k with 1 ≤ a, l ≤ ], such that c is nonnegative, ^ is a non singular matrix and `% is the disease -free equilibrium point (DFE). Since c is nonnegative and ^is nonsingular, then ^= m is nonnegative and also c^= m is nonnegative. Hence the matrix of c^= m is called the next generation matrix for the model. Finally the basic reproduction number % is given by % = n c^= m (5) Here in (5), n o denotes the spectral radius of matrix o and the spectral radius is the biggest non -negative eigenvalue of the next generation matrix.
Rewriting the model equations (1) Now, with the partial derivatives of (6) with respect to , , , the Jacobian matrix of c _ is constructed and that at the disease -free equilibrium point (4) takes that form as Similarly, with the partial derivatives of (7) with respect to , , , the Jacobian matrix of ^_ is constructed and that at the disease -free equilibrium point (4) reduces to the form as In order to find inverse of the matrix ^ it is needed to find ~9l ^ From (12), it is now possible to calculate the eigenvalues to determine the basic reproduction number % by taking the spectral radius of the matrix c^= m .
Thus, the characteristic equation is computed by|c^= m % − | = 0, and it yields Therefore the eigenvalues are given by From the formula for % , it can be quantified that higher values of , , , and can allow the outbreak of the disease.
Conversely, for smaller values of , , , and the disease dies out. The reproduction number is a powerful parameter which measures the existence and stability of the disease in the human and mosquito population.

Theorem 3:
The DFE E % of the system (1) is locally asymptotically stable if % < 1 and unstable if % > 1. Proof: To construct Jacobian matrix of the model, consider that the right hand sides of the model equations (1) are represented by the following functions: Here • _ = • _ , E , , T , , , E , , ∀ a = 1, 2, … , 8. The Jacobian matrix is given by In •, the matrix elements are partial derivatives denoted by the notations• m¡ H = i• m i ⁄ and so on. Further, the Jacobian matrix of model (1) at the disease -free equilibrium % reduces to the form as The characteristic equation |J E % − | = 0of the matrix (14) is given by It is simple to observe that all the four eigenvalues m , " , 6 , -are negative quantities.
Therefore, it can be concluded that the disease freeequilibrium E % is locally asymptotically stable if % < 1 and unstable if % > 1.

Numerical Simulations
To investigate the role played by epidemiological parameters in the model dynamics, the numerical simulations are carried out using DEDiscover 2.6.4 software. The parameters used, their estimated values and appropriate sources are given in tables 3 and 4. Different values i.e., % < 1 and % > 1 are obtained for the reproduction number % using the parametric values given in the former and the latter tables. Numerical simulation results are given in the following figures together with detailed descriptions.
Using parametric values given in Table 3 together with the initial conditions % = 120, % = 60, % = 20, % = 10, % = 18, % = 160, % = 80 and % = 150 a simulation study is conducted and the results are given in Figure 2. Further, using these parametric values the reproduction number is computed and found that its value to be % = 0.00251 . Here it can be observed that % < 1 showing that the disease free eqilibrium point is stable.   Table 3.
In Figure  2, the  fractions  of  the  populations , , , , , , and are plotted versus time. The susceptible populations will initially decreases with time and then increases and the fractions of infected human populations decrease. More over because of treatment the recovered human population tremendously increases and also the susceptible and infected mosquito population decreases over time.  shows the dynamics with time of susceptible, exposed, infectious, treatment and recovered human populations.  [12] Assuming the parameter values of Table 4 together with the initial conditions % 120, % 60, % 20, % 10, % 18, % 160, % 80 and % 150 a simulation study is conducted and the results are as shown in Figure 2.
By using the parameter values of Table 4, the reproduction number is calculated and evaluated to be % 51.70145.
Here it can be observed that % ' 1showing that the disease free eqilibrium point is unstable and thus the disease persists in the population. In Figure 4, the fractions of the populations , , , , , , and are plotted versus time. It is shown that the treatment rate of the infected human population decreases the infected human populations over time.

Conclusion
In this paper, a model for malaria is formulated taking into account the compartmentalization of both the human and mosquito populations. An SEITRS model with the inclusion of treatment for infected human is formulated for humans and similarly an SEI model is formulated for mosquitoes. Constant recruitment rates for both human and mosquito populations are considered. Dynamics of mosquito population is studied along with that of human populations since the mosquito population determines to a large extent whether a malaria outbreak will occur or not.
Further, it is shown that the solution of the model is both positive and bounded. Hence, it is interpreted that the model equations are mathematically and epidemiologically well posed. The disease free equilibrium theory is applied to study the stability analysis. In particular, the stability properties were investigated by paying more attention to the basic reproduction number.
In this study the treatment to infected humans is incorporated in already existing SEIRS model and developed an improved SEITRS model. This improvement remains as a reasonable contribution for controlling malaria disease. From the numerical results, it is observed that treatment to infected human population has a strong control over the spread of malaria disease.
Moreover, the dynamics of an SEITRS model is studied and applied to malaria transmission between human and mosquito populations. The basic reproduction number is derived. Further it is shown that a disease -free equilibrium exists and its stability condition % Z 1 is proved.
The analysis shows that if the reproduction number is less than one unit then the DFE is stable. This implies that only susceptible human population is present and all the other human populations reduce to zero. Thus, it can be concluded that the disease dies out as shown in Figure 2.
On the other hand, if the reproduction number is greater than one unit then the DFE is unstable. This fact has been verified by numerical simulation and the result is illustrated in Figure 4.
Clearly, the numerical simulation study supports that the DFE is locally asymptotically stable whenever the reproduction number is less than one unit. It is also noticed that in order to reduce the basic reproduction number below one, it is very necessary to give a focus on the reduction of the rate of infected human population through treatment. Hence, for reduction of malaria infection it is recommended that humans must be tested for the disease and provide sufficient treatment.