Global Analysis of a Reaction-Diffusion Within-Host Malaria Infection Model with Adaptive Immune Response

Malaria is one of the most dangerous global diseases. This paper studies a reaction-diffusion model for the within-host dynamics of malaria infection with both antibody and cell-mediated immune responses. The model explores the interactions between uninfected red blood cells (erythrocytes), three types of infected red blood cells, free merozoites, CTLs and antibodies. It contains some parameters to measure the effect of antimalarial drugs and isoleucine starvation on the blood cycle of malaria infection. The basic properties of the model are discussed. All possible equilibrium points and the threshold conditions required for their existence are addressed. The global stability of all equilibria are proved by selecting suitable Lyapunov functionals and using LaSalle’s invariance principle. The characteristic equations are used to study the local instability conditions of the equilibria. Some numerical simulations are conducted to support the theoretical results. The results indicate that antimalarial drugs with high efficacy can clear the infection and take the system towards the disease-free state. Increasing the efficacy of isoleucine starvation has a similar effect as antimalarial drugs and can eliminate the disease. The presence of immune responses with low efficacy of treatments does not provide a complete protection against the disease. However, the immune responses reduce the concentrations of all types of infected cells and limit the production of malaria parasites.


Introduction
Malaria is an infectious diseases caused by Plasmodium parasite. In 2018, about 228 million cases of malaria occurred around the world [1]. Most of these cases were in the World Health Organization (WHO) African Region with 93%, followed by the South-East Asia Region with 3.4% and the Eastern Mediterranean Region with 2.1% [1]. The global malaria deaths in 2018 were estimated at 405,000 deaths, where children under age 5 years accounted for 67% of these deaths. In 2018, there were 11 million pregnant women exposed to malaria infections and they delivered about 872,000 children with low birth weight [1]. Plasmodium falciparum (P. falciparum) is the most popular malaria parasite and it is accounted for the highest number of malaria cases [1]. The other four Plasmodium species that infect humans are P. vivax, P. ovale, P. malariae and P. knowlesi [2,3]. The malaria infection within a host is composed of two phases, the liver stage and the blood stage [3]. The blood stage is associated with the clinical symptoms of the infection. The infection begins by the bite of an infected Anopheles mosquito to take a blood meal. During the bite, it releases the malaria parasite in the form of sporozoites into the bloodstream of the host. The sporozoites migrate to the liver where they infect the liver cells and develop into schizonts. After rupturing, the schizonts release the parasites in the form of merozoites into the bloodstream where the blood stage starts [2,4,5]. The merozoites invade red blood cells and develop within them through three morphological stages. The youngest stage is the ring stage which grows into the trophozoite stage, and finally into the schizont stage. The schizont ruptures to produce 8-32 daughter merozoites which invade healthy red blood cells and repeat the infection cycle [3]. The merozoites are the asexual form of the malaria parasite. Nevertheless, some of the merozoites develop into gametocytes which are the sexual form of the parasite [5]. Artemisinin-based combination therapies (ACTs) are the current standard treatment for the blood stage of P. falciparum malaria infection [6,7]. In this paper, we focus on the asexual blood-stage dynamics of P. falciparum.
Mathematical modeling has been used to understand the dynamics of malaria parasites during infection and to design more effective antimalarial drugs. Malaria models were built using similar ideas to those used in virus dynamics models [8][9][10][11][12][13][14][15][16]. The basic model for the within-host dynamics of malaria infection was proposed by Anderson et al. [17]. The model takes the form where g = 1. The variables U(t), Y(t), and M(t) denote the concentrations of uninfected red blood cells, infected red blood cells, and free merozoites, respectively. Uninfected red blood cells are recruited from the bone marrow at a constant rate µ, infected by free merozoites at rate βU M, and die at a natural death ratedU. Infected red blood cells die at rate eY and rupture to produce η merozoites per infected cell. The total production rate of merozoites is given by ηeY. Free merozoites either die at rate f M or infect red blood cells at rate βU M [17]. A similar model was investigated by Anderson [18] but with g = 0. The models in [17,18] were extended in multiple ways to include other components or processes. In [19], A. Saul argued that the basic models contain an error which leads to overestimation of parasite growth rates. To avoid this error, he suggested using the growth rate constant e(ln(η) + 1) instead of eη. However, Gravenor and Lloyd [20] mentioned that Saul's suggestion does not solve the problem of large growth rates. Alternatively, they proposed a model with n categories for the infected erythrocytes. The transition between these classes is governed by constant rates. Hoshen et al. [21] introduced a delay differential equation model of blood-stage malaria infection. The model contains a parameter τ to model the lifespan of the infected erythrocytes (τ = 48 h in the case of P. falciparum). They converted the nonlinear equations into linear equations and they got the analytic solution by using suitable biological analysis of non-severe malaria. Iggidr et al. [22] performed a global analysis of the model suggested in [20] with a general recruitment rate for the uninfected red blood cells. In addition, they proposed and analyzed a model with n strains and k classes of infected erythrocytes. They identified the conditions under which one strain persists while all other strains die out. Saralamba et al. [23] developed a model to explore parasite population dynamics in age and time. They examined the effect of antimalarial drug on the killing rate of parasites at the three sequential asexual blood stages; rings, trophozoites, and schizonts. They suggested that the drug resistance is concentrated at the ring stage. Demasse and Ducrot [24] extended the model developed in [22] and proposed an age-structured blood-stage malaria model. They showed the global stability of the disease-free equilibrium and the positive equilibrium associated with the strongest strain. The blood-stage infection stimulates different immune responses to fight the infection and limit the growth of parasites [25]. The adaptive immunity plays a major role in fighting the disease. There are two types of adaptive immunity which are antibody immunity and cellular immunity. Cellular immunity stimulates the production of cytotoxic T lymphocytes (CTLs) to kill infected cells. Antibody immunity is mediated by B cells which produce antibodies to clear free merozoites from the blood [25,26]. It was mentioned that antibody immunity is more effective than cellular immunity [26]. However, this depends on the magnitude of each one of the two immune responses and their ability to constrict parasite replication. For simplicity, many works have added only one component to their models to represent the total capacity of the immune response directed against both infected red blood cells and free merozoites. For example, Anderson et al. [17] added the immune response to model (1) as follows where Z(t) denotes the concentration of immune cells. The clearance rates of infected red blood cells and free merozoites by immune cells are given byα 1 YZ andα 2 MZ, respectively. The infected red blood cells and free merozoites stimulate the production of immune cells at ratesρ 1 YZ andρ 2 MZ, respectively. The death rate of immune cells is given by hZ. The other parameters have the same biological meanings as those in model (1). Hetzel and Anderson [27] studied an extended version of model (1) with an additional equation for an immune response. They performed experiments to estimate the invasion rate of red blood cells. Tumwiine et al. [26] explored an extended version of model (2), where they established the existence and global stability of equilibrium points. Niger and Gumel [4] extended the models in [20] and [22] by adding the effects of immune response and malaria vaccines. They proved the global stability of the parasite-free equilibrium. Li et al. [25] expanded the model studied by Anderson et al. [17] by using nonlinear Michaelis-Menten-Monod functions to describe the proliferation rate of the immune response and the removal rates of both the infected cells and free merozoites by immune effectors. They investigated the existence and local stability of the equilibrium points. They also proved that the model undergoes Hopf bifurcation at the positive equilibrium and shows periodic oscillations, where the proliferation rate of immune cells is taken as a bifurcation parameter. Tewa et al. [28] analyzed a general within-host model of malaria infection with immune response. The production rate of healthy erythrocytes and the stimulation rate of immune response are given by general functions with suitable hypotheses. Tchinda et al. [29] considered a particular form of the model studied in [28]. The dynamics of the immune response is modeled by "Allee effect", where the immune response cannot persist under certain thresholds. They explored the existence and local stability of model's equilibria. On the other hand, some works have differentiated between the immune response that targets infected red blood cells, and the immune response that targets free merozoites. For instance, Chen et al. [2] studied a blood-stage malaria infection model with two competitive strains of P. falciparum and two types of immune responses, namely cellular immunity and antibody immunity. The infection rates and the stimulation rates of the immune responses are provided by saturated rates. They identified the global stability conditions of the equilibrium point at which two strains of malaria parasite coexist. Cai et al. [30] proposed within-host models for symptomatically and asymptomatically infected hosts. The infected cells are divided into immature trophozoites and mature trophozoites. The model contains an immune response directed against mature trophozoites. The dynamics of this model was linked to a between-host model. Orwa et al. [31] formulated and analyzed a within-host model for malaria infection with malaria vaccines and CTL immune response. The infected red blood cells are classified into blood trophozoites and schizonts. They showed that the parameter which represents the average number of merozoites released per schizont infected red blood cell is the most sensitive parameter in their model. Song et al. [6] developed and analyzed a model to study the evolution of drug resistance parasites in the presence of a competition with drug sensitive parasites, immune response, and antimalarial drug. The effects of CTL and antibody immune responses on malaria infection are examined in this paper. The previous works did not take into account the spatial effects, and they assumed that the spatial distribution of cells and merozoites are homogeneous. The homogeneity assumption is unrealistic for many reasons: (i) many parasites during malaria infection accumulate in organs such as the brain and lungs [3]; (ii) there might be local restriction on the availability of uninfected red blood cells, which can reduce their contact with free merozoites [27]; (iii) the interaction between pathogens and the immune response tends to be local within the body of infected hosts [32]; (iv) red blood cells and merozoites are in motion [21]. Hence, more realistic models should take into consideration the non-uniform distribution of cells and merozoites across the body of the host in addition to their ability to move within the body [3,27]. In [33], Takoutsing et al. investigated a reaction-diffusion model for the within-host dynamics of malaria with a periodic antimalarial drug and immune factors. They examined the effect of treatment duration and efficacy on the distribution of parasites in space and time. Another factor that was not considered in the previous models is the perturbation of the blood-stage asexual life-cycle of malaria due to nutrient deprivation [3,34]. The degradation of host erythrocyte hemoglobin releases amino acids which represent a nutrient source for the malaria parasite P. falciparum. Nevertheless, human hemoglobin does not contain the single amino acid isoleucine which exists in more than 99% of the proteins encoded by the parasite [34]. Thus, the parasite must get isoleucine from an extracellular source for its survival. It has been shown that when the concentration of isoleucine is low, the parasite slows its metabolism and progresses from the ring stage to the trophozoite stage at a reduced rate. This helps the parasite to save energy and resources during isoleucine starvation. However, this starvation can protect the parasite for certain periods and the parasite may not be able to recover after long periods of starvation [34]. The spatial effects and the impact of isoleucine starvation on the growth of malaria parasites are considered in this paper.
The aim of this paper is to study a reaction-diffusion model for the within-host dynamics of malaria infection with both CTL and antibody immune responses. The model includes the effects of blood-stage antimalarial drugs and isoleucine starvation on parasite dynamics. Furthermore, the model captures the progression of infected red blood cells through the three developmental stages: ring stage, trophozoite stage, and schizont stage. Therefore, the model is considered to be an extension of the models studied in [23,31], where the effect of isoleucine starvation, spatial factors, or the ring stage are not involved.

