The blood-stage dynamics of malaria infection with immune response

This article is concerned with the dynamics of malaria infection model with diffusion and delay. The disease free threshold and the immune response threshold value of the malaria infection are obtained, which characterize the stability of the disease free equilibrium and infection equilibrium (with or without immune response). In addition, fluctuations occur when the model undergoes Hopf bifurcation as the delay passes through a certain critical value . In this case, periodic oscillation appears near the positive steady state, which implies the recurrent attacks of disease. Finally, numerical simulations are provided to illustrate the theoretical results.


Introduction
Malaria is one of the top three infectious diseases in the world and poses a great threat to human health. Each year, approximately 2.2 billion people are affected by Plasmodium falciparum malaria worldwide, resulting in approximately 515 million endemic cases [18]. Malaria causes millions of infections every year, and more than 1.1 million people die of the disease, especially infants, young children and pregnant women [20].
Malaria infections are mainly caused by P. falciparum, P. malariae, P. ovale and P. vivax. In particular, P. falciparum is the most virulent, causing almost as many deaths as all other types of malaria infection combined [15]. It is well known that anopheles are the carriers of malaria. When an infected anopheles mosquito bites a host, the malaria parasites first penetrate liver cells of the host and then move into the blood, where they multiply and undergo replication cycles in the red blood cells. After 10-14 days or more, the malaria parasite matures inside the biting mosquitos until they can infect others. As malaria parasites reproduce in their hosts, they can stimulate the production of immune cells, leading to an immune response in the host.
In the past few years, many efforts have been made to control and prevent malaria [19]. However, it remains difficult for those people living in the areas with high malaria prevalence to have access to effective prevention, diagnosis and treatment. Moreover, with the emergence of malaria-resistant strains and insecticide-resistant mosquitoes, finding methods to control malaria infection remains a major challenge. The main harm caused by malaria parasites in the host occurs during the blood stage of development and reproduction. Therefore, in order to elucidate this mechanism, a lot of mathematical models have been used to describe the dynamic process of malaria infection. Up to now, a great deal of works have also been studied on the dynamics of plasmodium in the red blood cells of the host [3,6,11,21,24]. Spatial diffusion of human and vector populations may be important to the malaria dynamics [1,4,12,13,25].
In 2016, Chang et al. [2] made a significant progress in the study of the innate immune mechanism of P. falciparum escaping from the immune clearance in the host. The host's first line of defense against pathogen invasion is the innate immune response. When stimulated by infected red blood cells (iRBC), host's neutrophils release chromatin and lysosomal granules in the cytoplasm, forming a reticular structure (NETs) in a way of active death (Netosis). NETs can capture and kill infected red blood cells. Meanwhile, plasmodium can degrade NETs by secreting DNA enzymes (TatD-like DNase) to escape from host's immune clearance. Chang [2] pointed out that TatD-like DNase was an important factor to the survival of malaria parasites and a potential candidate for malaria vaccine.
In 2019, Ding et al. [5] studied the dynamics of malaria parasites infection in the host blood-stage with time delay. They established the following model: where u 1 , u 2 and u 3 are populations of uninfected red blood cells (RBCs), infected red blood cells (iRBCs) and immune factors (NETs), respectively. is the production rate of RBCs. d 1 , d 2 and d 3 represent the decay rates of RBCs, iRBCs and NETs, respectively. δ, α and β are the infection rate of RBCs by plasmodium, removal rate of iRBCs by NETs, production rate of NETs stimulated by iRBCs. τ denotes the time lag of the reproduction of immune factors. It shows that there are two threshold values 0 and 1 such that the malaria-free equilibrium was globally asymptotically stable if 0 ≤ 1, and if 0 > 1, there exist two kinds of infection equilibria: malaria infection equilibrium (without specific immune response) and positive equilibrium (with specific immune response). Based on the threshold values 0 and 1 , they obtained the existence and stability of aforementioned equilibria. Moreover, the results show that time delay would destroy the stability of the equilibrium state, causing Hopf bifurcation and periodic oscillation near the equilibrium state. According to etiological literature, when an infected anopheles mosquito bites a host, the malaria parasite first penetrates the liver and then enters the bloodstream and invades red blood cells. They ingest nutrients, reproduce and develop in infected red blood cells, and then transfer and invade other normal red blood cells. When the concentration of infected red blood cells in the blood reaches a certain amount, the host will show various symptoms of discomfort. We note that both malaria parasites and red blood cells move into the bloodstream from one location to another, and it makes sense to analyse malaria's detrimental to the host by considering the rate at which the parasite moves from one infected red blood cell to another normal one.
Based on [2] and model (1), in this paper we will consider the following malaria infection in vivo model with dispersal and delay.
The rest of this paper is organized as follows. In Section 2, we present global existence and uniqueness of solutions of initial boundary value problem (2). The threshold properties and stability analysis are established in Section 3. Section 4 is devoted to exhibiting periodic oscillations due to the Hopf bifurcation at the positive equilibrium E * . In Section 5, numerical simulations are provided to demonstrate the theoretical results. A short discussion on research results is given in the last section.
To proceed further, we need some information on the following scalar reaction-diffusion equation, where D, d are positive constants, g(x) is continuous and positive function on¯ .

