Next Article in Journal
Sensitivity Analysis of Mathematical Model to Study the Effect of T Cells Infusion in Treatment of CLL
Previous Article in Journal
Sub-Mesoscale Frontal Instabilities in the Omani Coastal Current
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Department of Mathematics, Faculty of Science, King Abdulaziz University, P.O. Box 80203, Jeddah 21589, Saudi Arabia
*
Author to whom correspondence should be addressed.
Mathematics 2020, 8(4), 563; https://doi.org/10.3390/math8040563
Submission received: 9 February 2020 / Revised: 27 March 2020 / Accepted: 3 April 2020 / Published: 11 April 2020

Abstract

:
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.

1. 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
d U ( t ) d t = μ d ˜ U ( t ) β U ( t ) M ( t ) , d Y ( t ) d t = β U ( t ) M ( t ) e Y ( t ) , d M ( t ) d t = η e Y ( t ) f M ( t ) g β U ( t ) M ( t ) ,
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 rate d ˜ U . Infected red blood cells die at rate e Y and rupture to produce η merozoites per infected cell. The total production rate of merozoites is given by η e Y . 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
d U ( t ) d t = μ d ˜ U ( t ) β U ( t ) M ( t ) , d Y ( t ) d t = β U ( t ) M ( t ) e Y ( t ) α ˜ 1 Y ( t ) Z ( t ) , d M ( t ) d t = η e Y ( t ) f M ( t ) g β U ( t ) M ( t ) α ˜ 2 M ( t ) Z ( t ) , d Z ( t ) d t = ρ ˜ 1 Y ( t ) Z ( t ) + ρ ˜ 2 M ( t ) Z ( t ) h Z ( t ) ,
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 Y Z and α ˜ 2 M Z , respectively. The infected red blood cells and free merozoites stimulate the production of immune cells at rates ρ ˜ 1 Y Z and ρ ˜ 2 M Z , respectively. The death rate of immune cells is given by h Z . 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.

2. 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
U ( x , t ) t = D U Δ U ( x , t ) + μ d ˜ U ( x , t ) β ( 1 ε 1 ) U ( x , t ) M ( x , t ) , N ( x , t ) t = D N Δ N ( x , t ) + β ( 1 ε 1 ) U ( x , t ) M ( x , t ) γ 1 ( 1 b ) N ( x , t ) d 1 N ( x , t ) , T ( x , t ) t = D T Δ T ( x , t ) + γ 1 ( 1 b ) N ( x , t ) γ 2 T ( x , t ) d 2 T ( x , t ) , S ( x , t ) t = D S Δ S ( x , t ) + γ 2 T ( x , t ) α 1 S ( x , t ) W ( x , t ) d 3 S ( x , t ) , M ( x , t ) t = D M Δ M ( x , t ) + η d 3 ( 1 ε 2 ) S ( x , t ) α 2 M ( x , t ) B ( x , t ) d 4 M ( x , t ) , W ( x , t ) t = D W Δ W ( x , t ) + ρ 1 S ( x , t ) W ( x , t ) d 5 W ( x , t ) , B ( x , t ) t = D B Δ B ( x , t ) + ρ 2 M ( x , t ) B ( x , t ) d 6 B ( x , t ) ,
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 S W , while they are stimulated at rate ρ 1 S W . Free merozoites are neutralized by antibodies at rate α 2 M B , and antibodies grow at rate ρ 2 M B . 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
U ( x , 0 ) = ϖ 1 ( x ) , N ( x , 0 ) = ϖ 2 ( x ) , T ( x , 0 ) = ϖ 3 ( x ) , S ( x , 0 ) = ϖ 4 ( x ) , M ( x , 0 ) = ϖ 5 ( x ) , W ( x , 0 ) = ϖ 6 ( x ) , B ( x , 0 ) = ϖ 7 ( x ) , x Ω ¯ ,
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
U ς = N ς = T ς = S ς = M ς = W ς = B ς = 0 , t > 0 , x Ω ,
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.