A Reaction-Diffusion within-Host Malaria Model with both Cytotoxic T Lymphocyte and Antibody Immune Responses
In this section, we provide a complete description of the model under consideration. The proposed model captures the interaction of uninfected red blood cells, infected red blood cells, free blood merozoites, CTLs, and antibodies. In the formulation of the model, we make the following assumptions: (i) The dynamics of uninfected red blood cells is similar to the one given in model (1); (ii) The infected red blood cell population is divided into ring, trophozoite, and schizont infected red blood cells; (iii) The ring infected red blood cells increase due to the invasion of healthy red blood cells, and they decrease due to either the development into trophozoite infected red blood cells or natural death; (iv) The trophozoite infected red blood cells decrease due to either the development into schizont infected red blood cells or natural death; (v) The schizont infected red blood cells are depleted due to either the CTLs or natural death; (vi) Free merozoites in the blood increase as a result of schizonts bursting, and they decrease due to either the removal by antibodies or natural death; (vii) CTLs are stimulated by schizont infected red blood cells, while the production of antibodies is stimulated by free merozoites; (viii) CTLs and antibodies decrease duo to natural decay rates.
Hence, motivated by [23,31,33], we study the following reaction-diffusion malaria model for x ∈ Ω and t > 0, where Ω is a connected and bounded spatial domain with smooth boundary ∂Ω.
The variables U(x, t), N(x, t), T(x, t), S(x, t), M(x, t), W(x, t), and B(x, t) denote the concentrations of uninfected red blood cells, ring infected red blood cells, trophozoite infected red blood cells, schizont infected red blood cells, free merozoites, CTLs, and antibodies at position x and time t, respectively. We assume that red blood cells, free merozoites, and immune cells diffuse freely in the domain with different abilities. In the diffusion terms, D Λ is the diffusion coefficient of the component Λ, and ∆ = ∂ 2 ∂x 2 is the Laplacian operator. Ring infected red blood cells develop into trophozoite infected cells at rate γ 1 N, while trophozoites are developed into schizont infected cells at rate γ 2 T. CTLs attack and kill schizont infected red blood cells at rate α 1 SW, while they are stimulated at rate ρ 1 SW. Free merozoites are neutralized by antibodies at rate α 2 MB, and antibodies grow at rate ρ 2 MB. The death rates of rings, trophozoites, schizonts, free merozoites, CTLs, and antibodies are given by d 1 N, d 2 T, d 3 S, d 4 M, d 5 W, and d 6 B, respectively. The symbols ε 1 and ε 2 measure the efficacy of antimalarial drugs that affect the infection rate of healthy cells and the production rate of free merozoites, where 0 ≤ ε 1 < 1 and 0 ≤ ε 2 < 1. The symbol b measures the efficacy of isoleucine starvation, where 0 ≤ b < 1. The other parameters have the same biological meanings as those in model (1). The full description of all parameters is given in Table 1. We consider model (3) with the following initial conditions where i (x), for i=1,...,7, are non-negative and continuous functions. The boundary conditions of model (3) are reflected by the homogeneous Neumann boundary conditions where ∂ ∂ ς is the outward normal derivative on ∂Ω. This type of boundary conditions means that cells and free merozoites do not migrate from the isolated boundary. Efficacy of blood-stage treatment Varied b Efficacy of isoleucine starvation Varied ε 2 Efficacy of blood-stage treatment Varied - Diffusion coefficient of antibodies 0.2 Assumed (3) In the current section, we prove that all solutions of model (3) are non-negative and bounded. Furthermore, we address all equilibrium points and derive the threshold conditions required for their existence.

