Parameter estimates of the 2016-2017 Zika outbreak in Costa Rica: An Approximate Bayesian Computation (ABC) approach

1.   Introduction

Malaria is one of the most dangerous global health problems. It causes millions of infections every year, and more than 1.1 million people died from the disease, especially infants, young children and pregnant women (WHO ^[1]). The mainly affected group is children under the age of five, accounting for 82% of all malaria deaths. It is endemic in the tropical and subtropical regions of the world, and people in Africa are at the risk of such disease accounting for 27% of the world's malaria infections (WHO ^[1]). It is reported that about every 30 seconds, a child's life will be threatened by such disease. Although some children are survived, they still suffer miserable physical problems, such as hearing impairment and brain damage. Meanwhile, pregnant women are also susceptible to malaria infection, which can cause maternal mortality, low birth weight in infants as well as maternal anaemia. According to the statistics, it caused an estimated 243 million malaria cases led to an estimated 863,000 deaths in 2008 (WHO ^[1]).

Malaria infection is mainly induced by the following four malaria parasites: P. falciparum, P. malariae, P. ovale and P. vivax. Among them, P. falciparum is responsible for almost all of the deaths attributed to malaria ^[2]. And anopheles is a vector of malaria. Once the infected anopheles bites the 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 parasites develop and mature in the biting mosquitos until they can infect others. When malaria parasites evolve in the host, they can stimulate the activity of immune cells which produce an immune response.

In the past few years, lots of efforts have been carried out to control and prevent such disease ^[3]. However, people living in high epidemic areas are still hardly to obtain effective malaria prevention, diagnosis and treatment. As the emergence of malaria drug resistant strains and insecticide resistant mosquitoes, how to prevent and control malaria infection is still a problem. The main detriment of plasmodium in the host occurs in the blood-stage, where the parasites will develop and reproduce. In order to clarify the mechanism, many mathematical models have been employed to describe the dynamics of malaria infection.

The first model was proposed by Hetzel and Anderson ^[4] which studied the performance of a mathematical model for the blood phase of malaria infection with a single strain. Analyzing the cells population dynamics in the absence of a host immune response, the authors demonstrated a relationship between host and parasite parameters, which defined a criterion for the successful invasion and persistence of the parasite. Hetzel and Anderson ^[4] indicated some important parameters, such as the production rates of merozoite and erythrocyte, the mortality rates of merozoite and erythrocyte as well as the invasion rate of merozoite.

In 2011, Li et al. ^[5] studied the dynamics for malaria parasites infection in the host blood-stage. Considered that the clearance of host immune cells was restricted by concentration, the following model was established:

where $H, I, M$ and $E$ are population of red blood cells, infected red blood cells, malaria parasites and immunity effectors, respectively. $d_{1}, \delta, \mu$ and $d_{2}$ represent the decay rate of $H(t)$, $I(t)$, $M(t)$ and $E(t)$ respectively. $\lambda$, $\alpha$ and $r$ are production rate of $H(t)$, infection of $H(t)$ by $M(t)$, product rate of $M(t)$ respectively. $p_{1}, p_{2}, k_{1}, k_{2}, \beta$ and $\gamma$ represent the removal rate of $I(t)$ and $M(t)$ by $E(t)$, proliferation rate of $E(t)$ by $I(t)$ and $M(T)$, $\frac{1}{\beta}$ means half saturation constant for $I(t)$ and $\frac{1}{\gamma}$ is half saturation constant for $M(t)$ respectively. In ^[5], the authors showed that there existed a threshold value $\Re_{0}$ and the malaria-free equilibrium was globally asymptotically stable if $\Re_{0} < 1$, if $\; \Re_{0} > 1$, there existed two kinds of infection equilibria: malaria infection equilibrium (without specific immune response) and positive equilibrium (with specific immune response). Conditions on the existence and stability of both infection equilibria are given. Moreover, it indicated that the model could undergo Hopf bifurcation at the positive equilibrium and exhibit periodic oscillations.

A malaria model between human and mosquitoes was established by Xiao and Zou ^[6] which presented an important conclusion that two malaria parasites may co-exist in the same region but not in the same host.

In 2016, Chang ^[7] made a great progress in the study of the innate immune mechanism of plasmodium falciparum escaping the immune clearance in the host. Innate immune response is the first line for host to defense pathogen invasion. After being stimulated by the red blood cells, the host neutrophils release chromatin and lysosomal granules in the cytoplasm to form a reticular structure (NETs) in a way of active death (Netosis). NETs can capture and kill infected red blood cells. At the same time, plasmodium degrades NETs by secreting DNA enzymes (TatD-like DNase) to escape host immune elimination. Chang ^[7] indicated that TatD-like DNase was an important factor in the survival of malaria parasites as well as a potential candidate for malaria vaccine.