3. Basic Properties of Model (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.
Let Y = B U C Ω ¯ , R 7 be the set of bounded continuous functions from Ω ¯ to R 7 . Let Y + = B U C Ω ¯ , 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].
Theorem 1.
Suppose that D U = D N , D T = D S = D W , and D M = D B . For any initial conditions satisfying (4), there exists a unique solution of system (3)–(5) which is non-negative and bounded on Ω ¯ × [ 0 , + ) .
Proof. 
For any ϖ = ( ϖ 1 , ϖ 2 , ϖ 3 , ϖ 4 , ϖ 5 , ϖ 6 , ϖ 7 ) T Y + , we define Q = ( Q 1 , Q 2 , Q 3 , Q 4 , Q 5 , Q 6 , Q 7 ) T : Y + Y by
Q 1 ( ϖ ) ( x ) = μ d ˜ ϖ 1 ( x ) β ( 1 ε 1 ) ϖ 1 ( x ) ϖ 5 ( x ) , Q 2 ( ϖ ) ( x ) = β ( 1 ε 1 ) ϖ 1 ( x ) ϖ 5 ( x ) γ 1 ( 1 b ) ϖ 2 ( x ) d 1 ϖ 2 ( x ) , Q 3 ( ϖ ) ( x ) = γ 1 ( 1 b ) ϖ 2 ( x ) γ 2 ϖ 3 ( x ) d 2 ϖ 3 ( x ) , Q 4 ( ϖ ) ( x ) = γ 2 ϖ 3 ( x ) α 1 ϖ 4 ( x ) ϖ 6 ( x ) d 3 ϖ 4 ( x ) , Q 5 ( ϖ ) ( x ) = η d 3 ( 1 ε 2 ) ϖ 4 ( x ) α 2 ϖ 5 ( x ) ϖ 7 ( x ) d 4 ϖ 5 ( x ) , Q 6 ( ϖ ) ( x ) = ρ 1 ϖ 4 ( x ) ϖ 6 ( x ) d 5 ϖ 6 ( x ) , Q 7 ( ϖ ) ( x ) = ρ 2 ϖ 5 ( x ) ϖ 7 ( x ) d 6 ϖ 7 ( x ) .
We can see that Q is locally Lipschitz on Y + . We rewrite system (3)–(5) as the abstract differential equation
d G d t = D G + Q ( G ) , t > 0 , G ( 0 ) = ϖ Y + ,
where G = ( U , N , T , S , M , W , B ) T and D G = ( D U Δ U , D N Δ N , D T Δ T , D S Δ S , D M Δ M , D W Δ W , D B Δ B ) T . 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
Π 1 ( x , t ) = U ( x , t ) + N ( x , t ) .
As D U = D N , we get
Π 1 ( x , t ) t D U Δ Π 1 ( x , t ) = μ d ˜ U ( x , t ) γ 1 ( 1 b ) + d 1 N ( x , t ) μ ω 1 U ( x , t ) + N ( x , t ) = μ ω 1 Π 1 ( x , t ) ,
where ω 1 = min d ˜ , γ 1 ( 1 b ) + d 1 . Thus, Π 1 ( x , t ) fulfills the following system
Π 1 ( x , t ) t D U Δ Π 1 ( x , t ) μ ω 1 Π 1 ( x , t ) , Π 1 ς = 0 , Π 1 ( x , 0 ) = ϖ 1 ( x ) + ϖ 2 ( x ) 0 .
Let Π ˜ 1 ( t ) be a solution to the ordinary differential equation system
d Π ˜ 1 ( t ) d t = μ ω 1 Π ˜ 1 ( t ) , Π ˜ 1 ( 0 ) = max x Ω ¯ { ϖ 1 ( x ) + ϖ 2 ( x ) } .
This implies that Π ˜ 1 ( t ) max μ ω 1 , max x Ω ¯ { ϖ 1 ( x ) + ϖ 2 ( x ) } . According to the comparison principle [38], we have Π 1 ( x , t ) Π ˜ 1 ( t ) . Hence, we get
Π 1 ( x , t ) max μ ω 1 , max x Ω ¯ { ϖ 1 ( x ) + ϖ 2 ( x ) } : = K 1 .
As a result, U ( x , t ) and N ( x , t ) are bounded. Next, we define a function
Π 2 ( x , t ) = T ( x , t ) + S ( x , t ) + α 1 ρ 1 W ( x , t ) .
As D T = D S = D W and N ( x , t ) K 1 , we obtain
Π 2 ( x , t ) t D T Δ Π 2 ( x , t ) = γ 1 ( 1 b ) N ( x , t ) d 2 T ( x , t ) d 3 S ( x , t ) α 1 d 5 ρ 1 W ( x , t ) γ 1 ( 1 b ) K 1 d 2 T ( x , t ) d 3 S ( x , t ) α 1 d 5 ρ 1 W ( x , t ) γ 1 ( 1 b ) K 1 ω 2 Π 2 ( x , t ) ,
where ω 2 = min d 2 , d 3 , d 5 . By using the comparison principle [38], we get
Π 2 ( x , t ) max γ 1 ( 1 b ) K 1 ω 2 , max x Ω ¯ { ϖ 3 ( x ) + ϖ 4 ( x ) + α 1 ρ 1 ϖ 6 ( x ) } : = K 2 .
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
Π 3 ( x , t ) = M ( x , t ) + α 2 ρ 2 B ( x , t ) .
As D M = D B and S ( x , t ) K 2 , we have
Π 3 ( x , t ) t D M Δ Π 3 ( x , t ) = η d 3 ( 1 ε 2 ) S ( x , t ) d 4 M ( x , t ) α 2 d 6 ρ 2 B ( x , t ) η d 3 ( 1 ε 2 ) K 2 d 4 M ( x , t ) α 2 d 6 ρ 2 B ( x , t ) η d 3 ( 1 ε 2 ) K 2 ω 3 Π 3 ( x , t ) ,
where ω 3 = min d 4 , d 6 . We deduce from the comparison principle [38] that
Π 3 ( x , t ) max η d 3 ( 1 ε 2 ) K 2 ω 3 , max x Ω ¯ { ϖ 5 ( x ) + α 2 ρ 2 ϖ 7 ( x ) } .
We conclude from the above analysis that U ( x , t ) , N ( x , t ) , T ( x , t ) , S ( x , t ) , M ( x , t ) , W ( x , t ) , and B ( x , t ) are bounded on [ 0 , σ max ) . Therefore, the solutions are bounded on Ω ¯ × [ 0 , + ) depending on the standard theory for semi-linear parabolic systems [39]. □
Theorem 2.
There exist threshold parameters R i , R j , R k , and R l such that model (3) has five equilibrium points:
(a)
The disease-free equilibrium E 0 = ( U 0 , 0 , 0 , 0 , 0 , 0 , 0 ) always exists.
(b)
The immune-free equilibrium E 1 = ( U 1 , N 1 , T 1 , S 1 , M 1 , 0 , 0 ) exists if R i > 1 .
(c)
The antibody-activated equilibrium E 2 = ( U 2 , N 2 , T 2 , S 2 , M 2 , 0 , B 2 ) exists if R j > 1 .
(d)
The CTL-activated equilibrium E 3 = ( U 3 , N 3 , T 3 , S 3 , M 3 , W 3 , 0 ) exists if R k > 1 .
(e)
The CTL/antibody coexistence equilibrium E 4 = ( U 4 , N 4 , T 4 , S 4 , M 4 , W 4 , B 4 ) exists if R l > 1 and R j R l > 1 .
Proof. 
Any equilibrium point E = ( U , N , T , S , M , W , B ) of model (3) satisfies
μ d ˜ U β ( 1 ε 1 ) U M = 0 ,
β ( 1 ε 1 ) U M γ 1 ( 1 b ) N d 1 N = 0 ,
γ 1 ( 1 b ) N γ 2 T d 2 T = 0 ,
γ 2 T α 1 S W d 3 S = 0 ,
η d 3 ( 1 ε 2 ) S α 2 M B d 4 M = 0 ,
ρ 1 S W d 5 W = 0 ,
ρ 2 M B d 6 B = 0 .
From Equation (11) we get W = 0 or S = d 5 ρ 1 . In addition, from Equation (12) we obtain B = 0 or M = d 6 ρ 2 . Thus, we have four cases:
(i)
If W = 0 and B = 0 , then the equilibrium conditions are shortened to
μ d ˜ U β ( 1 ε 1 ) U M = 0 ,
β ( 1 ε 1 ) U M γ 1 ( 1 b ) N d 1 N = 0 ,
γ 1 ( 1 b ) N γ 2 T d 2 T = 0 ,
γ 2 T d 3 S = 0 ,
η d 3 ( 1 ε 2 ) S d 4 M = 0 .
From Equations (13)–(17) we get two equilibrium points. The first one is the disease-free equilibrium E 0 = ( U 0 , 0 , 0 , 0 , 0 , 0 , 0 ) , where U 0 = μ d ˜ . As U 0 > 0 , the equilibrium E 0 always exists. The second one is the immune-free equilibrium E 1 = ( U 1 , N 1 , T 1 , S 1 , M 1 , 0 , 0 ) , where
U 1 = d 4 γ 2 + d 2 γ 1 ( 1 b ) + d 1 β η γ 1 γ 2 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) , N 1 = d ˜ d 4 γ 2 + d 2 R i 1 β η γ 1 γ 2 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) , T 1 = d ˜ d 4 R i 1 β η γ 2 ( 1 ε 1 ) ( 1 ε 2 ) , S 1 = d ˜ d 4 R i 1 β η d 3 ( 1 ε 1 ) ( 1 ε 2 ) , M 1 = d ˜ R i 1 β ( 1 ε 1 ) ,
and
R i = β η μ γ 1 γ 2 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) d ˜ d 4 γ 2 + d 2 γ 1 ( 1 b ) + d 1 .
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
U 2 = μ ρ 2 d ˜ ρ 2 + β d 6 ( 1 ε 1 ) , N 2 = β μ d 6 ( 1 ε 1 ) γ 1 ( 1 b ) + d 1 d ˜ ρ 2 + β d 6 ( 1 ε 1 ) , T 2 = β μ γ 1 d 6 ( 1 b ) ( 1 ε 1 ) γ 2 + d 2 γ 1 ( 1 b ) + d 1 d ˜ ρ 2 + β d 6 ( 1 ε 1 ) , S 2 = β μ γ 1 γ 2 d 6 ( 1 b ) ( 1 ε 1 ) d 3 γ 2 + d 2 γ 1 ( 1 b ) + d 1 d ˜ ρ 2 + β d 6 ( 1 ε 1 ) , B 2 = d 4 R j 1 α 2 ,
where
R j = β η μ γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) d 4 γ 2 + d 2 γ 1 ( 1 b ) + d 1 d ˜ ρ 2 + β d 6 ( 1 ε 1 ) .
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.
(iii)
If B = 0 and S = S 3 = d 5 ρ 1 , then we obtain the CTL-activated equilibrium E 3 = ( U 3 , N 3 , T 3 , S 3 , M 3 , W 3 , 0 ) . The other components of E 3 are given by
U 3 = μ ρ 1 d 4 ρ 1 d ˜ d 4 + β η d 3 d 5 ( 1 ε 1 ) ( 1 ε 2 ) , N 3 = β η μ d 3 d 5 ( 1 ε 1 ) ( 1 ε 2 ) γ 1 ( 1 b ) + d 1 ρ 1 d ˜ d 4 + β η d 3 d 5 ( 1 ε 1 ) ( 1 ε 2 ) , T 3 = β η μ γ 1 d 3 d 5 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) γ 2 + d 2 γ 1 ( 1 b ) + d 1 ρ 1 d ˜ d 4 + β η d 3 d 5 ( 1 ε 1 ) ( 1 ε 2 ) , M 3 = η d 3 d 5 ( 1 ε 2 ) ρ 1 d 4 , W 3 = d 3 R k 1 α 1 ,
where
R k = β η μ ρ 1 γ 1 γ 2 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) γ 2 + d 2 γ 1 ( 1 b ) + d 1 ρ 1 d ˜ d 4 + β η d 3 d 5 ( 1 ε 1 ) ( 1 ε 2 ) .
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.
(iv)
If S = S 4 = d 5 ρ 1 and M = M 4 = d 6 ρ 2 , then we get the CTL/antibody coexistence equilibrium E 4 = ( U 4 , N 4 , T 4 , S 4 , M 4 , W 4 , B 4 ) . The other components of E 4 are given by
U 4 = μ ρ 2 d ˜ ρ 2 + β d 6 ( 1 ε 1 ) , N 4 = β μ d 6 ( 1 ε 1 ) γ 1 ( 1 b ) + d 1 d ˜ ρ 2 + β d 6 ( 1 ε 1 ) , T 4 = β μ γ 1 d 6 ( 1 b ) ( 1 ε 1 ) γ 2 + d 2 γ 1 ( 1 b ) + d 1 d ˜ ρ 2 + β d 6 ( 1 ε 1 ) , W 4 = d 3 R l 1 α 1 , B 4 = d 4 α 2 R j R l 1 ,
where
R l = β μ ρ 1 γ 1 γ 2 d 6 ( 1 b ) ( 1 ε 1 ) d 3 d 5 γ 2 + d 2 γ 1 ( 1 b ) + d 1 d ˜ ρ 2 + β d 6 ( 1 ε 1 ) ,
and
R j R l = η ρ 2 d 3 d 5 ( 1 ε 2 ) ρ 1 d 4 d 6 .
We note that W 4 > 0 if R l > 1 , and B 4 > 0 if R j R l > 1 . Thus, the equilibrium E 4 exists if R l > 1 and R j R l > 1 . 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
U ( x , t ) U , N ( x , t ) N , T ( x , t ) T , S ( x , t ) S , M ( x , t ) M , W ( x , t ) W , B ( x , t ) B .