Basic Properties of Model
Let Y = BUC Ω , R 7 be the set of bounded continuous functions fromΩ to R 7 . Let Y + = BUC Ω , R 7 + be defined such that Y + ⊂ Y. The positive cone Y + induces a partial order on Y. Suppose that the norm is defined by z Y = sup x∈Ω |z(x)|, where | · | is the Euclidean norm on R 7 . Accordingly, the space (Y, · Y ) is a Banach lattice [35,36].

Proof. For any
We can see that Q is locally Lipschitz on Y + . We rewrite system (3)-(5) as the abstract differential equation It follows from [35][36][37] that system (3)-(5) has a unique non-negative solution on [0, σ max ) for any ∈ Y + , where [0, σ max ) is the maximal existence time interval on which the solution exists. To show boundedness, we define a function where Let Π 1 (t) be a solution to the ordinary differential equation system As a result, U(x, t) and N(x, t) are bounded. Next, we define a function .
where ω 2 = min {d 2 , d 3 , d 5 }. By using the comparison principle [38], we get This indicates that T(x, t), S(x, t), and W(x, t) are bounded. The final step is to prove the boundedness of M(x, t) and B(x, t). We define the function As where ω 3 = min {d 4 , d 6 }. We deduce from the comparison principle [38] that We conclude from the above analysis that U( . Therefore, the solutions are bounded onΩ × [0, +∞) depending on the standard theory for semi-linear parabolic systems [39].
The antibody-activated equilibrium The CTL/antibody coexistence equilibrium From Equation (11) we get W = 0 or S = d 5 ρ 1 . In addition, from Equation (12) From Equations (13)- (17) we get two equilibrium points. The first one is the disease-free Hence, the equilibrium E 1 exists if R i > 1. Here, R i is the within-host reproductive number and it represents the average number of new infected cells generated by one infected cell in a healthy cell population [18]. It follows that the condition R i > 1 is needed for successful establishment of malaria infection in the absence of CTL and antibody immune responses.
(ii) If W = 0 and M = M 2 = d 6 ρ 2 , then we get the antibody-activated equilibrium E 2 = (U 2 , N 2 , T 2 , S 2 , M 2 , 0, B 2 ). The other components of E 2 are given by Consequently, the equilibrium E 2 exists if R j > 1. The threshold condition R j > 1 is required for activating the antibody immune response during malaria infection.
. The other components of E 3 are given by , Thus, the equilibrium E 3 exists if R k > 1. The threshold condition R k > 1 is needed for activating the CTL immune response during malaria infection.
ρ 2 , then we get the CTL/antibody coexistence equilibrium . The other components of E 4 are given by , We note that These conditions are required for activating both the CTL and antibody immune responses during malaria infection.
In the next sections, we will use the following abbreviations (3) In this section, we show the global stability of the five equilibria of model (3) and the corresponding local instability conditions. For this purpose, we use the method developed in [40,41] for constructing appropriate Lyapunov functionals. It is worth noting that the function g(ϑ) = ϑ − 1 − ln ϑ ≥ 0 for ϑ > 0, and g(ϑ) = 0 ⇔ ϑ = 1.

Theorem 3.
The disease-free equilibrium E 0 is globally asymptotically stable if R i ≤ 1. It loses its stability when R i > 1.

Proof. Consider the Lyapunov functional candidate
Then, we get Accordingly, the time derivative of Θ 0 (t) is given by By using the Divergence theorem and the boundary conditions stated in (5), we get By using relation (19), the time derivative in (18) is reduced to It is easy to conclude from model (3) that the singleton {E 0 } is the largest invariant subset of According to LaSalle's invariance principle [42], we see that E 0 is globally asymptotically stable if R i ≤ 1. To show that E 0 loses its stability when R i > 1, we compute the characteristic equation. Let 0 = ν 1 < ν 2 < ... < ν n < .. be the complete set of eigenvalues of the Laplace operator −∆ with the homogeneous Neumann boundary conditions. Assume that E (ν i ) is the eigenfunction space corresponding to ν i (i = 1, 2, ...). Let P = dim E (ν i ) be the dimension of E (ν i ) and {ζ ij : j = 1, 2, ..., P} be an orthonormal basis of E (ν i ). Define H = {(U, N, T, S, M, W, B) ∈ C 1 (Ω) 7 : Thus, we get Define LF = D∆F + J (E l )F. Then, H i is invariant under the operator L for i ≥ 1. Moreover, λ is an eigenvalue of L if and only if it is an eigenvalue of −Dν i + J (E l ) for some i ≥ 1, for which there is an eigenvector in H i . Specifically, λ is a root of the characteristic equation where I denotes the identity matrix. To show the instability of any equilibrium point, it is enough to find i such that Equation (20) has a positive eigenvalue.

Theorem 4.
Suppose that R i > 1. Then, the immune-free equilibrium E 1 is globally asymptotically stable if R j ≤ 1 and R k ≤ 1. It is unstable if R j > 1 or R k > 1.

Proof. Take the following Lyapunov functional
Thus, we obtain The equilibrium conditions at E 1 are given by After using (23), Equation (22) is transformed into From the values of S 1 and S 3 , we find that From the values of M 1 and M 2 , we have (25) By using (19), (24) and (25), the time derivative of Θ 1 (t) is computed as The arithmetical and geometrical means relationship implies that 5 − where Two roots of Equation (26) are given by We see that Thus, we have at least one positive root if R k > 1 or R j > 1. As a consequence, E 1 is unstable if R k > 1 or R j > 1.

Theorem 5.
Suppose that R j > 1. Then, the antibody-activated equilibrium E 2 is globally asymptotically stable if R l ≤ 1. It loses its stability when R l > 1.

Proof. Define the following Lyapunov functional
Thus, we get At the equilibrium state, E 2 fulfills the following conditions After using (28), Equation (27) is transformed into From the equilibrium values, we get After using (19) and (29), the time derivative of Θ 2 (t) is given by We see that The global asymptotic stability of E 2 follows from LaSalle's invariance principle [42] when R j > 1 and R l ≤ 1. The characteristic equation at E 2 is provided as follows where One root of Equation (30) is given by This implies that E 2 becomes unstable when R l > 1.
Theorem 6. Suppose that R k > 1. Then, the CTL-activated equilibrium E 3 is globally asymptotically stable if Proof. Take the following Lyapunov functional Thus, we get At the equilibrium state, E 3 fulfills the following conditions After using (32), Equation (31) is simplified to By applying (19), the time derivative of Θ 3 (t) is given by We see that The characteristic equation at E 3 is computed as follows where

One root of Equation (33) is given by
Consider the case when i = 1, we get This implies that Equation (33) has a positive root and E 3 becomes unstable when R j R l > 1.
The CTL/antibody coexistence equilibrium E 4 is globally asymptotically stable if Proof. Define the following Lyapunov functional By taking the time partial derivative of Θ 4 (x, t), we get At the equilibrium state, E 4 satisfies the following conditions By using (19) and (34), the time derivative of Θ 4 (t) is given by It is easy to see that principle [42] that E 4 is globally asymptotically stable if R l > 1 and R j R l > 1.

Stability of Equilibria
In this subsection, we fix the values ε 1 = b = ε 2 = 0.1. We divide the simulations into five cases corresponding to the global stability of each one of the equilibrium points as follows: Case 1. We consider the values β = 2 × 10 −9 , ρ 1 = 2 × 10 −8 , and ρ 2 = 3 × 10 −7 . This gives R i = 0.1373 < 1. The equilibrium E 0 = (10 10 , 0, 0, 0, 0, 0, 0) is globally asymptotically stable in this situation (see Figure 1), which agrees with Theorem 3. This represents the case when the malaria infection is eliminated from the host's body. Thus, this point is the optimal objective for antimalarial drugs.
It is worth mentioning that the global stability of the equilibrium points assures that the long time behavior of solutions is not affected by initial condition [43], i.e., for any initial conditions satisfying (4), the long time behavior of solutions will be similar to the long time behavior of solutions obtained with the initial conditions (35) and shown in Figures 1-5.

Effect of Isoleucine Starvation and Drugs on the Malaria Dynamics
To see the effect of isoleucine starvation on the stability of equilibrium points, we consider the values of the parameters β = 2 × 10 −8 , ρ 1 = 2 × 10 −9 , ρ 2 = 3 × 10 −9 , ε 1 = 0.1, b = 0.8, and ε 2 = 0.1. In this situation, we get R i = 0.3462 < 1 which implies that the disease-free equilibrium E 0 is globally asymptotically stable. Mathematically, this means that increasing the value of b can decrease the value of R i to less than one, which destabilizes E 1 and stabilizes E 0 . The relation between b and R i is depicted in Figure 6a, where R i is drawn as a function of b. The dotted line denotes the values of b for which E 0 is unstable, while the solid line denotes the values of b for which E 0 is stable. Biologically, increasing isoleucine starvation efficacy forces the system to switch from the infection state to the malaria-free state. Accordingly, the impact of isoleucine starvation on system's behavior can help design more effective treatments.
We can see the effect of antimalarial drugs against infection by considering the values β = 2 × 10 −8 , ρ 1 = 2 × 10 −9 , ρ 2 = 3 × 10 −9 , ε 1 = 0.8, b = 0.1, and ε 2 = 0.1. This gives R i = 0.3051 < 1, where the disease-free equilibrium E 0 becomes globally asymptotically stable. Therefore, increasing the efficacy of treatment to large values decreases the infection rate and, as a result, clears the disease (see Figure 6b). A similar result can be obtained by increasing the value of ε 2 . Therefore, the infection can be cleared with highly effective treatments.

Discussion
In this paper, we investigated a reaction-diffusion model for the blood-stage dynamics of malaria infection with CTL and antibody immune responses. The model studies the interactions between uninfected erythrocytes, three types of infected red blood cells, free merozoites, CTL immune response and antibody immune response. The model contains parameters to measure the efficacy of antimalarial treatments and amino acid starvation. It has five equilibrium points: (a) The disease-free equilibrium E 0 is usually defined and globally asymptotically stable if R i ≤ 1.
This point corresponds to the elimination of malaria infection. (b) The immune-free equilibrium E 1 is defined if R i > 1, and it is globally asymptotically stable if R j ≤ 1 and R k ≤ 1. At this point, the malaria parasite P. falciparum succeeds in establishing the infection when the immune responses are not active. (c) The antibody-activated equilibrium E 2 is defined if R j > 1, and it is globally asymptotically stable if R l ≤ 1. At this point, the antibody immune response is activated to attack the free blood merozoites. (d) The CTL-activated equilibrium E 3 is defined if R k > 1, and it is globally asymptotically stable if R j R l ≤ 1. At this point, the CTL immune response is activated to kill the schizont infected red blood cells. (e) The CTL/antibody coexistence equilibrium E 4 is defined and globally asymptotically stable if R l > 1 and R j R l > 1. Here, the malaria infection stimulates the CTL and antibody immune responses to fight the infection.
We found that increasing isoleucine starvation efficacy b in model (3) drives the system towards the malaria-free equilibrium and, in general, affects the stability of equilibria. Thus, the perturbation in the asexual malaria life-cycle duo to isoleucine starvation can be an important target for antimalarial drugs [34]. The effectiveness of these drugs depends on the value the isoleucine starvation efficacy b could reach when these drugs are used. Isoleucine starvation was produced in vitro by exposing the parasites to isoleucine-free medium [34]. Also, protein malnutrition was produced in animal studies by feeding protein-restricted diets to malaria-infected rats [44,45]. However, the induction of isoleucine starvation in humans and its complete effect on malaria patients is still under investigation. In addition to isoleucine starvation, we found that blood-stage treatments with high efficacy (ε 1 and ε 2 ) are needed to clear the infection. Hence, the values of these parameters can have a crucial impact on the disease dynamics. Furthermore, the immune responses at the CTL/antibody coexistence equilibrium E 4 do not clear the infection, but they reduce the concentrations of free merozoites and infected cells at the three developmental stages (rings, trophozoites, and schizonts). The presence of diffusion in model (3) does not affect the global stability of equilibria. However, it affects the behavior of solutions at the beginning of infection before reaching the equilibrium points (see Figure 7). This may have a substantial impact on the estimates of parameters describing the early infection [32]. Thus, considering spatial effects is important to get a good approximation for infection dynamics. The spatial distribution of solutions at different progressive time points can be seen more clearly by considering two-dimensional square domain [46]. Compared to the existing models of blood-stage malaria infection, our model is the first model that addresses with a detailed mathematical analysis the development of the infected red blood cells through the three main sequential stages using a reaction-diffusion model with adaptive immune response. Our model also includes the effects of antimalarial drugs and isoleucine starvation on malaria dynamics. Previous models have ignored some of these factors. The results can help design more effective treatments to fight malaria infection by taking into account some important factors that can affect the disease progression. Model (3) uses bilinear rates and ignores the effects of delays or immune responses against the ring and schizont stages. Thus, model (3) can be improved in many ways: (i) by considering more general infection and stimulation rates (see, e.g. [47,48]); (ii) by assuming CTL immune response against the three types of infected cells; (iii) by taking into account the delays that may occur in some processes during malaria infection (see, e.g. [49]); (iv) by performing a bifurcation analysis of the model (see, e.g. [50]) . We leave these improvements as possible future works.

Conflicts of Interest:
The authors declare no conflict of interest.