In terms of a great deal of literatures in malaria infection ^{[8,9,10,11,12]}, we focus on some more important factors and hope to establish new models that can be used to study malaria infection in a host from different perspectives. Liu ^[13] considered three dynamical variables of populations: RBCs $T(t)$, iRBCs $I(t)$ and the immunity effectors (NETs) $M(t)$. Here the model reads as

The variables and parameters are given in the table below.

The main results of model (1.1) have been shown in ^[13]. And here, we just simply introduce the results as below.

(ⅰ) If $\Re_{0} < 1$, then the infection-free equilibrium is globally asymptotically stable.

(ⅱ) If $1 < \Re_{0} < 1+\frac{d_3\mu}{d_1\beta}$, then the infection equilibrium (without specific immune response) is globally asymptotically stable.

(ⅲ) If $\Re_{0} > 1+\frac{d_3\mu}{d_1\beta}$, then the infection equilibrium (with specific immune response) is globally asymptotically stable.

However, in real life, malaria infection does not tend to a fixed state over time, it turns out that some diseases will recur over a period of time or even change periodically. That means, it is not a good way to study the dynamics of malaria in host by using an ordinary differential equation model. Hence, in order to describe the time lag of the generation of immune factors, we consider adding a delay term to model (1.1). The delay differential equation model is as follows.

Compared with model (1.1), this model adds a time delay $\tau$ in the production item of immune factors $\beta I(t)M(t)$. In the following analysis, we calculated the characteristic root of the characteristic transcendental equation and found that the stability of positive equilibrium is changed from stable to unstable, and thus the Hopf bifurcation occurs.

This paper is organized as follows. In section 2, we propose a delay differential equation model for the within-host dynamics of malaria infection based on the basic reproduction number, existence of equilibria, stability of infection-free equilibrium $E_{0}$ and infection equilibria $E_{1}, E^*$ as well as persistence of positive equilibrium $E^{*}$. Section 3 is devoted to exhibit periodic oscillations due to the Hopf bifurcation at the positive equilibrium $E^{*}$. Specifically, by choosing immune response delay as bifurcation parameter, we can demonstrate that a limit cycle occurs via Hopf bifurcation when the time delay $\tau$ passes through the critical value. The direction and stability of Hopf bifurcation are derived by applying the center manifold method and the normal form theory. In section 4, some numerical simulations are presented to interpret our main results biologically. A short discussion on research results is also given in this section.

2.   Global dynamical properties of the model

Consider the delay differential equation model

The initial conditions of (2.1) are given as

where $(\varphi_{1}, \varphi_{2}, \varphi_{3}) \in \mathbf{C([-\tau, 0], \mathrm{R}_+^3)}$, the Banach space of continuous functions map the interval $[-\tau, 0]$ into $\mathrm{R}_+^3$. By fundamental theory of functional differential equation (Kuang ^[14]), there exists a unique solution $(T(t), I(t), M(t))$ of (2.1) with initial conditions (2.2).

Upon the positivity and boundedness of solutions for model (2.1), we claim the following result.

Theorem 2.1. Let $(T(t), I(t), M(t))$ be the solution of model (2.1) with initial conditions (2.2), then $T(t), I(t), M(t)$ are positive and ultimately bounded for all $t\geq0$.

Proof. The positivity of the model (2.1) follows directly from Theorem 3.4 ^[15]. As for the boundedness of solutions to the model (2.1) with initial conditions (2.2), we denote

By the positivity of solutions, we have

where $m = \min{\{d_{1}, d_{2}, d_{3}\}}$. Hence, it comes

This implies that $T(t), I(t), M(t)$ are ultimately bounded.

Next, we define threshold values $\Re_{0}$ and $\Re_{1}$. For model (2.1), the basic reproduction number of iRBCs, which describes 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

We will see that $\Re_{0}$ and $\Re_{1}$ play key roles in determining the existence and stability of equilibria of the model (2.1). Actually, let

The following results can be verified by direct calculations.

(ⅰ) If $\Re_{0}\leq1$, then (2.1) always has a unique infection-free equilibrium $E_0 = (T_0, 0, 0)$, where $T_0 = \frac{\lambda}{d_1}$.