4. Global Properties of Model (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
Θ 0 ( t ) = Ω Θ ˜ 0 ( x , t ) d x ,
where
Θ ˜ 0 ( x , t ) = U 0 U U 0 1 ln U U 0 + N + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) T + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) S + γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) M + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) B .
Then, we get
Θ ˜ 0 t = 1 U 0 U D U Δ U + μ d ˜ U β ( 1 ε 1 ) U M + D N Δ N + β ( 1 ε 1 ) U M γ 1 ( 1 b ) N d 1 N + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) D T Δ T + γ 1 ( 1 b ) N γ 2 T d 2 T + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) D S Δ S + γ 2 T α 1 S W d 3 S + γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) D M Δ M + η d 3 ( 1 ε 2 ) S α 2 M B d 4 M + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) D W Δ W + ρ 1 S W d 5 W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) D B Δ B + ρ 2 M B d 6 B = d ˜ U U 0 2 U + β ( 1 ε 1 ) U 0 d 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) M α 1 d 5 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) W α 2 d 6 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) B + 1 U 0 U D U Δ U + D N Δ N + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) D T Δ T + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) D S Δ S + γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) D M Δ M + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) D W Δ W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) D B Δ B .
Accordingly, the time derivative of Θ 0 ( t ) is given by
d Θ 0 d t = d ˜ Ω U U 0 2 U d x + d 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) R i 1 Ω M d x α 1 d 5 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) Ω W d x α 2 d 6 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) Ω B d x + D U Ω 1 U 0 U Δ U d x + D N Ω Δ N d x + D T γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) Ω Δ T d x + D S γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) Ω Δ S d x + D M γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) Ω Δ M d x + D W α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) Ω Δ W d x + D B α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) Ω Δ B d x .
By using the Divergence theorem and the boundary conditions stated in (5), we get
0 = Ω κ · ς d x = Ω div ( κ ) d x = Ω Δ κ d x , 0 = Ω 1 κ κ · ς d x = Ω div ( 1 κ κ ) d x = Ω Δ κ κ κ 2 κ 2 d x , for κ { U , N , T , S , M , W , B } .
This gives the following two relations
Ω Δ κ d x = 0 , Ω Δ κ κ d x = Ω κ 2 κ 2 d x , for κ { U , N , T , S , M , W , B } .
By using relation (19), the time derivative in (18) is reduced to
d Θ 0 d t = d ˜ Ω U U 0 2 U d x + d 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) R i 1 Ω M d x α 1 d 5 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) Ω W d x α 2 d 6 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) Ω B d x D U U 0 Ω U 2 U 2 d x .
We note that d Θ 0 d t 0 if R i 1 , and d Θ 0 d t = 0 if U = U 0 , M = 0 , W = 0 , and B = 0 . It is easy to conclude from model (3) that the singleton { E 0 } is the largest invariant subset of { ( U , N , T , S , M , W , B ) | d Θ 0 d t = 0 } . 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 ¯ = d i m E ( ν i ) be the dimension of E ( ν i ) and { ζ i j : j = 1 , 2 , , P ¯ } be an orthonormal basis of E ( ν i ) . Define
H = { ( U , N , T , S , M , W , B ) C 1 ( Ω ¯ ) 7 : U ς = N ς = T ς = S ς = M ς = W ς = B ς = 0 on Ω } , H i j = { a ζ i j : a R 7 } .
Thus, we get
H i = j = 1 P ¯ H i j and H = i = 1 H i .
Let E l = ( U l , N l , T l , S l , M l , W l , B l ) be any equilibrium point of model (3). The linearization of model (3) at E l is given by
F t = D Δ F + J ( E l ) F ,
where F = ( U , N , T , S , M , W , B ) T , D = d i a g ( D U , D N , D T , D S , D M , D W , D B ) and
J ( E l ) = d ˜ β ( 1 ε 1 ) M l 0 0 0 β ( 1 ε 1 ) U l 0 0 β ( 1 ε 1 ) M l γ 1 ( 1 b ) + d 1 0 0 β ( 1 ε 1 ) U l 0 0 0 γ 1 ( 1 b ) γ 2 + d 2 0 0 0 0 0 0 γ 2 α 1 W l d 3 0 α 1 S l 0 0 0 0 η d 3 ( 1 ε 2 ) α 2 B l d 4 0 α 2 M l 0 0 0 ρ 1 W l 0 ρ 1 S l d 5 0 0 0 0 0 ρ 2 B l 0 ρ 2 M l d 6 .
Define L F = 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
det ( λ I + D ν i J ( E l ) ) = 0 ,
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.
The characteristic equation at E 0 is computed as
λ + d ˜ + D U ν i λ + d 5 + D W ν i λ + d 6 + D B ν i f 0 ( λ ) = 0 ,
where
f 0 ( λ ) = λ 4 + λ 3 γ 1 ( 1 b ) + d 1 + D N ν i + γ 2 + d 2 + D T ν i + d 3 + D S ν i + d 4 + D M ν i + λ 2 ( γ 1 ( 1 b ) + d 1 + D N ν i + d 4 + D M ν i + γ 2 + d 2 + D T ν i d 3 + D S ν i + γ 1 ( 1 b ) + d 1 + D N ν i γ 2 + d 2 + D T ν i + d 3 + D S ν i d 4 + D M ν i ) + λ ( γ 1 ( 1 b ) + d 1 + D N ν i + d 4 + D M ν i γ 2 + d 2 + D T ν i d 3 + D S ν i + + γ 1 ( 1 b ) + d 1 + D N ν i γ 2 + d 2 + D T ν i + d 3 + D S ν i d 4 + D M ν i ) + γ 1 ( 1 b ) + d 1 + D N ν i γ 2 + d 2 + D T ν i d 3 + D S ν i d 4 + D M ν i β η U 0 γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) .
Equation (21) has the eigenvalues λ = d ˜ D U ν i < 0 , λ = d 5 D W ν i < 0 , and λ = d 6 D B ν i < 0 . The other roots are given by f 0 ( λ ) = 0 . Consider the situation when i = 1 corresponding to ν 1 = 0 , then we note
lim λ + f 0 ( λ ) = + , f 0 ( 0 ) | i = 1 = d 3 d 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 β η U 0 γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) = d 3 d 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 R i 1 .
Of note, f 0 ( 0 ) | i = 1 < 0 if R i > 1 . This means that Equation (21) has a positive root if R i > 1 . As a result, E 0 becomes unstable when R i > 1 . □
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
Θ 1 ( t ) = Ω Θ ˜ 1 ( x , t ) d x ,
where
Θ ˜ 1 ( x , t ) = U 1 U U 1 1 ln U U 1 + N 1 N N 1 1 ln N N 1 + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) T 1 T T 1 1 ln T T 1 + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) S 1 S S 1 1 ln S S 1 + γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) M 1 M M 1 1 ln M M 1 + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) B .
Thus, we obtain
Θ ˜ 1 t = 1 U 1 U D U Δ U + μ d ˜ U β ( 1 ε 1 ) U M + 1 N 1 N D N Δ N + β ( 1 ε 1 ) U M γ 1 ( 1 b ) N d 1 N + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) 1 T 1 T D T Δ T + γ 1 ( 1 b ) N γ 2 T d 2 T + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) 1 S 1 S D S Δ S + γ 2 T α 1 S W d 3 S + γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) 1 M 1 M D M Δ M + η d 3 ( 1 ε 2 ) S α 2 M B d 4 M + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) D W Δ W + ρ 1 S W d 5 W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) D B Δ B + ρ 2 M B d 6 B = 1 U 1 U μ d ˜ U + γ 1 ( 1 b ) + d 1 N 1 + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 ( 1 b ) T 1 + d 3 γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) S 1 + d 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) M 1 β ( 1 ε 1 ) U M N 1 N γ 1 ( 1 b ) + d 1 N T 1 T γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 ( 1 b ) T S 1 S d 3 γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) S M 1 M + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) S 1 W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) M 1 B α 1 d 5 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) W α 2 d 6 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) B + 1 U 1 U D U Δ U + 1 N 1 N D N Δ N + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) 1 T 1 T D T Δ T + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) 1 S 1 S D S Δ S + γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) 1 M 1 M D M Δ M + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) D W Δ W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) D B Δ B .
The equilibrium conditions at E 1 are given by
μ = d ˜ U 1 + β ( 1 ε 1 ) U 1 M 1 , β ( 1 ε 1 ) U 1 M 1 = γ 1 ( 1 b ) + d 1 N 1 = γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 ( 1 b ) T 1 = d 3 γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) S 1 , β ( 1 ε 1 ) U 1 M 1 = d 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) M 1 .
After using (23), Equation (22) is transformed into
Θ ˜ 1 t = 1 U 1 U d ˜ U 1 d ˜ U + β ( 1 ε 1 ) U 1 M 1 5 U 1 U N T 1 N 1 T T S 1 T 1 S S M 1 S 1 M U N 1 M U 1 N M 1 + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) S 1 S 3 W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) M 1 M 2 B + 1 U 1 U D U Δ U + 1 N 1 N D N Δ N + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) 1 T 1 T D T Δ T + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) 1 S 1 S D S Δ S + γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) 1 M 1 M D M Δ M + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) D W Δ W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) D B Δ B .
From the values of S 1 and S 3 , we find that
S 1 S 3 = β η μ ρ 1 γ 1 γ 2 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 d ˜ d 4 + β η d 3 d 5 ( 1 ε 1 ) ( 1 ε 2 ) β η ρ 1 d 3 ( 1 ε 1 ) ( 1 ε 2 ) γ 1 ( 1 b ) + d 1 γ 2 + d 2 = ρ 1 d ˜ d 4 + β η d 3 d 5 ( 1 ε 1 ) ( 1 ε 2 ) β η ρ 1 d 3 ( 1 ε 1 ) ( 1 ε 2 ) R k 1 .
From the values of M 1 and M 2 , we have
M 1 M 2 = β η μ γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) d 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d ˜ ρ 2 + β d 6 ( 1 ε 1 ) β ρ 2 d 4 ( 1 ε 1 ) γ 1 ( 1 b ) + d 1 γ 2 + d 2 = d ˜ ρ 2 + β d 6 ( 1 ε 1 ) β ρ 2 ( 1 ε 1 ) R j 1 .
By using (19), (24) and (25), the time derivative of Θ 1 ( t ) is computed as
d Θ 1 d t = d ˜ Ω U U 1 2 U d x + β ( 1 ε 1 ) U 1 M 1 Ω 5 U 1 U N T 1 N 1 T T S 1 T 1 S S M 1 S 1 M U N 1 M U 1 N M 1 d x + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 d ˜ d 4 + β η d 3 d 5 ( 1 ε 1 ) ( 1 ε 2 ) β η ρ 1 γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) R k 1 Ω W d x + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d ˜ ρ 2 + β d 6 ( 1 ε 1 ) β η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) R j 1 Ω B d x D U U 1 Ω U 2 U 2 d x D N N 1 Ω N 2 N 2 d x D T T 1 γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) Ω T 2 T 2 d x D S S 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) Ω S 2 S 2 d x D M M 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) Ω M 2 M 2 d x .
The arithmetical and geometrical means relationship implies that 5 U 1 U N T 1 N 1 T T S 1 T 1 S S M 1 S 1 M U N 1 M U 1 N M 1 0 . Thus, we see that d Θ 1 d t 0 if R j 1 and R k 1 . Furthermore, d Θ 1 d t = 0 if U = U 1 , N = N 1 , T = T 1 , S = S 1 , M = M 1 , W = 0 , and B = 0 . Hence, the singleton { E 1 } is the largest invariant subset of { ( U , N , T , S , M , W , B ) | d Θ 1 d t = 0 } . It follows from LaSalle’s invariance principle [42] that E 1 is globally asymptotically stable if R j 1 and R k 1 , where the point is defined for R i > 1 .
The characteristic equation at E 1 is computed as
λ ρ 1 S 1 + d 5 + D W ν i λ ρ 2 M 1 + d 6 + D B ν i f 1 ( λ ) = 0 ,
where
f 1 ( λ ) = λ + γ 1 ( 1 b ) + d 1 + D N ν i λ + γ 2 + d 2 + D T ν i λ + d 3 + D S ν i λ + d 4 + D M ν i × λ + β ( 1 ε 1 ) M 1 + d ˜ + D U ν i β η γ 1 γ 2 d 3 U 1 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) λ + d ˜ + D U ν i .
Two roots of Equation (26) are given by λ 1 = ρ 1 S 1 d 5 D W ν i and λ 2 = ρ 2 M 1 d 6 D B ν i . Consider i = 1 , then we get
λ 1 | i = 1 = ρ 1 S 1 S 3 = ρ 1 d ˜ d 4 + β η d 3 d 5 ( 1 ε 1 ) ( 1 ε 2 ) β η d 3 ( 1 ε 1 ) ( 1 ε 2 ) R k 1 . λ 2 | i = 1 = ρ 2 M 1 M 2 = d ˜ ρ 2 + β d 6 ( 1 ε 1 ) β ( 1 ε 1 ) R j 1 .
We see that λ 1 | i = 1 > 0 if R k > 1 , and λ 2 | i = 1 > 0 if R j > 1 . 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
Θ 2 ( t ) = Ω Θ ˜ 2 ( x , t ) d x ,
where
Θ ˜ 2 ( x , t ) = U 2 U U 2 1 ln U U 2 + N 2 N N 2 1 ln N N 2 + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) T 2 T T 2 1 ln T T 2 + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) S 2 S S 2 1 ln S S 2 + γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) M 2 M M 2 1 ln M M 2 + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) B 2 B B 2 1 ln B B 2 .
Thus, we get
Θ ˜ 2 t = 1 U 2 U D U Δ U + μ d ˜ U β ( 1 ε 1 ) U M + 1 N 2 N D N Δ N + β ( 1 ε 1 ) U M γ 1 ( 1 b ) N d 1 N + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) 1 T 2 T D T Δ T + γ 1 ( 1 b ) N γ 2 T d 2 T + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) 1 S 2 S D S Δ S + γ 2 T α 1 S W d 3 S + γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) 1 M 2 M D M Δ M + η d 3 ( 1 ε 2 ) S α 2 M B d 4 M + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) D W Δ W + ρ 1 S W d 5 W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) 1 B 2 B D B Δ B + ρ 2 M B d 6 B .
At the equilibrium state, E 2 fulfills the following conditions
μ = d ˜ U 2 + β ( 1 ε 1 ) U 2 M 2 , β ( 1 ε 1 ) U 2 M 2 = γ 1 ( 1 b ) + d 1 N 2 = γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 ( 1 b ) T 2 = d 3 γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) S 2 , β ( 1 ε 1 ) U 2 M 2 = α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) M 2 B 2 + d 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) M 2 .
After using (28), Equation (27) is transformed into
Θ ˜ 2 t = 1 U 2 U d ˜ U 2 d ˜ U + β ( 1 ε 1 ) U 2 M 2 5 U 2 U N T 2 N 2 T T S 2 T 2 S S M 2 S 2 M U N 2 M U 2 N M 2 + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) S 2 S 4 W + 1 U 2 U D U Δ U + 1 N 2 N D N Δ N + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) 1 T 2 T D T Δ T + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) 1 S 2 S D S Δ S + γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) 1 M 2 M D M Δ M + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) D W Δ W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) 1 B 2 B D B Δ B .
From the equilibrium values, we get
S 2 S 4 = β μ γ 1 γ 2 d 6 ( 1 b ) ( 1 ε 1 ) d 3 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d ˜ ρ 2 + β d 6 ( 1 ε 1 ) d 5 ρ 1 = d 5 ρ 1 R l 1 .
After using (19) and (29), the time derivative of Θ 2 ( t ) is given by
d Θ 2 d t = d ˜ Ω U U 2 2 U d x + β ( 1 ε 1 ) U 2 M 2 Ω 5 U 2 U N T 2 N 2 T T S 2 T 2 S S M 2 S 2 M U N 2 M U 2 N M 2 d x + α 1 d 5 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) R l 1 Ω W d x D U U 2 Ω U 2 U 2 d x D N N 2 Ω N 2 N 2 d x D T T 2 γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) Ω T 2 T 2 d x D S S 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) Ω S 2 S 2 d x D M M 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ( 1 b ) ( 1 ε 2 ) Ω M 2 M 2 d x D B B 2 α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 η γ 1 γ 2 ρ 2 ( 1 b ) ( 1 ε 2 ) Ω B 2 B 2 d x .
We see that d Θ 2 d t 0 if R l 1 , and d Θ 2 d t = 0 if U = U 2 , N = N 2 , T = T 2 , S = S 2 , M = M 2 , W = 0 , and B = B 2 . Thus, the singleton { E 2 } is the largest invariant subset of { ( U , N , T , S , M , W , B ) | d Θ 2 d t = 0 } . 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
λ ρ 1 S 2 + d 5 + D W ν i f 2 ( λ ) = 0 ,
where
f 2 ( λ ) = λ + γ 1 ( 1 b ) + d 1 + D N ν i λ + γ 2 + d 2 + D T ν i λ + d 3 + D S ν i λ + β ( 1 ε 1 ) M 2 + d ˜ + D U ν i × α 2 ρ 2 M 2 B 2 + λ + α 2 B 2 + d 4 + D M ν i λ ρ 2 M 2 + d 6 + D B ν i β η γ 1 γ 2 d 3 U 2 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) λ + d ˜ + D U ν i λ ρ 2 M 2 + d 6 + D B ν i .
One root of Equation (30) is given by λ = ρ 1 S 2 d 5 D W ν i . When i = 1 , we get
λ | i = 1 = ρ 1 S 2 S 4 = d 5 R l 1 > 0 if R l > 1 .
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 R j R l 1 . It is unstable if R j R l > 1 .
Proof. 
Take the following Lyapunov functional
Θ 3 ( t ) = Ω Θ ˜ 3 ( x , t ) d x ,
where
Θ ˜ 3 ( x , t ) = U 3 U U 3 1 ln U U 3 + N 3 N N 3 1 ln N N 3 + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) T 3 T T 3 1 ln T T 3 + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) S 3 S S 3 1 ln S S 3 + γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 3 η γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 2 ) M 3 M M 3 1 ln M M 3 + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) W 3 W W 3 1 ln W W 3 + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 3 η γ 1 γ 2 ρ 2 d 3 ( 1 b ) ( 1 ε 2 ) B .
Thus, we get
Θ ˜ 3 t = 1 U 3 U D U Δ U + μ d ˜ U β ( 1 ε 1 ) U M + 1 N 3 N D N Δ N + β ( 1 ε 1 ) U M γ 1 ( 1 b ) N d 1 N + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) 1 T 3 T D T Δ T + γ 1 ( 1 b ) N γ 2 T d 2 T + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) 1 S 3 S D S Δ S + γ 2 T α 1 S W d 3 S + γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 3 η γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 2 ) 1 M 3 M D M Δ M + η d 3 ( 1 ε 2 ) S α 2 M B d 4 M + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) 1 W 3 W D W Δ W + ρ 1 S W d 5 W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 3 η γ 1 γ 2 ρ 2 d 3 ( 1 b ) ( 1 ε 2 ) D B Δ B + ρ 2 M B d 6 B .
At the equilibrium state, E 3 fulfills the following conditions
μ = d ˜ U 3 + β ( 1 ε 1 ) U 3 M 3 , β ( 1 ε 1 ) U 3 M 3 = γ 1 ( 1 b ) + d 1 N 3 = γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 ( 1 b ) T 3 = γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 3 γ 1 γ 2 ( 1 b ) S 3 , β ( 1 ε 1 ) U 3 M 3 = d 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 3 η γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 2 ) M 3 .
After using (32), Equation (31) is simplified to
Θ ˜ 3 t = 1 U 3 U d ˜ U 3 d ˜ U + β ( 1 ε 1 ) U 3 M 3 5 U 3 U N T 3 N 3 T T S 3 T 3 S S M 3 S 3 M U N 3 M U 3 N M 3 + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 3 η γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 2 ) M 3 M 4 B + 1 U 3 U D U Δ U + 1 N 3 N D N Δ N + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) 1 T 3 T D T Δ T + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) 1 S 3 S D S Δ S + γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 3 η γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 2 ) 1 M 3 M D M Δ M + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) 1 W 3 W D W Δ W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 3 η γ 1 γ 2 ρ 2 d 3 ( 1 b ) ( 1 ε 2 ) D B Δ B .
By applying (19), the time derivative of Θ 3 ( t ) is given by
d Θ 3 d t = d ˜ Ω U U 3 2 U d x + β ( 1 ε 1 ) U 3 M 3 Ω 5 U 3 U N T 3 N 3 T T S 3 T 3 S S M 3 S 3 M U N 3 M U 3 N M 3 d x + α 2 d 6 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 3 η γ 1 γ 2 ρ 2 d 3 ( 1 b ) ( 1 ε 2 ) R j R l 1 Ω B d x D U U 3 Ω U 2 U 2 d x D N N 3 Ω N 2 N 2 d x D T T 3 γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) Ω T 2 T 2 d x D S S 3 γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) Ω S 2 S 2 d x D M M 3 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 3 η γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 2 ) Ω M 2 M 2 d x D W W 3 α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) Ω W 2 W 2 d x .
We see that d Θ 3 d t 0 if R j R l 1 , and d Θ 3 d t = 0 at E 3 . Hence, the largest invariant subset of { ( U , N , T , S , M , W , B ) | d Θ 3 d t = 0 } is the singleton { E 3 } . We deduce from LaSalle’s invariance principle [42] that E 3 is globally asymptotically stable if R j R l 1 , given that it is defined if R k > 1 .
The characteristic equation at E 3 is computed as follows
λ ρ 2 M 3 + d 6 + D B ν i f 3 ( λ ) = 0 ,
where
f 3 ( λ ) = λ + γ 1 ( 1 b ) + d 1 + D N ν i λ + γ 2 + d 2 + D T ν i λ + d 4 + D M ν i λ + β ( 1 ε 1 ) M 3 + d ˜ + D U ν i × α 1 ρ 1 S 3 W 3 + λ ρ 1 S 3 + d 5 + D W ν i λ + α 1 W 3 + d 3 + D S ν i β η γ 1 γ 2 d 3 U 3 ( 1 b ) ( 1 ε 1 ) ( 1 ε 2 ) λ + d ˜ + D U ν i λ ρ 1 S 3 + d 5 + D W ν i .
One root of Equation (33) is given by λ = ρ 2 M 3 d 6 D B ν i . Consider the case when i = 1 , we get
λ | i = 1 = ρ 2 M 3 M 4 = d 6 R j R l 1 > 0 if R j R l > 1 .
This implies that Equation (33) has a positive root and E 3 becomes unstable when R j R l > 1 . □
Theorem 7.
The CTL/antibody coexistence equilibrium E 4 is globally asymptotically stable if R l > 1 and R j R l > 1 .
Proof. 
Define the following Lyapunov functional
Θ 4 ( t ) = Ω Θ ˜ 4 ( x , t ) d x ,
where
Θ ˜ 4 ( x , t ) = U 4 U U 4 1 ln U U 4 + N 4 N N 4 1 ln N N 4 + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) T 4 T T 4 1 ln T T 4 + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) S 4 S S 4 1 ln S S 4 + γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 4 η γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 2 ) M 4 M M 4 1 ln M M 4 + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) W 4 W W 4 1 ln W W 4 + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 4 η γ 1 γ 2 ρ 2 d 3 ( 1 b ) ( 1 ε 2 ) B 4 B B 4 1 ln B B 4 .
By taking the time partial derivative of Θ ˜ 4 ( x , t ) , we get
Θ ˜ 4 t = 1 U 4 U D U Δ U + μ d ˜ U β ( 1 ε 1 ) U M + 1 N 4 N D N Δ N + β ( 1 ε 1 ) U M γ 1 ( 1 b ) N d 1 N + γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) 1 T 4 T D T Δ T + γ 1 ( 1 b ) N γ 2 T d 2 T + γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) 1 S 4 S D S Δ S + γ 2 T α 1 S W d 3 S + γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 4 η γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 2 ) 1 M 4 M D M Δ M + η d 3 ( 1 ε 2 ) S α 2 M B d 4 M + α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) 1 W 4 W D W Δ W + ρ 1 S W d 5 W + α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 4 η γ 1 γ 2 ρ 2 d 3 ( 1 b ) ( 1 ε 2 ) 1 B 4 B D B Δ B + ρ 2 M B d 6 B .
At the equilibrium state, E 4 satisfies the following conditions
μ = d ˜ U 4 + β ( 1 ε 1 ) U 4 M 4 , β ( 1 ε 1 ) U 4 M 4 = γ 1 ( 1 b ) + d 1 N 4 = γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 ( 1 b ) T 4 = γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 4 γ 1 γ 2 ( 1 b ) S 4 , β ( 1 ε 1 ) U 4 M 4 = α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 4 η γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 2 ) M 4 B 4 + d 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 4 η γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 2 ) M 4 .
By using (19) and (34), the time derivative of Θ 4 ( t ) is given by
d Θ 4 d t = d ˜ Ω U U 4 2 U d x + β ( 1 ε 1 ) U 4 M 4 Ω 5 U 4 U N T 4 N 4 T T S 4 T 4 S S M 4 S 4 M U N 4 M U 4 N M 4 d x D U U 4 Ω U 2 U 2 d x D N N 4 Ω N 2 N 2 d x D T T 4 γ 1 ( 1 b ) + d 1 γ 1 ( 1 b ) Ω T 2 T 2 d x D S S 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 γ 1 γ 2 ( 1 b ) Ω S 2 S 2 d x D M M 4 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 4 η γ 1 γ 2 d 3 ( 1 b ) ( 1 ε 2 ) Ω M 2 M 2 d x D W W 4 α 1 γ 1 ( 1 b ) + d 1 γ 2 + d 2 ρ 1 γ 1 γ 2 ( 1 b ) Ω W 2 W 2 d x D B B 4 α 2 γ 1 ( 1 b ) + d 1 γ 2 + d 2 d 3 + α 1 W 4 η γ 1 γ 2 ρ 2 d 3 ( 1 b ) ( 1 ε 2 ) Ω B 2 B 2 d x .
It is easy to see that d Θ 4 d t 0 and d Θ 4 d t = 0 at E 4 . Thus, the largest invariant subset of { ( U , N , T , S , M , W , B ) | d Θ 4 d t = 0 } is the singleton { E 4 } . It is deduced from LaSalle’s invariance principle [42] that E 4 is globally asymptotically stable if R l > 1 and R j R l > 1 . □

