The in-Human Host and in-Mosquito Dynamics of Malaria Parasites With Immune Responses

In this study, a mathematical model for the in-human host and in-mosquito dynamics of malaria parasite with immune responses was formualeted and analyzed. A positive invariant region of the model was established, and a basic reproduction number R0, of the model was computed. Existence and stability of two non-negative equilibrium points: malaria free equilibrium (MFE) and malaria infection equlibrium (MIE) were established. We, also proved that MFE is locally asymptotically stable if R0 < 1 and globally asymptotically stable (GAS) if R0 ≤ 1. Numerical simulations prove that MIE exists and is GAS. Moreover, our results revealed that immunity has significant influence on lowering malaria infection at blood and mosquito stages. However, an insignificant effect of immunity on both cells and parasites at liver stage infection was observed. Furthermore, the model depicts that infection decreases as lifespan of immune cells increases. The impact of immune cells to suppress production of merozoites is noted to be higher than that of antibodies to block invasion of sporozoites and merozoites.


Introduction
In human, malaria infection is initiated by a bite from a parasite-carrying mosquito, where sporozoites are injected into bloodstream and quickly migrate to the liver, where they penetrate hepatic liver cells (HLCs) and develop into schizonts that give rise to thousands of merozoites into bloodstream [7]. Merozoites invade red blood cells (RBCs), that develop into erythrocytic schizonts which finally burst and release an average of 16 merozoites [20] which either re-invade new RBCs or switch to sexual form termed gametocytes. Human-mosquito infection of malaria begins through ingestion of gametocytes into parasite-free mosquito during its blood-meal. Ingested gametocytes then fuse and develop into oocysts, where each oocyst ruptures and releases an average of 1000 sporozoites [18] that are responsible for new mosquito-human infection.
A major role of the human immune system is to defend a body against the infection-causing organisms, called pathogens. Human immune system has two main components: innate immunity and adaptive immunity. Innate immunity defends the body against any pathogenic invasion, while adaptive immunity provides protection against a specific pathogen, and usually comes into action after the infections outrun the innate immunity. Innate immune responses include macrophages, interferon and natural killer (NK) cells, while T-lymphocyte (cytotoxic T and helper T) and B-lymphocyte (B-cells) are some elements of adaptive immune responses. Cellular-mediated responses involves cell effectors such as cytotoxic T (CD8 + T) and NK cells to kill intracellular pathogens while humoral responses involve effector molecules such as antibodies (secreted by B-cells) to clear free pathogens in body's fluid such blood.
In malaria infection, both innate and adaptive immune responses are stimulated [16,15] to obstruct parasites by either preventing the re-invasion of parasite or increasing the death rate of infected cells [12,16], and both humoral and cell-mediated immune effector mechanisms are involved in immunology of malaria [15,14]. Antibodies neutralize the sporozoites and merozoites and inhibit sporozoites' invasion to HLCs [14,15] and merozoites' invasion to RBCs [20,14]. They also restrain the parasite growth [6,14]. Macrophages are activated by NK cells to intensify phagocytosis and clearance of intra-erythrocytic parasites [3]. The IFN-γ produced by CD8 + T cells (in help of CD4 + ) inhibit growth of, and kill intrahepatic parasites [14,3,15]. Moreover, antibodies and complement system that are ingested by mosquito during blood meal mediate the lysis of gametocytes and inhibit development of parasite in the mosquito [15].
In recent years, there has been increasing interest in mathematical models on intra host dynamics of malaria with immune responses [16,20,6] and references therein. Some of these studies ignored the absorption effect of merozoites into RBC [16]. However, some authors incorporate this effect since during malaria infection parasite penetrates into healthy cells. Therefore, in this case both populations (parasites and uninfected cells) decrease [2,20,6]. Moreover, the clearance of free parasites or infected cells by immune responses has been modelled either as simple mass-action [20,6] which is unbounded function or using the Michaelis-Menten-Monod function (MMMF) which is nonlinear-bounded [16]. Furthermore, [6] incorporate the effects of antibodies to inhibit parasite's growth or block invasion of host's cells by parasites using MMMF.
However, none of these studies discussed the liver stage dynamics of malaria parasites. As it has been stated in previous paragraph, some studies did not incoporate the absorption effect of parasites into uninfected cells. Moreover, some of these models ignored either the saturation effect on cell proliferation and/or suppression of parasites replication. In this study, we formulate a mathematical model for the in-human host and in mosquito dynamics of malaria with immune responses, where the MMMF is used to describe the effect of both immune responses on clearance of infected cells and free parasites. The effect of antibodies on supressing the replication of parasites at liver and blood stages of malaria infection is included. Lastly, we incoroprate the effect of antibodies picked-up by mosquito in mediating the lysis of gametocytes and preventing parasite's development in the mosquito.