Lemma 2.2:
System (6) admits a unique steady state ω * (x), which is globally attractive in The proof is similar to that of in [13] (Lemma 1) with some minor modifications, so we omit here.

Threshold dynamics
In this section, we introduce the threshold 0 for the time-delayed and diffusion model (2). To do so, we first study the infection-free steady state E 0 . As we know, an equilibrium E 0 is independent of t and the densities of the diseased compartments u 2 be zero, we consider the following problem: By Lemma 2.2, it is obvious that (3.1) admits a positive steady state u 0 1 , which is globally asymptotically stable. Consider the linearizing system of (2) at the infection-free equilibrium E 0 = ( d 1 , 0, 0), we define two threshold values 0 and 1 as follows. For model (2), the basic reproduction number of iRBCs, describing the average number of newly infected cells generated by one infected cell during its lifespan, is defined by and the immune response threshold value is defined as below The following results can be verified directly.
Next, we study the stability of the equilibrium E 0 , E 1 and E * .
Proof: Define a functional L 0 : C + τ → R as follows: Let Then By Green formula and Neumann boundary conditions, we obtain Similarly, Using the Divergence Theorem and Neumann boundary conditions, we have Obviously, if 0 ≤ 1, then g 2 (x, t) ≤ 0, so d dt L 0 (u 1t , u 2t , u 3t ) ≤ 0. In addition, the largest invariant set of system (2) on the region 0 : According to LaSalles invariance principle [10] (Theorem 4.3.4), it follows that equilibrium E 0 is globally asymptotically stable when 0 ≤ 1. This completes the proof. (2) is globally asymptotically stable.
Proof: Define a functional L (1) : C + → R as follows: We calculate derivatives of Then By Green formula and Neumann boundary conditions, we obtain that Using the Divergence Theorem and the Neumann boundary conditions, we have On the other hand, since . It can be verified that the equilibrium {E 1 } is the only largest invariant set in 1 Therefore, by LaSalle's invariance principle, we conclude that the endemic equilibrium E 1 is globally asymptotically stable for 1 ≥ 0 > 1.

Proof: Let
Differentiating V * with respect to t along model (2), we get Let Then By a similar argument as in the proof of Theorem 3.2, we have m 1 (x, t) ≤ 0 and m 2 (x, t) ≤ 0. Therefore, dV * dt ≤ 0 for all u i > 0, i = 1, 2, 3. Again by LaSalle's invariance principle, we conclude that if 0 > 1 , then the endemic equilibrium E * is globally asymptotically stable in X + .
In the following, we will study local stability of the endemic equilibrium E * when τ > 0. The linearization of system (2) at E * reads as where for any ϕ ∈ C τ , L E * (ϕ) is defined as The characteristic equation of (10) is as follows: It is well known that the eigenvalue problem − φ = μφ, φ x | x=0,lπ = 0 has eigenvalues μ k = k 2 l 2 , k ∈ N {0} with corresponding eigenfunctions cos kx l . By using Fourier expansion φ(x) = ∞ k=0 ccca k b k c k cos kx l , a k , b k , c k ∈ R, Equation (12) is equivalent to where and I is the identity matrix. By a direct computation, (13) is reduced to where Now we consider the stability of the equilibrium E * if 0 > 1 and τ > 0. By Theorem 4.4 in [17], for sufficiently small delay, the characteristic roots of (12) are either very near the eigenvalues of A + B or having more negative real parts than any of the eigenvalues of A + B. Hence, when the delay is small enough, the equilibrium E * is locally asymptotically stable.
Since (12). Note that any complex roots to the equations (13) appear in pairs, and all roots have negative real parts if τ = 0. Therefore, any root of (13) has negative real part for sufficiently small τ . Assume that there exists τ = τ 0 , such that (14) has a pair of pure imaginary roots for some k ∈ N ∪ {0}, denoted by λ = ±ω k i(ω k > 0). Substituting λ = ω k i into (14), we have