5. Numerical Simulations

In this section, we conduct some numerical simulations to illustrate the results obtained in the preceding sections. We select the spatial domain as Ω = [ 0 , 2 ] . The space and time step sizes are selected as Δ x = 0.02 and Δ t = 0.1 , respectively. The initial conditions of model (3) are considered to be the following
U ( x , 0 ) = 0.5 × 10 10 ( 1 + 0.5 cos 2 ( π x ) ) cells mL 1 , N ( x , 0 ) = 0 cells mL 1 , T ( x , 0 ) = 0 cells mL 1 , S ( x , 0 ) = 0 cells mL 1 , M ( x , 0 ) = 10 6 ( 1 + 0.5 cos 2 ( π x ) ) cells mL 1 , W ( x , 0 ) = 0.0001 ( 1 + 0.5 cos 2 ( π x ) ) cells mL 1 , B ( x , 0 ) = 0.0001 ( 1 + 0.5 cos 2 ( π x ) ) cells mL 1 , x [ 0 , 2 ] .
These values are mainly based on the values given in [26,27]. For numerical simulations, we vary the values of ε 1 , b, ε 2 , β , ρ 1 , and ρ 2 . The values of the rest parameters are fixed, and they are listed in Table 1. The diffusion coefficients values are taken as assumptions, while the values of the other parameters are based on values considered in the literature.