(ⅱ) If $1 < \Re_{0}\leq\Re_{1}$, then (2.1) has two equilibria, $E_0$ and infection equilibrium (without specific immune response) $E_1 = (T_{1}, I_{1}, 0)$, where $T_1 = \frac{d_{2}}{\mu}$ and $I_1 = \frac{\lambda \mu-d_{1}d_{2}}{d_{2}\mu} = \frac{d_1(\Re_{0}-1)}{\mu}$.

(ⅲ) If $\Re_{0} > \Re_{1}$, then (2.1) has three equilibria, $E_0, E_1$ and infection equilibrium (with specific immune response) $E^* = (T^*, I^*, M^*)$ respectively, where $T^* = \frac{\lambda\beta}{d_{1}\beta+d_{3}\mu}, I^* = \frac{d_{3}}{\beta}$ and $M^* = \frac{\lambda\mu\beta-d_{1}d_{2}\beta-d_{2}d_{3}\mu}{d_{1}\alpha\beta+d_{3}\alpha\mu} = \frac{d_1d_2\beta(\Re_0-\Re_1)}{d_1\alpha\beta+d_3\alpha\mu}$.

Now, we investigate the stability of the equilibria of model (2.1). Consider the linearization of the model (2.1) at $E = (T(t), I(t), M(t))$,

where $x(t) = (T(t), I(t), M(t))^{T}$,

The characteristic equation of (2.5) is formulated as (Hale ^[16])

Theorem 2.2. If $\Re_{0} < 1$, then the infection-free equilibrium $E_{0}$ of model (2.1) is locally asymptotically stable. Moreover, $E_{0}$ is unstable if $\Re_{0} > 1$.

Proof. By substituting $E_0$ in (2.5), (2.6), (2.7) and (2.8), we obtain the following characteristic equation

Clearly, $\Lambda_1 = -d_{1}$, $\Lambda_2 = -d_{3}$ and

which dominates the stability of $E_0$. If $\Re_0 < 1$, then $\Lambda_3 < 0$. Therefore, when $\Re_0 < 1$, all eigenvalues of model (2.1) have negative real parts and hence, equilibrium $E_0$ is locally asymptotically stable. In addition, when $\Re_0 > 1$, we have $\Lambda_3 > 0$, which means $E_0$ is unstable.

Theorem 2.3. If $\Re_{0}\leq1$, then the infection-free equilibrium $E_{0}$ of model (2.1) is globally asymptotically stable.

Proof. Define Lyapunov functional $V_0$ as follows

Differentiating $V_0$ with time $t$ along the model (2.1), then