Separating the real and imaginary parts leads to
It follows that Let z k = ω 2 k . Then where Equation (15) has no positive roots if it satisfies p k > 0, r k > 0 and p k q k − r k > 0 (Hurwitz criterion). Hence, the threshold value τ 0 doesn't exist. It follows that all solutions of Equation (14) have negative real parts and E * is locally asymptotically stable for any τ > 0. Therefore, we have the following theorem.
Proof: Recall that r k = C 2 k − C 2 k1 . Unfortunately, we can't guarantee that r k is always greater than 0 (see Figure 1). Note that C k + C k1 = D 3 is not guaranteed to be greater than 0. Fortunately, by the expression of C k − C k1 and D k = Dk 2 l 2 , it is clear that for any given parameter, there exists a finite number of k such that C k − C k1 < 0. Hence, r k < 0 for some k ∈ N ∪ {0}. Theorem 3.6: If 0 > 1 and for some k ∈ N ∪ {0}, r k < 0, then there exists τ 0 > 0, such that Equation (14) admits a pair of conjugate purely imaginary roots ±ωi when τ = τ 0 . Moreover, the infection equilibrium E * is locally asymptotically stable if τ < τ 0 .
Proof: We know that if r k < 0, then Equation (15) has at least one positive root and at most three positive roots. Assume that Equation (15) has three positive roots, denoted by On the left-hand side, for the sake of convenience, we take β = 7 × 10 −6 , D = l, such that R 0 = 150, R 1 = 26.8176. The variable relationship between r k and k is considered, where k = 0, r k = −0.0097 < 0. Hence, r k is not always greater than 0.
where n ∈ N. Suppose that From theorem 3.5, the number of k satisfying r k < 0 is finite, and Equation (15) z k has at most three solutions, we assert that such a τ 0 is attainable. When τ = τ 0 , Equation (14) has a pair of conjugate purely imaginary roots. Note that all roots of Equation (14) have negative real parts if τ = 0, and all roots of Equation (14) have negative real parts if τ < τ 0 , which follows from the continuous dependence of the solution on τ . This completes the proof.

Properties of Hopf bifurcations
By Theorem 3.6, we know that τ 0 exists, which corresponds to a certain k 0 ∈ N ∪ {0}. For the simplicity of notations, we still denote it by k in the following analysis.
In the following, we shall determine when the equilibrium E * becomes unstable and Hopf bifurcation occurs. We seek conditions which guarantee that characteristic equation (3.6) has a pair of conjugate purely imaginary roots.

Theorem 4.1:
Suppose that 0 > 1 , r k < 0 and p 2 k < 3q k , then there exists τ 0 > 0 such that Hopf bifurcation occurs at τ = τ 0 . That is, a periodic solution to (2) appears bifurcating from equilibrium E * as τ passes through the critical value τ 0 . Proof: Consider the derivative of the root of (14) with respect to τ at τ = τ 0 . Assume that λ(τ ) = ξ(τ ) + iω(τ ) is the eigenvalue of Equation (14) at τ near τ 0 . Directly calculating the differential of Equation (14) with respect to τ , we obtain This gives Note that ξ(τ 0 ) = 0 and ω(τ 0 ) = ω * , then 2 ]ω 2 * > 0. Since z 3 * + p k z 2 * + q k z * + r k = 0 and r k < 0, we have 3z 2 We have already shown the existence of Hopf bifurcation, which implies there exists a family of periodic solutions bifurcating from the steady state solution E * of the system (2) when the parameter τ crosses through the critical values τ 0 . Next, we will give a formula by using the center manifold theorem and the normal form theory [8,22] of partial functional differential equations to determine the direction of Hopf bifurcation and stability of periodic solutions.
Then (16) can be rewritten asU Define cos kx/l cos kx/l , cos kx/l = lπ 0 cos 2 kx/ldx 1 2 . Denote From the Riesz representation theorem, there exists a matrix η k (θ , μ), whose components are bounded variation functions, such that Let C * = C([0, τ ], X) and a bilinear form (·, ·) on C * × C is defined as follows: It follows from the definition of b k and l 0 b k b j = 0 for k = j that the above bilinear form is equivalent to where (·, ·) c is bilinear form defined on C * × C and satisfies Now we can define the adjoint the operator A * as below, Let q(θ )b k = (1, α 1 , β 1 ) e iω k θ b k and q * (s)b k = M(1, α 2 , β 2 ) e iω k s b k be the eigenvectors of A and A * corresponding to iω k and −iω k respectively. By some tedious calculations, it has w(t) = zα 1 +zᾱ 1 + W (2) Direct substitution and comparing the coefficients with (22) show that 20 (−τ ) 11 (−τ ))]Mb k .
In addition, we obtain that From the define of A μ and (24), it has Therefore, Furthermore, Consequently, g 21 can be expressed in terms of the parameters. Then the following values can be calculated explicitly, (i) If μ 2 > 0(< 0), then the Hopf bifurcation is supercritical (subcritical) at τ = τ 0 .