5.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.
Case  2.
We select the values β = 2 × 10 8 , ρ 1 = 2 × 10 9 , and ρ 2 = 3 × 10 9 . The threshold conditions are given by R i = 1.3729 > 1 , R j = 0.1056 < 1 , and R k = 0.371 < 1 . With these values, the equilibrium E 1 = ( 7.284 × 10 9 , 1.151 × 10 8 , 1.726 × 10 7 , 3.453 × 10 6 , 5.179 × 10 5 , 0 , 0 ) is globally asymptotically stable (see Figure 2). This result coincides with the result of Theorem 4. When the CTL and antibody immune responses are not active, the malaria parasites succeed in infecting red blood cells and establishing the infection.
Case  3.
We take the values β = 2 × 10 7 , ρ 1 = 2 × 10 9 , and ρ 2 = 3 × 10 7 . The corresponding thresholds are R j = 6.2404 > 1 and R l = 0.2496 < 1 . It follows that the equilibrium E 2 = ( 4.546 × 10 9 , 2.308 × 10 8 , 3.463 × 10 7 , 6.925 × 10 6 , 1.665 × 10 5 , 0 , 2.515 × 10 10 ) is globally asymptotically stable (see Figure 3), which supports Theorem 5. In these circumstances, the malaria infection stimulates the production of antibodies to clear free merozoites from the blood. Nevertheless, the malaria infection persists despite the presence of both antibody immune response and blood-stage drugs.
Case  4.
We select the values β = 2 × 10 7 , ρ 1 = 2 × 10 8 , and ρ 2 = 3 × 10 9 . This gives R k = 3.7107 > 1 and R j R l = 0.0225 < 1 . In agreement with Theorem 6, the equilibrium E 3 = ( 2.703 × 10 9 , 3.092 × 10 8 , 4.638 × 10 7 , 2.5 × 10 6 , 3.75 × 10 5 , 1.355 × 10 8 , 0 ) is globally asymptotically stable (see Figure 4). Here, the malaria infection activates the CTL immune response to kill the schizont infected red blood cells, which decreases the count of free merozoites in the blood.
Case  5.
We take β = 2 × 10 7 , ρ 1 = 2 × 10 8 , and ρ 2 = 3 × 10 7 . This gives R l = 2.4961 > 1 and R j R l = 2.25 > 1 . Thus, the equilibrium E 4 = ( 4.545 × 10 9 , 2.311 × 10 8 , 3.467 × 10 7 , 2.5 × 10 6 , 1.667 × 10 5 , 8.867 × 10 7 , 6 × 10 9 ) is globally asymptotically stable (see Figure 5), which agrees with Theorem 7. At this point, the malaria infection stimulates both the antibody immune response and CTL immune response. These responses work together to increase the concentration of uninfected cells and reduce the concentration of free merozoites and infected cells at all stages. With low treatment efficacy, the immune responses are not able to eliminate the infection.
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 Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5.