Model desription
In development of this model, we extend the model by [19] by incorporating the effect of immune system. The dynamics of interactions of malaria parasites in-mosquito vector and in-human host with immune responses are described using a system of nonlinear ordinary differential equations. Variables involved in this model are: the uninfected hepatocytic liver cells (uHLCs), H; infected hepatocytic liver cells (iHLCs), I h ; hepatic schizonts, T h ; merozoites, M; uninfected red blood cells (uRBCs), R; infected red blood cells (iRBCs), I r ; erythrocytic schizonts, T r ; gametocytes, G b ; gametes, G m ; and oocysts, C. Others are sporozoites in mosquito's salivary gland, S m ; sporozoites in human, S h ; Immune cells (which includes IFN-γ, CD8 + and CD4 + T, B and NK cells, macrophages etc) against liver stage and blood stage infections, Z 1 and Z 2 respectively; and antibodies, B.
Sporozoites, S h are injected into uninfected human host at a constant rate abν, during a blood meal of infected mosquito, where a is number of mosquito bites per individual, b is number of sporozoites injected per bite and ν is probability that a mosquito bite is infective to human. Sporozoites attack the uHLCs at the rate β 1 S h H/(1 + k 1 B), where β 1 is sporozoites' infection rate to uHLCs and k 1 is efficiency of antibodies to block invasion of uHLCs.
The uHLCs are constantly recruited at rate Λ h , die naturally at rate µ h H and reduced at rate β 1 S h H/(1 + k 1 B) due to infection by sporozoites, and they are killed by immune cells (macrophages) at a rate σ sh Z 1 S h /(1 + π sh S h ), where σ sh is rate of successful removal of intra-human sporozoites by immune cells and 1/π sh is a half saturation constant of sporozoites. The iHLCs increase at a rate β 1 S h H/(1 + k 1 B) due to infection of uHLCs and die at a rate µ ih I h . They progress to schizonts at a rate α 1 I h and killed by immune cells (IFN-γ, CD8 + , NK) at a rate σ ih Z 1 I h /(1 + π ih I h ), where σ ih is rate of successful removal of iHLCs by immune cells and 1/π ih is a half saturation constant of I h . The hepatic schizonts die naturally at a rate µ th T h , rupture to release merozoites a rate δ 1 T h , and cleared by immune cells (IFN-γ, CD8 + , NK) at a rate σ th Z 1 T h /(1 + π th T h ), where σ th is rate of successful removal of hepatic schizonts by immune cells and 1/π th is a half saturation constant of hepatic schizonts.
Merozoites are released from the hepatic schizonts at a rate r 1 δ 1 T h /(1 + c 1 Z 1 ), where c 1 is efficiency of immune cells (IFN-γ and CD8 + ) to inhibit the merozoites' production. They invade uRBCs at a rate β 2 RM/(1 + k 2 B) and they die naturally at µ m M. The parameter β 2 is merozoites' infection rate to uRBCs and k 2 is efficiency of antibodies to inhibit the infection of uRBCs by merozoites. They are cleared by immune cells (macrophages activated by IFN-γ) at a rate σ m Z 2 M/(1 + π m M), where σ m is rate of successful removal of merozoites by immune cells and 1/π m is a half saturation constant of merozoites.
The uRBCs are constantly recruited at a rate Λ r from the bone marrow. Their density is reduced by natural death at a rate µ r R and due to the infection by merozoites at a rate β 2 RM/(1 + k 2 B). The iRBCs increases at a rate β 2 RM/(1 + k 2 B) and decreases due to death at a rate µ ir I r and due to progression to erythrocytic schizonts at a rate α 2 T r . The immune cells (macrophages activated by IFN-γ) phagocytize the iRBCs at a rate σ ir Z 2 I r /(1 + π ir I r ), where σ ir is rate of successful removal of iRBCs by immune cells and 1/π ir is a half saturation constant of iRBCs. The erythrocytic schizonts die at rate µ tr T r and rupture to release new merozoites at rate δ 2 T r , which they initiate a series of repetitive cycles to infect other uRBCs. Proportion p, of these newly released merozoites procceds with asexual replication cycle at a rate pr 2 δ 2 T r /(1 + c 2 Z 2 ), while the other proportion 1 − p switch to sexual form of parasites called gametocytes at a rate (1 − p)r 2 δ 2 T r /(1 + c 2 Z 2 ). The parameter c 2 is efficiency of immune cells to inhibit the production of intra-erythrocytic merozoites or gametocytes. They are also killed by immune cells (macrophages) at a rate σ tr Z 2 I r /(1 + π tr T r ). The parameter σ tr is rate of successful removal of erythrocytic schizonts by immune cells and 1/π tr is a half saturation constant of erythrocytic schizonts.
The gametocytes in blood stream increases at a rate (1 − p)r 2 δ 2 T r /(1 + c 2 Z 2 ) and decrease by natural death at a rate µ gb G b , for being ingested by mosquito at a rate qωG b and being cleared by immune cells at a rate σ gb Z 2 G b /(1 + π gb G b ). The parameters σ gb and 1/π gb are respectively, the rate of successful removal of gametocytes by immune cells and the half saturation constant of gametocytes. The gametes in mosquitoes are recruited at a rate ρqωG b /(1 + k 3 B), where ρ is number of bites a mosquito can make during its lifetime and k 3 is efficiency of antibodies picked up by mosquito during its blood meal to mediate lysis of gametocytes and prevent parasite's development in mosquito [15]. Gametes decrease by natural death at rate µ gm G m and progression to oocysts at a rate, α 3 G m . The oocysts rupture to release an avarage of r 3 sporozoites per oocyst at a rate, δ 3 C and die at a rate µ c C. The released sporozoites, S m migrate to salivary glands where they either naturally die at a rate µ sm S m or injected into a new host at rate abν. Table 1 below summarizes the variables of the model and their biological descriptions.