Numerical simulations
In this section, we will give some numerical simulations to analyse the dynamic properties of the malaria model.   . Furthermore, for the purpose of illustration, we choose τ = 1 and initial values (u 1 , u 2 , u 3 ) = (3.975 × 10 6 , 7000, 1.9 × 10 7 ). Then 0 = δ d 1 d 2 = 0.0025 < 1. By Theorem 3.1, the infection-free equilibrium E 0 = (5 × 10 6 , 0, 0) is globally asymptotically stable for any τ . Meanwhile, this result shows that iRBCs can be eliminated by immune response in host infected with malaria, such that malaria infection cannot be established within a host if 0 ≤ 1 (see Figure 2).
If we choose δ = 3 × 10 −7 , β = 1.5 × 10 −8 and the remaining parameter values remain the same as in Figure 1, then 0 = 1.5 > 1, 1 ≈ 121. The infection equilibrium (without specific immune response) E 1 is globally asymptotically stable. It means that the quantity of immune factors in the host is degraded gradually. That is, immune factors do not remove the iRBCs in the host completely, and malaria infection can be established in the host if 1 < 0 < 1 (see Figure 3).
In the case of delay τ > 0, taking τ = 0, δ = 3 × 10 −6 , β = 7 × 10 −6 , initial values (u 1 , u 2 , u 3 ) = (4.975 × 10 6 , 2.4 × 10 5 , 2.55 × 10 9 ) and the remaining parameters are the same as in Figure 2, then 0 = 15 and 1 ≈ 1.26. The infection equilibrium (with specific immune response) E * is globally asymptotically stable. It shows that immune factors play roles in the host, but the number of immune factors and iRBCs will eventually tend to a dynamic balance. Immune factors cannot completely eliminate the iRBCs in the host, and   malaria infection can be established if τ = 0 and 0 ≥ 1 (see Figure 4). Furthermore, we tried to change the diffusion parameters, where D = (D 1 , D 2 , D 3 ) = (10 −4 , 10 −2 , 1), but it seems to be no significant effect on the dynamics of the model (see Figure 5). Moreover, when τ < τ 0 ≈ 0.238, it follows from the fact that small delays are harmless, E * is locally asymptotically stable.
When τ increases and passes through τ 0 ≈ 0.238, we fix δ = 3 × 10 −6 , β = 7 × 10 −4 , τ = 1 > τ * ≈ 0.238 and the remaining parameter values are the same as in Figure  2, then 0 = 15 and 1 ≈ 3.6. E * becomes unstable and a periodic oscillation occurs near E * . Theorem 4.1 implies that system undergoes Hopf bifurcation and a periodic solution appears (see Figure 6). Furthermore, we change the diffusion parameters, where D = (D 1 , D 2 , D 3 ) = (10 −5 , 10 −1 , 1), but it did not play a decisive role (see Figure 7). It shows that in the previous period of relative stability, the number of iRBCs is very small. Meanwhile, the number of immune factors is small as well, the iRBCs in the host has not been cleared completely. We consider that the incubation period of malaria leads to the amount of iRBCs rising again after a period of time. Maybe that is why malaria might rekindle in the host.

Discussion
In this paper, we studied the dynamics of malaria infection model with diffusion and delay. For this biological mathematical model, we give two biologically meaningful threshold  values, the basic reproduction number 0 and the immune response threshold value 1 . The uninfected equilibrium E 0 is globally asymptotically stable for any delay and diffusion coefficient when 0 ≤ 1. It turns out that the malaria is cleared and immune response is not active. When 1 < 0 ≤ 1 , the infection equilibrium E 1 without immune response exists and it is globally asymptotically stable for any delay and diffusion. It signifies that the immune response would not be activated and malaria infection can be established in the host.
When 0 > 1 , a set of sufficient conditions, guaranteeing that malaria infection can be established in the host, are given as below: (i) If τ = 0, the infection equilibrium with immune response E * is globally asymptotically stable. (ii) If τ < τ 0 , the infection equilibrium with immune response E * is locally asymptotically stable. (iii) If τ > τ 0 , the dynamical behaviours of equilibrium E * will occur and locally asymptotically stable becomes unstable and Hopf bifurcation appears. It turns out that disease recurs in the host. The direction and stability of Hopf bifurcation is derived by applying the center manifold method and the normal form theory.
Finally, numerical simulations are provided to demonstrate the theoretical results.

Disclosure statement
No potential conflict of interest was reported by the author(s).