5.2. 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.

6. 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.

Author Contributions

Conceptualization, A.A.A.; Formal analysis, A.E.; Investigation, A.A.A.; Methodology, A.E.; Supervision, A.E. All authors have read and agreed to the published version of the manuscript.

Funding

This project was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, under grant No. (D-662-130-1441). The authors, therefore, gratefully acknowledge the DSR technical and financial support.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. WHO. World Malaria Report 2019; License: CC BY-NC-SA 3.0 IGO; World Health Organization: Geneva, Switzerland, 2019. [Google Scholar]
  2. Chen, H.; Wang, W.; Fu, R.; Luo, J. Global analysis of a mathematical model on malaria with competitive strains and immune responses. Appl. Math. Comput. 2015, 259, 132–152. [Google Scholar] [CrossRef]
  3. Khoury, D.; Aogo, R.; Randriafanomezantsoa-Radohery, G.; McCaw, J.; Simpson, J.; McCarthy, J.; Haque, A.; Cromer, D.; Davenport, M. Within-host modeling of blood-stage malaria. Immun. Rev. 2018, 285, 168–193. [Google Scholar] [CrossRef]
  4. Niger, A.; Gumel, A. Immune response and imperfect vaccine in malaria dynamics. Math. Popul. Stud. 2011, 18, 55–86. [Google Scholar] [CrossRef]
  5. Agusto, F.; Leite, M.; Orive, M. The transmission dynamics of a within-and between-hosts malaria model. Ecol. Complex. 2019, 38, 31–55. [Google Scholar] [CrossRef]
  6. Song, T.; Wang, C.; Tian, B. Mathematical models for within-host competition of malaria parasites. Math. Biosci. Eng. 2019, 16, 6623–6653. [Google Scholar] [CrossRef] [PubMed]
  7. Zaloumis, S.; Humberstone, A.; Charman, S.; Price, R.; Moehrle, J.; Gamo-Benito, J.; McCaw, J.; Jamsen, K.; Smith, K.; Simpson, J. Assessing the utility of an anti-malarial pharmacokinetic-pharmacodynamic model for aiding drug clinical development. Malaria J. 2012, 11, 303. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  8. Nowak, M.A.; May, R.M. Virus Dynamics: Mathematical Principles of Immunology and Virology; Oxford University: Oxford, UK, 2000. [Google Scholar]
  9. Ren, X.; Tian, Y.; Liu, L.; Liu, X. A reaction-diffusion within-host HIV model with cell-to-cell transmission. J. Math. Biol. 2018, 76, 1831–1872. [Google Scholar] [CrossRef] [PubMed]
  10. Elaiw, A.M.; Hobiny, A.D.; Al Agha, A.D. Global dynamics of reaction-diffusion oncolytic M1 virotherapy with immune response. Appl. Math. Comput. 2020, 367, 1–21. [Google Scholar] [CrossRef]
  11. Elaiw, A.M.; AlShamrani, N.H. Stability of an adaptive immunity pathogen dynamics model with latency and multiple delays. Math. Meth. Appl. Sci. 2018, 41, 6645–6672. [Google Scholar] [CrossRef]
  12. Miao, H.; Teng, Z.; Abdurahman, X.; Li, Z. Global stability of a diffusive and delayed virus infection model with general incidence function and adaptive immune response. Comput. Appl. Math. 2017, 37, 3780–3805. [Google Scholar] [CrossRef]
  13. Elaiw, A.M.; Elnahary, E.K. Analysis of general humoral immunity HIV dynamics model with HAART and distributed delays. Mathematics 2019, 7, 157. [Google Scholar] [CrossRef] [Green Version]
  14. Elaiw, A.M.; Alshehaiween, S.F.; Hobiny, A.D. Global properties of delay-distributed HIV dynamics model including impairment of B-cell functions. Mathematics 2019, 7, 837. [Google Scholar] [CrossRef] [Green Version]
  15. Elaiw, A.M.; Alshehaiween, S.F. Global stability of delay-distributed viral infection model with two modes of viral transmission and B-cell impairment. Math. Methods Appl. Sci. 2020. [Google Scholar] [CrossRef]
  16. Elaiw, A.M.; AlShamrani, N.H. Stability of a general adaptive immunity virus dynamics model with multi-stages of infected cells and two routes of infection. Math. Methods Appl. Sci. 2020, 43, 1145–1175. [Google Scholar] [CrossRef]
  17. Anderson, R.; May, R.; Gupta, S. Non-linear phenomena in host-parasite interactions. Parasitology 1989, 99, S59–S79. [Google Scholar] [CrossRef]
  18. Anderson, R. Complex dynamic behaviours in the interaction between parasite populations and the host’s immune system. Int. J. Parasitol. 1998, 28, 551–566. [Google Scholar] [CrossRef]
  19. Saul, A. Models for the in-host dynamics of malaria revisited: Errors in some basic models lead to large over-estimates of growth rates. Parasitology 1998, 117, 405–407. [Google Scholar] [CrossRef]
  20. Gravenor, M.; Lloyd, A. Reply to: Models for the in-host dynamics of malaria revisited: Errors in some basic models lead to large over-estimates of growth rates. Parasitology 1998, 171, 409–410. [Google Scholar] [CrossRef] [Green Version]
  21. Hoshen, M.; Heinrich, R.; Stein, W.; Ginsburg, H. Mathematical modelling of the within-host dynamics of Plasmodium falciparum. Parasitology 2000, 121, 227–235. [Google Scholar] [CrossRef]
  22. Iggidr, A.; Kamgang, J.; Sallet, G.; Tewa, J. Global analysis of new malaria intrahost models with a competitive exclusion principle. SIAM J. Appl. Math. 2006, 67, 260–278. [Google Scholar] [CrossRef] [Green Version]
  23. Saralamba, S.; Pan-Ngum, W.; Maude, R.; Lee, S.J.; Tarning, J.; Lindegårdh, N.; Chotivanich, K.; Nosten, F.; Day, N.; Socheat, D.; et al. Intrahost modeling of artemisinin resistance in Plasmodium falciparum. Proc. Natl. Acad. Sci. USA 2011, 108, 397–402. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Demasse, R.; Ducrot, A. An age-structured within-host model for multistrain malaria infections. SIAM J. Appl. Math. 2013, 73, 572–593. [Google Scholar] [CrossRef] [Green Version]
  25. Li, Y.; Ruan, S.; Xiao, D. The within-host dynamics of malaria infection with immune response. Math. Biosci. Eng. 2011, 8, 999–1018. [Google Scholar] [PubMed]
  26. Tumwiine, J.; Mugisha, J.; Luboobi, L. On global stability of the intra-host dynamics of malaria and the immune system. J. Math. Anal. Appl. 2008, 341, 855–869. [Google Scholar] [CrossRef] [Green Version]
  27. Hetzel, C.; Anderson, R. The within-host cellular dynamics of bloodstage malaria: Theoretical and experimental studies. Parasitology 1996, 113, 25–38. [Google Scholar] [CrossRef]
  28. Tewa, J.; Fokouop, R.; Mewoli, B.; Bowong, S. Mathematical analysis of a general class of ordinary differential equations coming from within-hosts models of malaria with immune effectors. App. Math. Comput. 2012, 2018, 7347–7361. [Google Scholar] [CrossRef]
  29. Tchinda, P.; Tewa, J.; Mewoli, B.; Bowong, S. Mathematical analysis of a general class of intra-host model of malaria with “allee effect“. J. Nonl. Syst. Appl. 2013, 36, 52. [Google Scholar]
  30. Cai, L.; Tuncer, N.; Martcheva, M. How does within-host dynamics affect population-level dynamics? Insights from an immuno-epidemiological model of malaria. Math. Methods Appl. Sci. 2017, 40, 6424–6450. [Google Scholar] [CrossRef]
  31. Orwa, T.; Mbogo, R.; Luboobi, L. Mathematical model for the in-host malaria dynamics subject to malaria vaccines. Lett. Biomath. 2018, 5, 222–251. [Google Scholar] [CrossRef] [Green Version]
  32. Funk, G.; .Jansen, V.A.A.; Bonhoeffer, S.; Killingback, T. Spatial models of virus-immune dynamics. J. Theor. Biol. 2005, 233, 221–236. [Google Scholar] [CrossRef]
  33. Takoutsing, E.; Temgoua, A.; Yemele, D.; Bowong, S. Dynamics of an intra-host model of malaria with periodic antimalarial treatment. Int. J. Nonl. Sci. 2019, 27, 148–164. [Google Scholar]
  34. Babbitt, S.; Altenhofen, L.; Cobbold, S.; Istvan, E.; Fennell, C.; Doerig, C.; Llinas, M.; Goldberg, D. Plasmodium falciparum responds to amino acid starvation by entering into a hibernatory state. Proc. Natl. Acad. Sci. USA 2012, 109, E3278–E3287. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Zhang, Y.; Xu, Z. Dynamics of a diffusive HBV model with delayed Beddington-DeAngelis response. Nonl. Anal. Real World Appl. 2014, 15, 118–139. [Google Scholar] [CrossRef]
  36. Xu, Z.; Xu, Y. Stability of a CD4+ T cell viral infection model with diffusion. Int. J. Biomath. 2018, 11, 1–16. [Google Scholar] [CrossRef]
  37. Smith, H.L. Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems; American Mathematical Society: Providence, RI, USA, 1995. [Google Scholar]
  38. Protter, M.H.; Weinberger, H.F. Maximum Principles in Differential Equations; Prentic Hall: Englewood Cliffs, NJ, USA, 1967. [Google Scholar]
  39. Henry, D. Geometric Theory of Semilinear Parabolic Equations; Springer: New York, NY, USA, 1993. [Google Scholar]
  40. Hattaf, K.; Yousfi, N. Global stability for reaction-diffusion equations in biology. Comput. Math. Appl. 2013, 66, 1488–1497. [Google Scholar] [CrossRef]
  41. Hsu, S. A survey of constructing lyapunov functions for mathematical models in population biology. Taiwan. J. Math. 2005, 9, 151–173. [Google Scholar] [CrossRef]
  42. Hale, J.K.; Verduyn Lunel, S.M. Introduction to Functional Differential Equations; Springer: New York, NY, USA, 1993. [Google Scholar]
  43. Khalil, H.K. Nonlinear Systems; Prentice-Hall: New Jersey, NJ, USA, 1996. [Google Scholar]
  44. Bakker, N.; Eling, W.; De Groot, A.; Sinkeldam, E.; Luyken, R. Attenuation of malaria infection, paralysis and lesions in the central nervous system by low protein diets in rats. Acta Trop. 1992, 50, 285–293. [Google Scholar] [CrossRef]
  45. Edirisinghe, J.; Fern, E.; Targett, G. Resistance to superinfection with Plasmodium berghei in rats fed a protein-free diet. Trans. R Soc. Trop. Med. Hyg. 1982, 76, 382–386. [Google Scholar] [CrossRef]
  46. Sun, H.; Wang, J. Dynamics of a diffusive virus model with general incidence function, cell-to-cell transmission and time delay. Comput. Math. Appl. 2019, 77, 284–301. [Google Scholar] [CrossRef]
  47. Elaiw, A.M.; Raezah, A.A. Stability of general virus dynamics models with both cellular and viral infections and delays. Math. Methods Appl. Sci. 2017, 40, 5863–5880. [Google Scholar] [CrossRef]
  48. Elaiw, A.M.; AlShamrani, N.H. Global stability of a delayed adaptive immunity viral infection with two routes of infection and multi-stages of infected cells. Commun. Nonl. Sci. Num. Simul. 2020, 86, 1–27. [Google Scholar] [CrossRef]
  49. Elaiw, A.M.; AlAgha, A.D. Analysis of a Delayed and Diffusive Oncolytic M1 Virotherapy Model with Immune Response; Nonlinear Analysis: Real World Applications; Article 103116; Elsevier: Amsterdam, The Netherlands, 2020. [Google Scholar] [CrossRef]
  50. Luo, J.; Wang, W.; Chen, H.; Fu, R. Bifurcations of a mathematical model for HIV dynamics. J. Math. Anal. Appl. 2016, 434, 837–857. [Google Scholar] [CrossRef]