Model Assumptions
In development of this model, we assume the following (i) Both HLCs and RBCs are constantly recruited from bone marrow and they die naturally.

Compartmental diagram
Based on the dynamics described in Section 2.1 and the assumptions described in Section 2.2, the proposed model for the in-human host and in-mosquito dynamics of entire life cycle of malaria parasites with immune responses is shown in Figure 1, in which the variables and parameters are described in Table 1 and Table 2, respectively.

Model Equations
Based on the compartmental diagram illustrated in Figure 1 above, the in-human host and in-mosquito dynamics for the entire life cycle of malaria parasite with immune responses are governed by the following system of ordinary differential equations.
3 Analysis of the Model

Wellposedness of the model
In this section we assess the wellposedness of the model by investigating the existence and feasibility of its solution. That is, to test whether the solutions are epidemiologically (variables have biological interpretation) and mathematically (a unique bounded solution exists for all the time) wellposed. The model system (1a)-(1o) can be expressed in the compact form [11,19] dX , a 23 = ε ir Z 2 1 + π ir I r , a 24 = ε tr Z 2 1 + π tr T r , a 25 = ε gb Z 2 1 + π gb G b , a 26 = η 2 MZ 2 1 + π m M , and F is a column vector given by is Meltzer matrix since all its off diagonal elements are non negative, for x ∈ R 15 + and F ≥ 0. Therefore, the system (1a)-(1o) is positively invariant in R 15 + , meaning that an arbitrary trajectory of the sytstem started in R 15 + remains there forever. Also F is Lipschitz continous. Hence, a unique maximal solution exists and so is the feasible region for the model. Thus, the model (1a)-(1o) is epidemiologically and mathematically wellposed in the region D.

Malaria Free Equilibrium (MEF) and Basic Reproduction Number, R 0
In absence of malaria infection, all variables representing infectious classes are zeros. Thus the MFE of the system (1a)-(1o) is We calculate basic reproduction number, R 0 for the model (1a)-(1o) using the next generation matrix method [21]. In this method, R 0 is given by the spectral radius ρ(FV −1 ), of next generation matrix FV −1 where F and V are transmission and transition matrices given by and c ⃝ 2017 BISKA Bilisim Technology From (4) we obtained the inverse, V −1 , of V given by Hence, from (3) and (6), we have The basic reproduction number, R 01 , is the dorminant eigenvalue of F 1 V −1 1 . Hence using equations (7) and (8), we obtain only one nonzero eigenvalue which is A 5 . Thus, the basic reproduction number, R 0 is given by where r 0 = Λ r µ r and z 0 = Λ z 2 µ z 2 are values of R and Z 2 at malaria-free equilibrium, respectively.
From equation (9), the term is the proportion of RBCs that a merozite introduced into entirely susceptile RBC population infects before it dies (either naturally or cleared by immune cells), while the term c ⃝ 2017 BISKA Bilisim Technology represents the proportion of infected RBCs that progress to schizonts before dying and the term ( is an average duration a schizont spends before it burst or cleared by immune cells and is number of merozoites produced by a schizont when it bursts.

Local Stability of MFE
The Jacobian matrix of the system (1a)-(1o) evaluated at malaria-free equilibrium, E 0 , is Local stability of a MFE, E 0 , is determined by using the signs of real part of eigenvalues of J(E 0 ). The MFE is locally assympotically stable if and only if all eigenvalues of J(E 0 ) have negative real parts. Clearly, six eigenvalues −µ h , −µ r , −u 14 , −µ z 1 , −µ z 2 and − µ b are negative. The other nine eigenvalues are obtained from a reduced 9 × 9 submatrix, J 1 (E 0 ), given by From eighth column of J 1 (E 0 ) we observe that other eigenvalue is −u 13 which is also negative. Further reduction of this, leads us to 8 × 8 submatrix given by from which another eigenvalue −u 12 is obtained. Another reduction leads a 7 × 7 submatrix given And from J 3 (E 0 ) above we obtained another eigenvalue −u 11 , and a new reduced matrix is We investigate the signs of other remaining six eigenvalues using the trace-determinant technique. If the trace and determinat of J 4 (E 0 ) are strictly negative and positive respectively, then all eigenvalues of J 4 (E 0 ) have negative real parts. The following results were obtained using MAPLE 12, where R 0 = α 2 u 6 u 7 u 5 u 8 u 9 since the values of u ′ s given in equation (10) are all positive, then equation (12) is true only if Thus, all eigenvalues of J(E 0 ) have negative real parts only R 0 < 1. Hence, MFE is locally assymptotically stable provided R 0 < 1, which leads us to the following theorem.

Global Stability of MFE
We applied the Meltzer matrix theory to establish the global stability of MFE as used in [5], [17] and [19] by expressing the model (1a)-(1o) in the form where and where , w 8 = α 2 + µ ir + σ ir Z 2 1 + π ir I r , It can easily be observed from (15) that, all eigenvalues of A 1 are real and negative. So, the system dX n dt = A 1 (x)(X n − X E 01 ,n ) + A 12 (x)X e is globally assymptotically stable at E 01 . It can be observed from (18) and (19) that all off diagonal elements of A 2 are non-negative. Therefore, A 2 is a Metzler matrix. To establish the global asymptotic stability of MFE, we have to show that A 2 is Metzler stable matrix (all its diagonal elements are negative) by proving the following lemma as described in [11] and [13].
It is observed that equation (20) holds only when Substituting the expression of w 5 , w 6 , w 7 , w 8 and w 9 given in equation (19) evaluated at MFE, into equation (21) we get which leads us to the following theorem.

Existence and Stability of Malaria Infection Equilibrium
If R 0 > 1, then Theorem 2 suggests the existence of malaria-infection equilibrium (MIE), which is given by Due to the complexity of the model system (1a)-(1o), it is found awkward to express MIE explicitly. Therefore, the existence and stability are established numerically in the next section.

Numerical Simulations and Discussion
In this section, we perform some numerical simulations of the model (1a)-(1o), to illustrate the dynamics of model using MATLAB symbolic package. In the simulation of this model, initial values are almost assummed to allow computer executions, and their values are as listed in Table 3 c ⃝ 2017 BISKA Bilisim Technology Initial values 3000 0 0 2000 500000 0 1000 3000 1500 1000 2000 2000 10 10 0 The numerical values of parameters used for simulation of our model are listed in Table 2. These values are either assumed or taken from some related studies among existing literature. The reason as to why some parameters values are assumed is that, mathematical modelling on liver and/or mosquito stage dynamics of malaria infection have not yet been done or the values found on existing literatures are not suitable for this model. However, our main concern is not on accuracy of these parameter values, rather is on the effect of these parameters on basic reproduction number, which alerts on how and where to target to eradicate or control the disease [6].
We performed numerical simulations to establish the existence of malaria infection equilibrium (MIE) as stated earlier. Figures 2a , 2b, 2c and 2d indicate that each variable varies with time and reaches a constant value (i.e., a value at MIE). Therefore, Figure 2  for stability of E * . Figure 3 depicts that with different initial values, each variable in exo-erythrocytic cycle converges to certain value (value at MIE). Figure 4 shows that each variable of erythrocytic cycle converges to a certain steady value irrespective of initial value it takes. The case is similar for variables of sporogonic cycle and for immune responses as indicated in Figure 5 and Figure  6 respectively. Therefore, with reference to Figures 3, 4, 5 and 6, we conclude that, malaria infection equilibrium, E * is, globally asymptotically stable. However, the immune responses have shown to have little or almost no effect on attacking liver stage malaria infection. This is because there is unnoticable change in number of uHLCs and iHLCs even after immune system being included the model as indicated in Figure 8. It can be observed that despite of very big change on value of k 1 (efficiency of antibodies to block invasion of hepatocytes by sporozoites) from 0.075/day to 0.9/day, but still the change on number of uHLCs and iHLCs is almost negligible. The case is the same to sporozoites injected and hepatic schizonts. These results suggest that at the liver stage is not a good target for intervetion using immune responses.
This finding supports the [15]'s argument that the function of antibodies to sporozoites is thought to be insignificant. Furthermore, the immune responses has shown positive results on the sporogonic stages of malaria parasites as shown in Figure 9. This implies that, the antibodies taken up by mosquito during the blood meal results on lysis of gametocytes and prevent the development of parasites within mosquito. Consequently, it can lead us to a population of non-infectious mosquitoes, and reduce further mosquito-human malaria transmission.
Since  influence on the number of uninfected RBCs. It shows that the number of uRBCs increases with immune cells' lifespan. Hence, using Figure 10 we argue that introducing a long-term immunity may significantly reduce the infection.
With reference to basic reproduction number, R 0 ; as expressed in equation (9) it is noted that R 0 can be made less than or equal to unity by: decreasing the infection rate, β 2 of uRBCs by merozoites; increasing death rates of iRBCs, blood-schizonts and merozoites; increasing the rates at which iRBCs, blood-schizonts and merozoites are cleared by immune cells; and increasing the rate at which immune cells suppress the production of merozoites from blood-schizonts. Therefore, any biological means that can be implemented to facilitate these may have impact on development of control strategies.

Conclusion
In this work, we have formulated and analyzed a mathematical model for the in-human host and in-human dynamics of malaria parasite with immune responses. In this model, we include the effect of immune responses to block invasion of sporozoites and merozoites on hepatic liver cells and red blood cells respectively. The effect of immune responses to inhibit the production of merozoites from both liver and blood cells was included. Additionally, the model incudes terms    for influence of immune responses on clearance of both parasites and infected cells. Finally, we include the term for effect of antibodies on gametocytes picked-up during the blood meal. In all cases immune response are described using the nonlinear-bounded Michaelis-Menten-Monod function.
A positive invariant region, where the model is epidemiologically (variables biological interpretation meaningful) and mathematically (always a unique bounded solution exists) well-posed was established. Using the next generation method, basic reproduction number R 0 , of the model was computed. Also, existence and stability of two non-negative equilibrium points: malaria free equilibrium (MFE) and malaria infection equlibrium (MIE) were established. Furthermore, we proved that MFE is locally asymptotically stable if R 0 < 1 and globally asymptotically stable (GAS) if R 0 ≤ 1. In addition to that, we noted that the impact of immune cells to suppress production of merozoites is higher than that of antibodies to block invasion of sporozoites and merozoites into liver and blood cells respectively. This is because none of k ′ s, efficiency of antibodies to inhibit invasion appeared into the expression of R 0 , while one of the c ′ s (c 2 ), efficiency of immune cells to suppress production of merozoites does appear.
Numerical simulations prove the existence of MIE, which is GAS irrespective of the initial values do state variables have. Moreover, in comparison with the results of [19], our results revealed that including immunity has significant c ⃝ 2017 BISKA Bilisim Technology  influence on lowering infection at blood and mosquito stages. An increase on number of uninfected cells, and a decrease on number of infected cells and free parasite were noted. Also, antibodies picked-up by mosquito seems to have an effect on reducing the number of parasites within mosquito.
Based on the above results, this study suggests that the a combined drug therapy should be boosted so as to improve their ability to suppress parasite's production in bloodstream. Also, causing lysis to gametocytes is of great importance in preventing parasite' development in mosquito as it may lead us to a population of non-infectious mosquito. Hence, reduce mosquito-human infection.
Our work is an example on how mathematical models can be used to get insight of complex systems like malaria life cycle. Therefore, it provides basis for mathematical models to study the in-host and in-mosquito dynamics malaria and immunity. In future work we propose the extension this model by incorporating the effect of antimalaria drugs. c ⃝ 2017 BISKA Bilisim Technology