Obviously, if $\Re_0\leq1$, then $\frac{dV_{0}}{dt}\leq0$ for any $T(t), I(t)$ and $M(t)$. In addition, $\frac{dV_{0}}{dt} = 0$ if and only if $T(t) = T_0, I(t) = 0$ and $M(t) = 0$. Let $\Gamma$ be the largest invariant set of $\{(T(t), I(t), M(t) \in \mathrm{R}_+^3\mid \frac{dV_{0}}{dt} = 0$}, we easily get $\Gamma = \{E_0\}$. According to LaSalles invariance principle (Kuang ^[14]), it follows that equilibrium $E_0$ is globally asymptotically stable when $\Re_0\leq1$. This completes the proof.

Theorem 2.4. If $1 < \Re_{0} < \Re_{1}$, then the infection equilibrium (without specific immune response) $E_{1}$ of model (2.1) is locally asymptotically stable.

Proof. By substituting $E_1$ into (2.5), (2.6), (2.7) and (2.8), we have the following characteristic equation

It is easy to check that $\Lambda^2 +d_{1}\Re_{0}\Lambda+d_{1}d_{2}(\Re_{0}-1) = 0$ has two eigenvalues $\Lambda_1$ and $\Lambda_2$ with negative real parts if $1 < \Re_{0} < \Re_{1}$. We assume that $\Lambda = a+ib$ satisfies

Substituting $\Lambda = a+ib$ into (2.10) and separating the real and imaginary parts, we obtain

Assume that $a\geq0$, when $1 < \Re_{0} < \Re_{1}$, then $a < d_{3}[e^{-a\tau}\cos(b\tau)-1]$ from the right hand side of the first equation of (2.11). This contradicts the assumption. So we get $a < 0$, i.e. $Re\Lambda_{3} < 0$. Thus, all eigenvalues of equation (2.9) have negative real parts and equilibrium $E_1$ is locally asymptotically stable.

Theorem 2.5. If $1 < \Re_{0}\leq\Re_{1}$, then the infection equilibrium (without specific immune response) $E_{1}$ of model (2.1) is globally asymptotically stable.

Proof. Define Lyapunov functional $V_1$ as follows

Differentiating $V_1$ with time $t$ along model (2.1), then

Obviously, if $1 < \Re_{0}\leq\Re_{1}$, then $\frac{dV_{1}}{dt}\leq0$ for any $T(t), I(t)$ and $M(t)$. In addition, $\frac{dV_{1}}{dt} = 0$ if and only if $T(t) = T_1, I(t) = I_1$ and $M(t) = 0$. By LaSalles invariance principle (Kuang ^[14]), equilibrium $E_1$ is globally asymptotically stable when $1 < \Re_{0}\leq\Re_{1}$. The proof is completed.

Next, we study the stability of the infected equilibrium (with specific immune response), which is denoted by positive equilibrium $E^{*}$. Recall that the positive equilibrium $E^{*}$ exists if and only if $\Re_{0} > \Re_{1}$. At equilibrium $E^{*}$, the characteristic equation for the corresponding linearized model of (2.1) is

where

When $\tau = 0$, equation (2.12) becomes

where

Clearly, $(A+A_1)(B+B_1) > C+C_1$. By Hurwitz criteria, all eigenvalues of equation (2.13) have negative real parts, which leads to the local stability of equilibrium $E^*$ when $\tau = 0$. In fact, in this case, $E^*$ is also globally stable.

Theorem 2.6. If $\Re_{0} > \Re_{1}$ and $\tau = 0$, then the infection equilibrium (with specific immune response) $E^{*}$ of model (2.1) is globally asymptotically stable.

Proof. Define Lyapunov functional $V^*$ that

Differentiating $V^*$ with respect to $t$ along model (2.1), we get

If $\Re_{0} > \Re_{1}$, then $\frac{dV^*}{dt}\leq0$ for any $T(t), I(t)$ and $M(t)$. In addition, $\frac{dV^*}{dt} = 0$ if and only if $T(t) = T^*, I(t) = I^*$ and $M(t) = M^*$. By using LaSalles invariance principle (Kuang ^[14]), equilibrium $E^*$ is globally asymptotically stable when $\Re_{0} > \Re_{1}$. This completes the proof.

Now we consider the stability of the equilibrium $E^*$ if $\Re_{0} > \Re_{1}$ and $\tau > 0$. By theorem $4.4$ of Hal Smith ^[15], for small enough delay, the characteristic roots of (2.12) are either very near the eigenvalues of (2.13) or having more negative real parts than any of the eigenvalues of (2.13). Hence, when the delay is small the equilibrium $E^*$ is locally asymptotically stable.

When $\Re_{0} > \Re_{1}$, for any $\tau > 0$, zero is not a root of (2.12). Note that any complex roots to the equations (2.12) appear in pairs, and all roots of (2.12) have negative real parts if $\tau = 0$. Therefore, any root of (2.12) has negative real part for sufficiently small $\tau$. Assume that there exists $\tau = \tau^*$, such that (2.12) has a pair of pure imaginary roots, denoted by $\Lambda = \pm\omega i, (\omega > 0)$. Substituting $\Lambda = \omega i$ into (2.12), we have,

Separate the real and imaginary parts,

It follows that

Let $z = \omega^2$, it follows that

where

By calculating $p_{1}, q_{1}$ and $r_{1}$, we obtain

Equation (2.15) has no positive roots if it satisfies $p_1 > 0, \; r_1 > 0, \; p_1q_1-r_1 > 0$ (Hurwitz criterion). Hence, $\tau^*$ doesn't exist. It follows that all solutions of equation (2.15) have negative real part and $E^*$ is locally asymptotically stable for any $\tau > 0$. Therefore, we have the following theorem.

Theorem 2.7. If $\Re_{0} > \Re_{1}$ and $p_1 > 0, \; r_1 > 0, \; p_1q_1-r_1 > 0$, then the infection equilibrium $E^*$ is locally asymptotically stable for any $\tau > 0$.

Theorem 2.8. If $\Re_{0} > \Re_{1}$ and $r_{1} < 0$, there exists $\tau_* > 0$, such that equation (2.15) has a pair of conjugate purely imaginary roots $\pm\omega_{*}i$ when $\tau < \tau_{*}$. Moreover, the infection equilibrium $E^*$ is locally asymptotically stable if $\tau < \tau_{*}$.

Proof. We know that if $r_{1} < 0$, then equation (2.15) has at least one positive root and not more than three positive roots. Assume that equation (2.15) has three positive roots, denoted by $z_{k}, k = 1, 2, 3$, it follows that

where $k = 1, 2, 3, \; n\in Z^{+}$. Suppose that

When $\tau = \tau_{*}$, equation (2.15) has a pair of conjugate purely imaginary roots. Note that all roots of equation (2.15) have negative real part if $\tau = 0$, and all roots of equation (2.15) have negative real part if $\tau < \tau^*$, which follows from the continuous dependence of the solution on $\tau$.

In what follows, we study the persistence of equilibrium $E^*$, following the analysis of Wang et al. ^[17].

Lemma 2.1. Let $X$ be a locally compact metric space and it is the union of two disjoint subsets $X_1$ and $X_2$, with $X_2$ be compact in $X$, $X_1$ be open and forward invariant under the continuous semiflow $\Phi$ on $X$. Assume that

has an acyclic isolated covering $M = \cup_{k = 1}^{m} M_k$. If each part $M_k$ of $M$ is a weak repeller for $X_1$, then $X_2$ is a uniform strong repeller for $X_1$.

Theorem 2.9. If $\Re_{0} > \Re_{1}$, then $M$ is uniformly persistent, that is, there exists a positive constant $\sigma$, such that the positive solutions of model (2.1) satisfy $\liminf\limits_{t\rightarrow \infty} M(t) > \sigma$.

Proof. By Theorem 2.1, model (2.1) is point dissipative. Let

To prove Theorem 2.9, we need to show that $X_2$ is uniformly strong repeller for $X_1$. It suffice to verify the conditions of Lemma 2.1. Obviously, $X$ is locally compact, $X_2$ is compact and $X_1$ is forward invariant. In addition, there are two equilibria $E_{0} = (\frac{\lambda}{d_{1}}, 0, 0)$ and $E_{1} = (\frac{d_{2}}{\mu}, \frac{\lambda\mu-d_{1}d_{2}}{d_{2}\mu}, 0)$ in $X_2$. Consider

By Theorem 2.5, $E_{11} = (T_{1}, I_{1}) = (\frac{d_{2}}{\mu}, \frac{\lambda\mu-d_{1}d_{2}}{d_{2}\mu})$ is globally asymptotically stable for solutions $(T(t), I(t))$ to (2.16) with positive initial value if $1 < \Re_{0}\leq\Re_{1}$. In addition, let $Y = \{(T, 0)|T\geq0\}$, then $Y$ is positively invariant for (2.16) and the solutions to (2.16) in $Y$ approach $E_{00} = (T_{0}, 0) = (\frac{\lambda}{d_{1}}, 0)$ as $t\rightarrow \infty$, it becomes zero or unbounded as $t\rightarrow -\infty$. Therefore, it follows that

Thus, we know that $M$ is acyclic in $X_2$. From theorem 2.4, 2.5 and the Jacobian matrix at $E_{00}$, it is easy to check that $E_{00}$ and $E_{11}$ are hyperbolic if $\Re_{0} > \Re_{1}$, which implies that $E_{00}$ and $E_{11}$ are isolated in $X_2$. From what has been discussed above, $\Omega_2$ has an acyclic isolated covering $M = E_{00}\cup E_{11}$.

Next, we need to show that each part $M_k (k = 1, 2)$ of $M$ is a weak repeller for $X_1$. It suffice to show that $W^s(E_{11})\bigcap X_1 = \emptyset$ and $W^s(E_{00})\bigcap X_1 = \emptyset$, where $W^s(E_{11})$ denotes the stable manifold of $E_{11}$, and $W^s(E_{00})$ denotes the stable manifold of $E_{00}$. Suppose $W^s(E_{11})\bigcap X_1\neq\emptyset$, then there is a solution of model (2.1) such that ${\lim_{t\rightarrow\infty}}(T(t), I(t), M(t))\rightarrow E_{11}$. Therefore, from the third equation of model (2.1) that for any $\varepsilon > 0$, there exists $t_0 > 0$, then we obtain

Note that

then for $t\ge t_0$ and $M'(t)\ge g(\varepsilon)M$. Since $\Re_{0} > \Re_{1}$, we can choose $\varepsilon$ sufficiently small such that $g(\varepsilon) > 0$. Thus, (2.17) shows that $M(t)\rightarrow\infty$ as $t\rightarrow\infty$. This contradicts the ultimate boundedness of Theorem 2.1. Hence, $W^{s}(E_{11})\bigcap X_{1} = \emptyset$. Similarly, we can prove $W^s(E_{00})\bigcap X_1 = \emptyset$.

3.   Hopf bifurcation

In the following, we shall determine when the equilibrium $E^*$ becomes unstable and Hopf bifurcation occurs. We will seek conditions which guarantee characteristic equation (2.15) to have a root with negative real part as well as a pair of conjugate purely imaginary roots.

Theorem 3.1. Suppose $\Re_{0} > \Re_{1}$ and $r_{1} < 0$, then there exists $\tau^* > 0$. If $\tau > \tau_{*}$, there is a Hopf bifurcation of model (2.1) from equilibrium $E^*$ as $\tau$ passes through the critical value $\tau^*$.

Proof. Consider the derivative of eigenvalue on real axis across section of positive equilibrium point $E^*$ on $\tau = \tau_{*}$ (^[18]). Assume that $\Lambda(\tau) = \xi(\tau)+i\omega(\tau)$ is the eigenvalue of equation (2.15) at $\tau$ near $\tau_{*}$ and calculate the derivative of equation (2.15) on $\tau$, we obtain (^[19])

This gives

For $\xi(\tau_{*}) = 0$ and $\omega(\tau_{*}) = \omega_{*}$, then

where $Z = [b_{1}^2\omega_{*}^2+(C_{1}-A_{1}\omega_{*}^2)^2]\omega_{*}^2 > 0$. Since $z_{*}^3+p_{1}z_{*}^2+q_{1}z_{*}+r_{1} = 0$ and $r_{1} < 0$, then $3z_{*}^2+2p_{1}z_{*}+q_{1} > 0$ if either $p_{1}\geq 0$ or $q_{1}\leq0$. So we have

if either $p_{1}\geq0$ or $q_{1}\leq0$.

Remark 3.1. By theorem 3.1, we can extend that if there exists $p_{1}, q_{1}$ and $r_{1}$, which makes $z_{*}^3+p_{1}z_{*}^2+q_{1}z_{*}+r_{1} = 0$ and $\frac{z_{*}}{Z}(3z_{*}^2+2p_{1}z_{*}+q_{1}) > 0$ hold. Then periodic solutions are bifurcated near the positive equilibrium $E^*$.

We have already shown the existence of Hopf bifurcation. Next, we will give a formula by using the center manifold theorem and the normal form theory to determine the direction of Hopf bifurcation and stability of periodic solutions.

Let $v(t) = T(\tau t)-T^*, x(t) = I(\tau t)-I^*$ and $y(t) = M(\tau t)-M^*$. Then (2.1) can be rewritten as

where

$A = \left(\begin{array}{ccc} -d_1-\mu I^* & -\mu T^* & 0 \\ \mu I^* & 0 & -\alpha I^* \\ 0 & 0 & -d_3 \\ \end{array} \right),$

$B = \left(\begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & \beta M^* & \beta I^* \\ \end{array} \right),$

$F(v_t, x_t, y_t, \tau) = \left(\begin{array}{c} -\tau \mu v(t)x(t) \\ \tau \mu v(t)x(t)-\tau \alpha x(t)y(t)\\ \tau \beta x(t-1)y(t-1)-\tau d_3 M^*\\ \end{array} \right).$

Let $\hat{\tau}$ be the critical value of $\tau$ where model (2.1) undergoes a Hopf bifurcation at $E^*$. Assume $\tau = \hat{\tau}+h$, then $h = 0$ is the Hopf bifurcation value of (2.1).

Choose the phase space $C = C([-1, 0], R^3)$. Define $L(h):C \rightarrow R^3$ by

From the Riesz representation theorem, there exists a matrix whose components are bounded variation functions $\eta (\theta, h):[-1, 0]\rightarrow R^3, \theta \in [-1, 0]$, such that

We select $\eta (\theta, h) = (\hat{\tau} +h)A\delta (\theta)-(\hat{\tau} +h)B\delta (\theta+1)$, where

For $\phi \in C^1 ([-1, 0], R^3)$, we give

such that

Then (3.1) can be rewritten as

For $\varphi \in C^1([0, 1], (C^3)^*)$, define

And for $\phi \in C([-1, 0], C^3)$ with $\varphi \in C([0, 1], (C^3)^*)$, define

then $A^*$ and $A(0)$ are adjoint operators, such that $\langle \varphi, A\phi\rangle = \langle A^*\varphi, \phi\rangle$. Let $q(\theta)$ and $q^*(s)$ be the eigenvectors of $A$ and $A^*$ corresponding to $iw\hat{\tau}$ and $-iw\hat{\tau}$ respectively. Via direct calculation, we obtain following result.

Lemma 3.1. $q(\theta) = (1, \alpha_1, \beta_1)^T e^{i\hat{w}\hat{\tau} \theta}$ is the eigenvector of operator A on $i\hat{w}\hat{\tau}$, $q^*(s) = \overline{D} (1, \alpha_2, \beta_2)^T e^{i\hat{w}\hat{\tau} s}$ is the eigenvector of operator $A^*$ on $-i\hat{w}\hat{\tau}$, and $\langle q^*, q \rangle = 1, \langle q^*, \bar{q} \rangle = 0$, where $\alpha_1 = -\frac{i\hat{w}+d_1+\mu I^*}{\mu T^*}$, $\beta_1 = -\frac{i\hat{w}\alpha_1-\mu I^*}{\alpha I^*}$, $\alpha_2 = \frac{-i\hat{w}+d_1+\mu I^*}{\mu I^*}$, $\beta_2 = \frac{-i\hat{w}\alpha_2+\mu T^*}{e^{i\hat{w}\hat{\tau}}\beta M^*}.$

Proof. Let $q(\theta) = (1, \alpha_1, \beta_1)^T e^{i\hat{w}\hat{\tau} \theta}$ be an eigenvector of operator A on $i\hat{w}\hat{\tau}$, it follows from $Aq(0) = i\hat{w}\hat{\tau} q(0)$, then

The calculation shows that $\alpha_1 = -\frac{i\hat{w}+d_1+\mu I^*}{\mu T^*}$, $\beta_1 = -\frac{i\hat{w}\alpha_1-\mu I^*}{\alpha I^*}$. Similarly, we can obtain $\alpha_2 = \frac{-i\hat{w}+d_1+\mu I^*}{\mu I^*}$, $\beta_2 = \frac{-i\hat{w}\alpha_2+\mu T^*}{e^{i\hat{w}\hat{\tau}}\beta M^*}.$

To ensure $\langle q^*, q \rangle = 1$, we only need

The details have been given in Hassard et al. ^[20] corroborated $\langle q^*, \bar{q} \rangle = 0$. Here we are not going to repeat it. The proof is completed.

Following the algorithm in Hassard et al. ^[20], we first compute the coordinates to describe the center manifold $C_0$ at $h = 0$. Let $\mu _t = (v_t, x_t, y_t)^T$ be the solution of (3.1) when $\tau = \hat{\tau}$, i.e. when $h = 0$. Define

On the center manifold $C_0$, it shows

where

$z$ and $\bar{z}$ represent the local coordinates for the center manifold $C_0$ in the direction of $q^*$ and $\bar{q}^*$. Note that $W$ is real if $u_t$ is real, and here, we only consider real solutions. For the solution $u_t \in C_0$ of (3.1), since $h = 0$,

We rewrite the above equation as

where

From (3.2) and (3.3), we have

where

Expanding the above series and comparing the coefficients, it yields

Note that

and

Direct substitution and comparing the coefficients with (3.4) shows that

We still need to compute $W_{20}(\theta)$ and $W_{11}(\theta)$. For $\theta \in [-1, 0]$, it comes

Comparing the coefficients with (3.5), we found that

and

It follows from (3.6) that

Then solving the above equation, we get

Similarly,

where $E_1$ and $E_2$ are both two-dimensional vectors and can be determined by setting $\theta = 0$ in $H$. In fact, since

we have

and

It follows from (3.6) and the definition of $A$ that

and

Substituting (3.7), (3.8) into the above two equations respectively yields

where $I$ is the $3\times3$ identity matrix. Consequently, $g_{21}$ can be expressed in terms of the parameters and delay $\hat{\tau}$. Then the following values can be computed,

Utilizing the results of Hassard et al. ^[20], we make a conclusion as follows.

Theorem 3.2.

(ⅰ) If $\mu_2 > 0(<0)$, then the Hopf bifurcation is supercritical (subcritical).

(ⅱ) If $\beta_2 < 0(>0)$, then the bifurcating periodic solutions are stable (unstable).

4.   Simulations, biological explanations and discussions

In order to interpret the conclusions from a quantitative perspective, the dynamics of malaria infection in the host by numerical simulations will be analyzed in the following. In this section, we use matlab to find the numerical solutions of the model (2.1) and analyze the effect of basic reproduction number $\Re_0$ and the immune response reproduction number $\Re_1$ on malaria infection.

With parameters values given in Table 1 and Table 2, we get $\Re_0 = 0.0025 < 1$. Hence, Theorem 2.3 indicates that the infection-free equilibrium $E_0$ is globally asymptotically stable. Meanwhile, this result shows that iRBCs can be eliminated by immune response in host infected with malaria, such that malaria infection can not be established within a host if $\Re_{0}\leq1$ (see Figure 1).

Taking $u = 3\times10^{-7}, \beta = 1.5\times10^{-8}$ and the other parameters values as in Table 1 and Table 2, we get $\Re_0 = 1.5$ and $\Re_1\approx121$. The infection equilibrium (without specific immune response) $E_1$ is globally asymptotically stable. It means that the generation of immune factors in the host are degraded immediately, which ultimately leads to the failure of the immune factors in the host, that is, immune factors do not remove the iRBCs in the host completely, such that malaria infection can be established in the host if $\Re_0 > 1$ (see Figure 2).

Taking $\tau = 0, u = 3\times10^{-7}, \beta = 7\times10^{-6}$ and the other parameters values as in Table 1 and Table 2, we get $\Re_0 = 1.5$ and $\Re_1\approx1.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 can not completely eliminate the iRBCs in the host, such that malaria infection can be established in the host if $\tau = 0$ and $\Re_0\geq\Re_1$ (see Figure 3).

When all parameters values are as in Figure 3 with $\tau = 0.1$, we can find the solution curve is the same as Figure 3. It means that there is no great change in the stability of the equilibrium $E^*$ when the time delay $\tau$ is small. Immune factors can not completely eliminate the iRBCs in the host, so malaria infection can be established in the host if $\tau = 0.1$ and $\Re_0 > \Re_1$ (see Figure 4).

When $\tau$ increases and passes through $\tau^*\approx0.238$. For example, taking $\tau = 2$, we get $\Re_0 = 1.5, \Re_1\approx1.26$ and $E^*$ becomes unstable. Theorem 3.1 implies that model (2.1) undergoes Hopf bifurcation and a periodic solution appears. It means 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 for malaria leads to the number of iRBCs rising again after a period of time. Maybe that is why malaria might rekindle in the host (see Figure 5).

In this paper, we have discussed a malaria infection model with time delay.

For the model (2.1), immune response delay $\tau$ plays a important role in the stability of the equilibrium $E^*$. We derive the basic reproduction number $\Re_0$ for the malaria infection and establish that the dynamics are completely determined by the basic reproduction number $\Re_0$. The results show that when $\Re_0\leq1$, the infection-free equilibrium $E_0$ is globally asymptotically stable for any delay $\tau$. It signifies that the malaria is cleared and immune response is not active. When $1 < \Re_0\leq\Re_1$, the infection equilibrium $E_1$ without immune response exists and it is globally asymptotically stable for any delay $\tau$, which means that the immune response would not be activated and malaria infection can be established in the host. When $\Re_0 > \Re_1$, some results are given as below:

(ⅰ) If $\tau = 0$, the infection equilibrium with immune response $E^*$ is globally asymptotically stable.

(ⅱ) If $\tau < \tau^*$, the infection equilibrium with immune response $E^*$ is locally asymptotically stable.

(ⅲ) If model (2.1) satisfies theorem 3.1 or remark 3.1, the dynamical behaviors of equilibrium $E^*$ will occur and locally asymptotically stable becomes unstable and Hopf bifurcation appears. By choosing immune response delay as bifurcation parameter, we have demonstrated that a limit cycle occurs via Hopf bifurcation, when the delay passes through the critical value $\tau^*$. This explains the fact that the immune response delay plays negative role in controlling disease progression. The direction and stability of Hopf bifurcation is derived by applying the center manifold method and the normal form theory.

Moreover, we study the uniform persistence of infection equilibrium $E^*$.

Numerical simulations are also provided to demonstrate these theoretical results. Finally, we hope that our work will be helpful to the study of malaria infection.

Acknowledgements

The authors wish to thank two anonymous reviewers for their very helpful comments. QD, JL and ZG were supported by National Science Foundation of China (No. 11771104), Program for Chang Jiang Scholars and Innovative Research Team in University (IRT-16R16).

Conflict of interest

All authors declare no conflicts of interest in this paper.