Figure 1. The numerical simulations of model (3) when R i 1 . The disease-free equilibrium E 0 is globally asymptotically stable. The subfigures show the spatiotemporal behaviors of (a) uninfected red blood cells U, (b) ring infected red blood cells N, (c) trophozoite infected red blood cells T, (d) schizont infected red blood cells S, (e) free merozoites M, (f) CTLs W, and (g) antibodies B.
Figure 1. The numerical simulations of model (3) when R i 1 . The disease-free equilibrium E 0 is globally asymptotically stable. The subfigures show the spatiotemporal behaviors of (a) uninfected red blood cells U, (b) ring infected red blood cells N, (c) trophozoite infected red blood cells T, (d) schizont infected red blood cells S, (e) free merozoites M, (f) CTLs W, and (g) antibodies B.
Mathematics 08 00563 g001
Figure 2. The numerical simulations of model (3) when R i > 1 , R j 1 , and R k 1 . The immune-free equilibrium E 1 is globally asymptotically stable. The subfigures show the spatiotemporal behaviors of (a) uninfected red blood cells U, (b) ring infected red blood cells N, (c) trophozoite infected red blood cells T, (d) schizont infected red blood cells S, (e) free merozoites M, (f) CTLs W, and (g) antibodies B.
Figure 2. The numerical simulations of model (3) when R i > 1 , R j 1 , and R k 1 . The immune-free equilibrium E 1 is globally asymptotically stable. The subfigures show the spatiotemporal behaviors of (a) uninfected red blood cells U, (b) ring infected red blood cells N, (c) trophozoite infected red blood cells T, (d) schizont infected red blood cells S, (e) free merozoites M, (f) CTLs W, and (g) antibodies B.
Mathematics 08 00563 g002
Figure 3. The numerical simulations of model (3) when R j > 1 and R l 1 The antibody-activated equilibrium E 2 is globally asymptotically stable. The subfigures show the spatiotemporal behaviors of (a) uninfected red blood cells U, (b) ring infected red blood cells N, (c) trophozoite infected red blood cells T, (d) schizont infected red blood cells S, (e) free merozoites M, (f) CTLs W, and (g) antibodies B.
Figure 3. The numerical simulations of model (3) when R j > 1 and R l 1 The antibody-activated equilibrium E 2 is globally asymptotically stable. The subfigures show the spatiotemporal behaviors of (a) uninfected red blood cells U, (b) ring infected red blood cells N, (c) trophozoite infected red blood cells T, (d) schizont infected red blood cells S, (e) free merozoites M, (f) CTLs W, and (g) antibodies B.
Mathematics 08 00563 g003
Figure 4. The numerical simulations of model (3) when R k > 1 and R j R l 1 The CTL-activated equilibrium E 3 is globally asymptotically stable. The subfigures show the spatiotemporal behaviors of (a) uninfected red blood cells U, (b) ring infected red blood cells N, (c) trophozoite infected red blood cells T, (d) schizont infected red blood cells S, (e) free merozoites M, (f) CTLs W, and (g) antibodies B.
Figure 4. The numerical simulations of model (3) when R k > 1 and R j R l 1 The CTL-activated equilibrium E 3 is globally asymptotically stable. The subfigures show the spatiotemporal behaviors of (a) uninfected red blood cells U, (b) ring infected red blood cells N, (c) trophozoite infected red blood cells T, (d) schizont infected red blood cells S, (e) free merozoites M, (f) CTLs W, and (g) antibodies B.
Mathematics 08 00563 g004
Figure 5. The numerical simulations of model (3) when R l > 1 and R j R l > 1 The coexistence equilibrium E 4 is globally asymptotically stable. The subfigures show the spatiotemporal behaviors of (a) uninfected red blood cells U, (b) ring infected red blood cells N, (c) trophozoite infected red blood cells T, (d) schizont infected red blood cells S, (e) free merozoites M, (f) CTLs W, and (g) antibodies B.
Figure 5. The numerical simulations of model (3) when R l > 1 and R j R l > 1 The coexistence equilibrium E 4 is globally asymptotically stable. The subfigures show the spatiotemporal behaviors of (a) uninfected red blood cells U, (b) ring infected red blood cells N, (c) trophozoite infected red blood cells T, (d) schizont infected red blood cells S, (e) free merozoites M, (f) CTLs W, and (g) antibodies B.
Mathematics 08 00563 g005
Figure 6. The effect of varying the efficacy of blood-stage treatment ε 1 and the efficacy of isoleucine starvation b on the within-host basic reproductive number R i 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. (a) R i is a nonlinear function of b, where increasing the value of b causes the disease-free equilibrium E 0 to become stable at approximately b = 0.375 . (b) R i is a linear function of ε 1 where increasing the value of ε 1 causes the disease-free equilibrium E 0 to become stable at approximately ε 1 = 0.344 .
Figure 6. The effect of varying the efficacy of blood-stage treatment ε 1 and the efficacy of isoleucine starvation b on the within-host basic reproductive number R i 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. (a) R i is a nonlinear function of b, where increasing the value of b causes the disease-free equilibrium E 0 to become stable at approximately b = 0.375 . (b) R i is a linear function of ε 1 where increasing the value of ε 1 causes the disease-free equilibrium E 0 to become stable at approximately ε 1 = 0.344 .
Mathematics 08 00563 g006
Figure 7. The effect of changing the diffusion coefficient on the concentration of free merozoites M at the beginning of infection. All parameters are identical to those used in Figure 5, while the diffusion coefficient D M is changed from (a) D M = 0.2 to (b) D M = 0.0002 .
Figure 7. The effect of changing the diffusion coefficient on the concentration of free merozoites M at the beginning of infection. All parameters are identical to those used in Figure 5, while the diffusion coefficient D M is changed from (a) D M = 0.2 to (b) D M = 0.0002 .
Mathematics 08 00563 g007
Table 1. Parameters values of model (3).
Table 1. Parameters values of model (3).
ParameterDescriptionValueReference
μ Recruitment rate of red blood cells 2.5 × 10 8 cell mL 1 day 1  [27]
d ˜ Death rate of uninfected red blood cells 0.025 day 1  [27]
β Infection rate of red blood cellsVaried
γ 1 Progression rate from rings to trophozoites 0.1 day 1  [31]
γ 2 Progression rate from trophozoites to schizonts 0.1 day 1  [31]
η Number of merozoites released per schizont infected cell16 [17]
α 1 Killing rate of infected cells by CTLs 10 8 mL cell 1 day 1  [27]
α 2 Removal rate of free merozoites by antibodies 10 8 mL cell 1 day 1  [27]
ρ 1 Stimulation rate of CTL immune responseVaried
ρ 2 Stimulation rate of antibody immune responseVaried
d 1 Death rate of ring infected red blood cells 0.5 day 1  [31]
d 2 Death rate of trophozoite infected red blood cells 0.5 day 1  [31]
d 3 Death rate of schizont infected red blood cells 0.5 day 1  [31]
d 4 Death rate of free merozoites48 day 1  [27]
d 5 Decay rate of CTLs 0.05 day 1  [17]
d 6 Decay rate of antibodies 0.05 day 1  [17]
ε 1 Efficacy of blood-stage treatmentVaried
bEfficacy of isoleucine starvationVaried
ε 2 Efficacy of blood-stage treatmentVaried
D U Diffusion coefficient of uninfected red blood cells 0.1 Assumed
D N Diffusion coefficient of ring infected red blood cells 0.1 Assumed
D T Diffusion coefficient of trophozoite infected red blood cells 0.01 Assumed
D S Diffusion coefficient of schizont infected red blood cells 0.01 Assumed
D M Diffusion coefficient of free merozoites 0.2 Assumed
D W Diffusion coefficient of CTLs 0.01 Assumed
D B Diffusion coefficient of antibodies 0.2 Assumed

Share and Cite

MDPI and ACS Style

Elaiw, A.; Al Agha, A. Global Analysis of a Reaction-Diffusion Within-Host Malaria Infection Model with Adaptive Immune Response. Mathematics 2020, 8, 563. https://doi.org/10.3390/math8040563

AMA Style

Elaiw A, Al Agha A. Global Analysis of a Reaction-Diffusion Within-Host Malaria Infection Model with Adaptive Immune Response. Mathematics. 2020; 8(4):563. https://doi.org/10.3390/math8040563

Chicago/Turabian Style

Elaiw, Ahmed, and Afnan Al Agha. 2020. "Global Analysis of a Reaction-Diffusion Within-Host Malaria Infection Model with Adaptive Immune Response" Mathematics 8, no. 4: 563. https://doi.org/10.3390/math8040563

APA Style

Elaiw, A., & Al Agha, A. (2020). Global Analysis of a Reaction-Diffusion Within-Host Malaria Infection Model with Adaptive Immune Response. Mathematics, 8(4), 563. https://doi.org/10.3390/math8040563

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop