A Malaria Model with Two Delays

Abstract

A transmission model of malaria with two delays is formulated. We calculate the basic reproduction number 𝑅₀ for the model. It is shown that the basic reproduction number is a decreasing function of two time delays. The existence of the equilibria is studied. Our results suggest that the model undergoes a backward bifurcation, which implies that bringing the basic reproduction number below 1 is not enough to eradicate malaria.

1. Introduction

Malaria is a mosquito-borne disease which is due to the four species of the genus Plasmodium, namely, Plasmodium falciparum, Plasmodium vivax, Plasmodium malariae, and Plasmodium ovale. These parasites are transmitted to the human host through a bite by an infected female anopheles mosquito.

Sporozoites are injected into a human host, which are carried through the blood to the liver within 30 min [1]. They invade hepatocytes and undergo a process of asexual replication (exoerythrocytic schizogony) to give rise to 10–40 thousand merozoites per sporozoite. Up to this point, the infection is nonpathogenic and clinically silent. After about 7–9 days, the liver schizonts rupture to release the merozoites into the blood and clinical symptoms, such as fever, pain, chills, and sweats, may develop. Each merozoite invades an erythrocyte and divides to form an erythrocytic schizont containing about 16 daughter merozoites [2]. These merozoites either reinfect fresh erythrocytes, giving rise to cyclical blood-stage infection with a periodicity of 48–72 h, depending on the Plasmodium species, or differentiate into sexual transmission stages called gametocytes. When a second mosquito bites the infected human, the gametocytes are ingested, giving rise to extracellular gametes. In the mosquito midgut, the gametes fuse to form a motile zygote (ookinete), which penetrates the mid-gut wall and forms an oocyst, within which meiosis takes place and haploid sporozoites develop. These sporozoites migrate to the salivary glands. The incubation period within the mosquito may last 8–22 days [3]. The variation in the length of time is due to the environmental temperature. For P. falciparum, the average time is 12 days [4, 5]. Malaria can also be transmitted through blood transfusion, organ transplantation and transplacental malaria (i.e., congenital malaria) can also be significant in populations which are partially immune to malaria.

The application of mathematics to study of infectious disease appear to have been initiated by Berniulli, 1760. Presently, a lot of mathematical models have been proposed to evaluate and compare control procedures and preventive strategies, and to investigate the relative effects of various sociological, biological, and environmental factors on the spread of diseases ([6–14], etc.). These models have played a very important role in the history and development of epidemiology.

The transmission process involves considerable time delay both in human and in mosquitoes due to the incubation periods of the several forms of the parasites [9]. Several works have been done to consider the effect of delay. Taking account of the incubation period in vector (mosquito), Takeuchi et al. [10] considered a differential-delay model for vector-borne diseases and get the global stability of the endemic equilibrium under appropriate conditions. While Wei et al. [15] proposed a similar model of vector-borne disease which has direct mode of transmission in addition to the vector-mediated transmission and proved that the introduction of a time delay in the host-to-vector transmission term can destabilize the system and periodic solutions can arise through Hopf bifurcation. In [9], Ruan et al. proposed a delayed Ross-Macdonald model with delays both in human and in mosquitoes. The effect of time delays on the basic reproduction number and the dynamics of the transmission was studied.

Both the models in [10, 15] adopted mass action formulation, where the contact rate depends on the size of the total host population while the model in [9] adopted standard incidence formulation, where the contact rate is assumed to be constant (see [16] for detailed derivation of these incidence functions). The choice of one formulation over the other really depends on the disease being modeled and, in some cases, the need for mathematical (analytical) tractability [17]. A crucial question is which of the two incidence functions is more suitable for modeling infectious disease in general, and malaria in particular? As noted in [16], data ([4, 18]) strongly suggests that standard incidence is better for modeling human diseases than mass action incidence. Consequently, standard incidence rate is adopted in this paper.

2. Derivation of the Model

In order to study the impact of incubation periods in both human and mosquitoes on the basic reproduction number and the transmission dynamics of malaria over long periods, we propose a model based on SIRS in human population and SI for the mosquito vector population. Since the mosquito dynamics operates on a much faster time scale than the human dynamics and the turnover of the mosquito population is very high and the total size of vector population is largely exceeding the human total size, the mosquito population can be considered to be at a equilibrium with respect to changes in the human population. Hence, the total number of mosquito population is assumed to be constant.

We formulate an SIRS-SI model for the spread of malaria in the human and mosquito population with the total population size at time 𝑡 is given by 𝑁_ℎ(𝑡) and 𝑁_𝑣(𝑡), correspondingly. For the human population, the three compartments represent individuals who are susceptible, infectious, and partially immune, with sizes at time 𝑡 denoted by 𝑆_ℎ(𝑡), 𝐼_ℎ(𝑡), and 𝑅_ℎ(𝑡), respectively. The incubation period in a human has duration 𝜏₁≥0. A susceptible individual infected by a mosquito will not become infectious until 𝜏₁ time units later and we assume that no one recovers during this incubation period. At the infectious stage a host may die from the disease, recover into the susceptible class, or may recover with acquired partial immunity. Infectious human hosts are those with infectious gametocytes in the blood stream. Partially immune hosts still have protective antibodies and other immune effectors at low levels. Λ_ℎ>0 is the human input (birth) rate. 𝜇_ℎ>0 and 𝛼_ℎ≥0 are the natural and disease-induced death rates, respectively. 𝛾_ℎ>0 is the rate at which human hosts acquire immunity. 𝜃≥0 is the per capita rate of recovery into the susceptible class from being infectious. 𝜌≥0 is the per capita rate of loss of immunity in human hosts. 𝛽_ℎ>0 is the proportion of bites on man that produce an infection. For the mosquito population, the two compartments represent susceptible and infectious mosquitoes, with sizes at time 𝑡 denoted by 𝑆_𝑣(𝑡) and 𝐼_𝑣(𝑡), respectively. The vector component of the model does not include immune class as mosquitoes never recover from infection that is, their infective period ends with their death due to their relatively short life cycle. Thus the immune class in the mosquito population is negligible and death occurs equally in all groups. Our model also excludes the immature mosquitoes since they do not participate in the infection cycle and are, thus, in the waiting period, which limits the vector population growth. We assume that the mosquito population is constant and equal to 𝑁_𝑣, with birth and death rate constants equal to 𝜇_𝑣. 𝛽_𝑣 is the probability that a mosquito becomes infectious when it bites an infectious human. The biting rate 𝑏 of mosquitoes is the average number of bites per mosquito per day. This rate depends on a number of factors, in particular, climatic ones, but for simplicity in this paper, we assume 𝑏 constant. The incubation interval in the mosquito has duration 𝜏₂≥0.

Considering the assumption made above, the interaction between human hosts and the mosquito vector population with standard incidence rate is described as shown below: 𝑑𝑆_ℎ𝑑𝑡=Λ_ℎ−𝑏𝛽_ℎ𝑁_ℎ𝑡−𝜏₁𝐼_𝑣𝑡−𝜏₁𝑆_ℎ𝑡−𝜏₁𝑒^{−𝜇ℎ𝜏1}−𝜇_ℎ𝑆_ℎ+𝜃𝐼_ℎ+𝜌𝑅_ℎ,𝑑𝐼_ℎ=𝑑𝑡𝑏𝛽_ℎ𝑁_ℎ𝑡−𝜏₁𝐼_𝑣𝑡−𝜏₁𝑆_ℎ𝑡−𝜏₁𝑒^{−𝜇ℎ𝜏1}−𝛾_ℎ+𝜇_ℎ+𝛼_ℎ+𝜃𝐼_ℎ,𝑑𝑅_ℎ𝑑𝑡=𝛾_ℎ𝐼_ℎ−𝜇_ℎ𝑅_ℎ−𝜌𝑅_ℎ,𝑑𝑆_𝑣𝑑𝑡=𝜇_𝑣𝑆_𝑣+𝜇_𝑣𝐼_𝑣−𝑏𝛽_𝑣𝑁_ℎ𝑡−𝜏₂𝐼_ℎ𝑡−𝜏₂𝑆_𝑣𝑡−𝜏₂𝑒^{−𝜇𝑣𝜏2}−𝜇_𝑣𝑆_𝑣,𝑑𝐼_𝑣=𝑑𝑡𝑏𝛽_𝑣𝑁_ℎ𝑡−𝜏₂𝐼_ℎ𝑡−𝜏₂𝑆_𝑣𝑡−𝜏₂𝑒^{−𝜇𝑣𝜏2}−𝜇_𝑣𝐼_𝑣,𝑑𝑁_ℎ𝑑𝑡=Λ_ℎ−𝜇_ℎ𝑁_ℎ−𝛼_ℎ𝐼_ℎ, (1) where 𝑆_ℎ(𝑡), 𝐼_ℎ(𝑡), and 𝑅_ℎ(𝑡) denote the number of susceptible, infective, and recovered in the human population, and 𝑆_𝑣(𝑡), 𝐼_𝑣(𝑡) denote the number of susceptible and infective in the mosquito vector population. In addition, 𝑆_ℎ+𝐼_ℎ+𝑅_ℎ=𝑁_ℎ and 𝑆_𝑣+𝐼_𝑣=𝑁_𝑣.

System (1) satisfies the initial conditions: 𝑆_ℎ(𝜃)=𝑆_ℎ0(𝜃), 𝐼_ℎ(𝜃)=𝐼_ℎ0(𝜃), 𝑅_ℎ(𝜃)=𝑅_ℎ0(𝜃), 𝑆_𝑣(𝜃)=𝑆_𝑣0(𝜃), 𝐼_𝑣(𝜃)=𝐼_𝑣0(𝜃), 𝜃∈[−𝜏,0], where 𝜏=max{𝜏₁,𝜏₂}. The nonnegative cone of 𝑅⁶ space 𝑆_ℎ𝐼_ℎ𝑅_ℎ𝑆_𝑣𝐼_𝑣𝑁_ℎ is positively invariant for system (1) since the vector field on the boundary does not point to the exterior. What is more, since 𝑑𝑁_ℎ/𝑑𝑡<0 for 𝑁_ℎ>Λ_ℎ/𝜇_ℎ and 𝑁_𝑣 is constant, all trajectories in the first quadrant enter or stay inside the region 𝑍₊=𝑆_ℎ+𝐼_ℎ+𝑅_ℎ=𝑁_ℎ⩽Λ_ℎ𝜇_ℎ,𝑆_𝑣+𝐼_𝑣=𝑁_𝑣. (2) The continuity of the right-hand side of (1) implies that unique solutions exist on a maximal interval. Since solutions approach, enter or stay in 𝐷₊, they are eventually bounded and hence exist for 𝑡>0. Therefore, the initial value problem for system (1) is mathematically well posed and biologically reasonable since all variables remain nonnegative.

In order to reduce the number of parameters and simplify system (1), we normalize the human and mosquito vector population 𝑠_ℎ=𝑆_ℎ/(Λ_ℎ/𝜇_ℎ), 𝑖_ℎ=𝐼_ℎ/(Λ_ℎ/𝜇_ℎ), 𝑟_ℎ=𝑅_ℎ/(Λ_ℎ/𝜇_ℎ), 𝑛_ℎ=𝑁_ℎ/(Λ_ℎ/𝜇_ℎ), 𝑠_𝑣= 𝑆_𝑣/𝑁_𝑣, 𝑖_𝑣=𝐼_𝑣/𝑁_𝑣, and 𝑚=𝑁_𝑣/(Λ_ℎ/𝜇_ℎ).

Since 𝑟_ℎ=𝑛_ℎ−𝑠_ℎ−𝑖_ℎ and 𝑠_𝑣=1−𝑖_𝑣, the dynamics of system (1) is qualitatively equivalent to the dynamics of system given by 𝑑𝑠_ℎ𝑑𝑡=𝜇_ℎ−𝑏𝛽_ℎ𝑚𝑛_ℎ𝑡−𝜏₁𝑖_𝑣𝑡−𝜏₁𝑠_ℎ𝑡−𝜏₁𝑒^{−𝜇ℎ𝜏1}−𝜇_ℎ𝑠_ℎ+𝜃𝑖_ℎ+𝜌𝑛_ℎ−𝑠_ℎ−𝑖_ℎ,𝑑𝑖_ℎ=𝑑𝑡𝑏𝛽_ℎ𝑚𝑛_ℎ𝑡−𝜏₁𝑖_𝑣𝑡−𝜏₁𝑠_ℎ𝑡−𝜏₁𝑒^{−𝜇ℎ𝜏1}−𝛾_ℎ+𝜇_ℎ+𝛼_ℎ+𝜃𝑖_ℎ,𝑑𝑖_𝑣=𝑑𝑡𝑏𝛽_𝑣𝑛_ℎ𝑡−𝜏₂𝑖_ℎ𝑡−𝜏₂1−𝑖_𝑣𝑡−𝜏₂𝑒^{−𝜇𝑣𝜏2}−𝜇_𝑣𝑖_𝑣,𝑑𝑛_ℎ𝑑𝑡=𝜇_ℎ−𝜇_ℎ𝑛_ℎ−𝛼_ℎ𝑖_ℎ, (3) and all trajectories in the nonnegative cone 𝑅⁴₊ enter or stay inside the region 𝐷={0⩽𝑠_ℎ, 0⩽𝑖_ℎ, 𝑠_ℎ+𝑖_ℎ⩽𝑛_ℎ⩽1, 0⩽𝑖_𝑣⩽1}.

Define the basic reproduction number by 𝑅₀=𝑚𝑏²𝛽_ℎ𝛽_𝑣𝑒^{−𝜇ℎ𝜏1}𝑒^{−𝜇𝑣𝜏2}𝜇_𝑣𝛾_ℎ+𝜇_ℎ+𝛼_ℎ.+𝜃 (4) It represents the expected number of secondary cases caused by a single infected individual introduced into an otherwise susceptible population of hosts and vectors. Take a primary case of host with a natural death rate of 𝜇_ℎ, the average time spent in an infectious state is 1/(𝛾_ℎ+𝜇_ℎ+𝛼_ℎ+𝜃). Since the incubation period in mosquitoes has duration 𝜏₂, during which period some mosquitoes may die, the average number of mosquito bites received from 𝑚 susceptible mosquitoes each with a biting rate 𝑏 gives a total of 𝑏𝛽_𝑣𝑚𝑒^{−𝜇𝑣𝜏2}/(𝛾_ℎ+𝜇_ℎ+𝛼_ℎ+𝜃) mosquitoes successfully infected by the primary human case. Each of these mosquitoes survives for an average time 1/𝜇_𝑣. Because of another incubation period 𝜏₁ in human, during which period some individuals may die, totally 𝑏𝛽_ℎ𝑒^{−𝜇ℎ𝜏1}/𝜇_𝑣 hosts will be infected by a single mosquito. Therefore, the total number of secondary cases is thus 𝑚𝑏²𝛽_ℎ𝛽_𝑣𝑒^{−𝜇ℎ𝜏1}𝑒^{−𝜇𝑣𝜏2}/ 𝜇_𝑣(𝛾_ℎ+𝜇_ℎ+𝛼_ℎ+𝜃), which is the basic reproduction number 𝑅₀.

3. Existence of Equilibria

Equating the derivatives on the left-hand side to zero and solving the resulting algebraic equations. The points of equilibrium 𝑃(𝑠_ℎ̂,𝑖_ℎ̂,𝑖_𝑣̂,𝑛_ℎ) satisfy the following relations: 𝑠_ℎ̂=𝜇_ℎ+𝜌−𝐻𝑖_ℎ̂𝜇_ℎ,𝑛+𝜌_ℎ=𝜇_ℎ−𝛼_ℎ𝑖_ℎ̂𝜇_ℎ,𝑖_𝑣̂=𝑏𝛽_𝑣𝜇_ℎ𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2}𝑏𝛽_𝑣𝜇_ℎ𝑒^{−𝜇𝑣𝜏2}−𝛼_ℎ𝜇_𝑣𝑖_ℎ̂+𝜇_𝑣𝜇_ℎ, (5) where 𝐻=𝛾_ℎ+𝜇_ℎ+𝛼_ℎ+𝜌+𝜌(𝛼_ℎ/𝜇_ℎ). Substituting (5) in the corresponding second equilibrium equation of (3), we obtain that the solutions are 𝑖_ℎ̂=0, and the zero points of the quadratic polynomial 𝑟𝑖_ℎ=𝑎₁𝑖²_ℎ+𝑎₂𝑖_ℎ+𝑎₃, (6) where 𝑎₁=𝑏𝛽_𝑣𝜇_ℎ𝑒^{−𝜇𝑣𝜏2}−𝛼_ℎ𝜇_𝑣𝛼_ℎ𝜇_ℎ,𝑎₂=2𝛼_ℎ𝜇_𝑣−𝑏𝛽_𝑣𝜇_ℎ𝑒^{−𝜇𝑣𝜏2}−𝑚𝑏²𝛽_ℎ𝛽_𝑣𝜇_ℎ𝐻𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝛼_ℎ𝑒+𝜃^{−𝜇ℎ𝜏1}𝑒^{−𝜇𝑣𝜏2},𝑎₃=𝜇_ℎ𝜇_𝑣𝑅₀,𝑅−1₀=𝑚𝑏²𝛽_ℎ𝛽_𝑣𝑒^{−𝜇ℎ𝜏1}𝑒^{−𝜇𝑣𝜏2}𝜇_𝑣𝛾_ℎ+𝜇_ℎ+𝛼_ℎ.+𝜃 (7)

The solution 𝑖_ℎ̂=0 gives the disease-free equilibrium point 𝑃₀(1,0,0,1). We are looking for nontrivial equilibrium solutions in the interior of Ω. From (5), in order to keep 𝑠_ℎ̂ and 𝑛_ℎ positive, 𝑖_ℎ̂∈(0,(𝜇_ℎ+𝜌)/𝐻) must be satisfied. Evaluating 𝑟(𝑖_ℎ) at the end points of the interval, we obtain 𝑟(0)=𝜇_ℎ𝜇_𝑣𝑅₀,𝑟𝜇−1_ℎ+𝜌𝐻<0. (8)

When 𝑅₀>1, then 𝑟(0)>0; therefore, there exists a unique root in the interval (0,(𝜇_ℎ+𝜌)/𝐻), which implies the existence of a unique equilibrium point 𝑃^∗(𝑠^∗_ℎ̂,𝑖^∗_ℎ̂,𝑖^∗_𝑣̂,𝑛^∗_ℎ), where 𝑖^∗_ℎ̂ satisfies: if 𝑎₁=0, 𝑖^∗_ℎ̂=−𝑎₃𝑎₂. (9) If 𝑎₁≠0, 𝑖^∗_ℎ̂=−𝑎₂−√Δ2𝑎₁, (10) where Δ=𝑎²₂−4𝑎₁𝑎₃.

When 𝑅₀=1, the roots of (6) are 0 and −𝑎₂/𝑎₁. It is easy to see that there exists a unique root in the interval (0,(𝜇_ℎ+𝜌)/𝐻) if and only if 𝑎₂>0, which implies the existence of a unique equilibrium point 𝑃^∗.

When 𝑅₀<1, then 𝑟(0)<0. The conditions to have at least one root in the mentioned interval are 𝑎₁<0, 0<−𝑎₂/2𝑎₁<(𝜇_ℎ+𝜌)/𝐻 and Δ=𝑎²₂−4𝑎₁𝑎₃⩾0. If Δ>0, there exist two roots in the interval mentioned above, which implies the existence of two equilibria 𝑃₁(𝑠̂_ℎ1,𝑖̂_ℎ1,𝑖̂_𝑣1,𝑛_ℎ1) and 𝑃₂(𝑠̂_ℎ2,𝑖̂_ℎ2,𝑖̂_𝑣2,𝑛_ℎ2), where 𝑖̂_ℎ1=−𝑎₂+√Δ2𝑎₁,𝑖̂_ℎ2=−𝑎₂−√Δ2𝑎₁. (11)

If Δ=0, there is a unique root in the interval mentioned above, which implies the existence of unique equilibrium.

Then, we can conclude the above results in the following theorem.

Theorem 1. (1) The disease-free equilibrium 𝑃₀(1,0,0,1) always exists.
(2) If 𝑅₀>1, there exists a unique positive equilibrium 𝑃^∗(𝑠^∗_ℎ̂,𝑖^∗_ℎ̂,𝑖^∗_𝑣̂,𝑛^∗_ℎ).
(3) If 𝑅₀=1, then there is a positive equilibrium 𝑃^∗(𝑠^∗_ℎ̂,𝑖^∗_ℎ̂,𝑖^∗_𝑣̂,𝑛^∗_ℎ) when 𝑎₂>0; otherwise, there is no positive equilibrium.
(4) If 𝑅₀<1, then (a) if 𝑎₁⩾0, there is no positive equilibrium; (b) if 𝑎₁<0, the system (3) has two positive equilibria 𝑃₁(𝑠̂_ℎ1,𝑖̂_ℎ1,𝑖̂_𝑣1,𝑛_ℎ1) and 𝑃₂(𝑠̂_ℎ2,𝑖̂_ℎ2,𝑖̂_𝑣2,𝑛_ℎ2) if and only if Δ>0 and 0<−𝑎₂/2𝑎₁<(𝜇_ℎ+𝜌)/𝐻, and these two equilibria merge with each other if and only if 0<−𝑎₂/2𝑎₁<(𝜇_ℎ+𝜌)/𝐻 and Δ=0; otherwise, there is no positive equilibrium.

4. Disease-Free Equilibrium

It is easy to see the disease-free equilibrium 𝑃₀(1,0,0,1) always exists. In this section, we study the stability of the disease-free equilibrium 𝑃₀.

Theorem 2. 𝑃₀ is locally stable if 𝑅₀<1. 𝑃₀ is a degenerate equilibrium of codimension one and is stable except in one direction if 𝑅₀=1. 𝑃₀ is unstable if 𝑅₀>1.

Proof. Linearizing the system (3) at 𝑃₀. The eigenvalues of the Jacobian matrix 𝐽(𝑃₀) are −𝜇_ℎ−𝜌, −𝜇_ℎ and the solutions of the following transcendental equation: 𝜆+𝜇_𝑣𝛾_ℎ+𝜇_ℎ+𝛼_ℎ+𝜃−𝑚𝑏²𝛽_ℎ𝛽_𝑣𝑒^{−𝜇ℎ𝜏1−𝜇𝑣𝜏2}𝑒^{−(𝜏1+𝜏2)𝜆}=0. (12) Let 𝐹𝜆,𝜏₁,𝜏₂=𝜆+𝜇_𝑣𝛾_ℎ+𝜇_ℎ+𝛼_ℎ+𝜃−𝑚𝑏²𝛽_ℎ𝛽_𝑣𝑒^{−𝜇ℎ𝜏1−𝜇𝑣𝜏2}𝑒^{−(𝜏1+𝜏2)𝜆}. (13) Clearly, 𝐹(𝜆,𝜏₁,𝜏₂) is an analytic function. 𝐹(0,𝜏₁,𝜏₂)=𝜇_𝑣(𝛾_ℎ+𝜇_ℎ+𝛼_ℎ+𝜃)(1−𝑅₀), and 𝐹(𝜆,0,0)= (𝜆+𝜇_𝑣)(𝛾_ℎ+𝜇_ℎ+𝛼_ℎ+𝜃)−𝑚𝑏²𝛽_ℎ𝛽_𝑣. In the following, we discuss the distribution of the solutions of (12) in three cases.(i)𝑅₀<1, then 𝐹(0,𝜏₁,𝜏₂)>0. Since 𝐹^′_𝜆(𝜆,𝜏₁,𝜏₂)>0 for 𝜆≥0, 𝜏₁>0 and 𝜏₂>0, (12) has no zero root and positive real roots for all positive 𝜏₁ and 𝜏₂. Now we claim that (12) does not have any purely imaginary roots. Suppose that (12) has a pair of purely imaginary roots 𝜔𝑖, 𝜔>0 for some 𝜏₁ and 𝜏₂. Then, by calculation, 𝜔 must be a positive real root of 𝜔⁴+𝛾_ℎ+𝜇_ℎ+𝛼_ℎ+𝜃²+𝜇²_𝑣𝜔²+𝛾_ℎ+𝜇_ℎ+𝛼_ℎ+𝜃²𝜇²_𝑣−𝑚𝑏²𝛽_ℎ𝛽_𝑣𝑒^{−𝜇ℎ𝜏1−𝜇𝑣𝜏22}=0. (14) However, it is easy to see that (14) does not have nonnegative real roots when 𝑅₀<1. Hence, (12) does not have any purely imaginary roots. On the other hand, one can easily get that the roots of 𝐹(𝜆,0,0)=0 all have negative real parts when 𝑅₀<1. By the implicit function theorem and the continuity of 𝐹(𝜆,𝜏₁,𝜏₂), we know that all roots of (12) have negative real parts for positive 𝜏₁ and 𝜏₂, which implies that 𝑃₀ is stable.(ii)𝑅₀=1, then 𝐹(0,𝜏₁,𝜏₂)=0. Since 𝐹^′_𝜆(𝜆,𝜏₁,𝜏₂)>0 for 𝜆≥0, 𝜏₁>0 and 𝜏₂>0, (12) has a simple zero root and no positive root for all positive 𝜏₁ and 𝜏₂. Using a similar argument as in (i), we can obtain that except a zero root, all roots of (12) have negative real parts for positive 𝜏₁ and 𝜏₂. Thus, 𝑃₀ is a degenerate equilibrium of codimension one and is stable except in one direction.(iii)𝑅₀>1, then 𝐹(0,𝜏₁,𝜏₂)<0. Since lim_𝜆→∞𝐹(𝜆,𝜏₁,𝜏₂)=∞ and 𝐹^′_𝜆(𝜆,𝜏₁,𝜏₂)>0 for 𝜆≥0, 𝜏₁>0 and 𝜏₂>0, (12) has a unique positive real root for all positive 𝜏₁ and 𝜏₂ and 𝑃₀ is unstable.

Theorem 3. If 𝑅₀≤1, the disease-free equilibrium 𝑃₀ is globally asymptotically stable when 𝛼=0.

Proof. We denote by 𝑥_𝑡 the translation of the solution of the system (3), that is, 𝑥_𝑡=(𝑠_ℎ(𝑡+𝜃), 𝑖_ℎ(𝑡+𝜃), 𝑖_𝑣(𝑡+𝜃)), 𝑛_ℎ(𝑡+𝜃) where 𝜃∈[−𝜏,0], 𝜏=max{𝜏₁,𝜏₂}. In order to prove the globally stability of 𝑃₀, we define a Lyapunov function: 𝐿𝑥_𝑡=𝑒^𝜇ℎ𝜏1𝑖_ℎ+𝑚𝑏𝛽_ℎ𝜇_𝑣𝑒^{𝜇𝑣𝜏2}𝑖_𝑣(𝑡)+𝑚𝑏𝛽_ℎ^𝑡_𝑡−𝜏1𝑖_𝑣(𝑡)𝑠_ℎ(𝑡)𝑛_ℎ(𝑡)+𝑑𝑡𝑚𝑏²𝛽_ℎ𝛽_𝑣𝜇_𝑣^𝑡_𝑡−𝜏2𝑖_ℎ(𝑡)1−𝑖_𝑣(𝑡)𝑛_ℎ(𝑡)𝑑𝑡. (15) The derivative of 𝐿 along the solutions of (3) is given by 𝑑𝐿𝑥_𝑡𝜇𝑑𝑡=−_ℎ+𝛾_ℎ+𝜃𝑒^𝜇ℎ𝜏1𝑖_ℎ(𝑡)1−𝑅₀𝑒^{−𝜇ℎ𝜏1}1−𝑖_𝑣(𝑡)𝑛_ℎ(𝑡)−𝑚𝑏𝛽_ℎ𝑒^{𝜇𝑣𝜏2}𝑖_𝑣(𝑡)𝑠1−^ℎ(𝑡)𝑛_ℎ(𝑡)𝑒^{−𝜇𝑣𝜏2}. (16) Since lim_𝑡→∞𝑛_ℎ=1 and 𝑅₀≤1, we can easily get that 𝑑𝐹/𝑑𝑡≤0 and the set {𝑃₀} is the largest invariant set within the set where 𝑑𝐿(𝑥_𝑡)/𝑑𝑡=0. By LaSalle-Lyapunov Theorem, all trajectories approach 𝑃₀ as 𝑡→∞.

Usually, the disease can be controlled if the basic reproduction number is smaller than one. Nevertheless, Theorem 1 indicates the possibility of backward bifurcation (where the locally asymptotically stable DFE coexists with endemic equilibrium when 𝑅₀<1) when 𝑎₁<0 for the model (3) and it becomes impossible to control the disease by just reducing the basic reproduction number below one. Therefore a further threshold condition beyond the basic reproduction is essential for the control of the spread of malaria.

To check for this, the discriminant Δ is set to zero and solved for the critical value of 𝑅₀, denoted by 𝑅₀, given by 𝑅₀𝑎=1+²²4𝑎₁𝜇_ℎ𝜇_𝑣. (17) Thus, backward bifurcation would occur for values of 𝑅₀ such that 𝑅₀<𝑅₀<1. This is illustrated by simulating the model with the following set of parameter values (it should be stated that these parameters are chosen for illustrative purpose only, and may not necessarily be realistic epidemiologically).

Let 𝜇_ℎ=1/22000, 𝜇_𝑣=1/14, 𝛽_ℎ=0.3, 𝛽_𝑣=0.95, 𝑚=0.5, 𝛼_ℎ=0.05, 𝛾_ℎ=0.26, 𝑏=0.5, 𝜌=0.01, 𝜃=0.08, 𝜏₁=15. With this set of parameters, 𝑅₀=0.093970<1. When 𝜏₂=10, 𝑅₀=0.62555<1, so that 𝑅₀<𝑅₀<1. There are two positive roots for (6), which are 𝑖̂_ℎ1=.3501527770𝑒−3 and 𝑖̂_ℎ2=.8864197396𝑒−3.

Thus, the following result is established.

Lemma 4. The model (3) undergoes backward bifurcation when 𝑎₁<0 and 0<−𝑎₂/2𝑎₁<(𝜇_ℎ+𝜌)/𝐻 holds and 𝑅₀<𝑅₀<1.

Epidemiologically, because of the existence of backward bifurcation, whether malaria will prevail or not depends on the initial states. 𝑅₀<1 is no longer sufficient for disease elimination. In order to eradicate the disease, we have a subthreshold number 𝑅₀ and it is now necessary to reduce 𝑅₀ to a value less than 𝑅₀.

The existence of the subthreshold condition 𝑅₀ also has some implications. In some areas where malaria dies out (𝑅₀<𝑅₀), it is possible for the disease to reestablish itself in the population because of a small change in environmental or control variables which increase the basic reproduction number above. On the other hand, in some areas with endemic malaria, it may be possible to eradicate the disease with small increasing in control programs such that the basic reproduction number less than 𝑅₀.

From the expression of 𝑅₀, we can see that the basic reproduction number is a decreasing function of both time delays. The time delays have important effect on the basic reproduction number. Especially, since global warming decreases the duration of incubation period and 𝑅₀ increases with the decreasing of 𝜏₁ or 𝜏₂, global warming will affect the transmission of malaria a lot.

From Theorem 3, if there is no disease-induced death, backward bifurcation will not occur and there is a unique endemic equilibrium when 𝑅₀>1. Therefore, the disease-induced death is the cause of the backward bifurcation phenomenon in the model (3).

5. Endemic Equilibria

To obtain precise results, we first assume that malaria does not produce significant mortality (𝛼=0). The assumption is not justifiable in all regions where malaria is endemic but it is a useful first approximation. In this part, we will determine the stability of the unique endemic equilibria 𝑃^∗ when 𝑅₀>1.

Linearizing the system (3) at 𝑃^∗. The eigenvalues of the Jacobian matrix 𝐽(𝑃^∗) are −𝜇_ℎ and the solutions of the following transcendental equation: 𝜆³+𝐴₁𝜆²+𝐴₂𝜆+𝐴₃=0, (18) where 𝐴₁=2𝜇_ℎ+𝛾_ℎ+𝜌+𝜃+𝜇_𝑣+𝑚𝑏𝛽_ℎ𝑖_𝑣̂𝑒^{−𝜇ℎ𝜏1}𝑒^−𝜆𝜏1+𝑏𝛽_𝑣𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2}𝑒^−𝜆𝜏2,𝐴₂=𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃+𝜇_ℎ+𝜌𝜇_𝑣+𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣+𝜇_ℎ+𝜌+𝛾_ℎ+𝜇_ℎ+𝜃𝑏𝛽_𝑣𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2}𝑒^−𝜆𝜏2−𝜇_𝑣𝛾_ℎ+𝜇_ℎ+𝜃𝑒^{−(𝜏1+𝜏2)𝜆}+𝛾_ℎ+𝜇_ℎ+𝜌+𝜇_𝑣+𝑏𝛽_𝑣𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2}𝑒^−𝜆𝜏2×𝑚𝑏𝛽_ℎ𝑖_𝑣̂𝑒^{−𝜇ℎ𝜏1}𝑒^−𝜆𝜏1,𝐴₃=𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣+𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃𝑏𝛽_𝑣𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2}𝑒^−𝜆𝜏2−𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣𝑒^{−(𝜏1+𝜏2)𝜆}+𝛾_ℎ+𝜇_ℎ+𝜌𝑚𝑏𝛽_ℎ𝑖_𝑣̂𝑒^{−𝜇ℎ𝜏1}𝑒^−𝜆𝜏1×𝜇_𝑣+𝑏𝛽_𝑣𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2}𝑒^−𝜆𝜏2. (19) For any nonnegative 𝜏₁ and 𝜏₂, we have the following proposition.

Proposition 5. For the endemic equilibrium 𝑃^∗(𝑠^∗_ℎ̂,𝑖^∗_ℎ̂,𝑖^∗_𝑣̂,𝑛^∗_ℎ) of the system with characteristic equation (18), one always has 𝐴₁>0,𝐴₂>0,𝐴₃>0,𝐴₁𝐴₂−𝐴₃>0. (20)

Proof. It is clear that 𝐴₁>0,  𝐴₂>0 and 𝐴₃>0. For 𝐴₁𝐴₂−𝐴₃, we have 𝐴₁𝐴₂−𝐴₃=𝛾_ℎ+𝜇_ℎ×𝜇+𝜌_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜇+𝜃_ℎ+𝜌𝜇_𝑣+𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣+𝜇_ℎ+𝜌+𝛾_ℎ+𝜇_ℎ+𝜃𝑏𝛽_𝑣𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2}𝑒^−𝜆𝜏2−𝜇_𝑣𝛾_ℎ+𝜇_ℎ+𝜃𝑒^{−(𝜏1+𝜏2)𝜆}+𝛾_ℎ+𝜇_ℎ+𝜌+𝜇_𝑣𝑚𝑏𝛽_ℎ𝑖_𝑣̂𝑒^{−𝜇ℎ𝜏1}𝑒^−𝜆𝜏1+𝑚𝑏𝛽_ℎ𝑖_𝑣̂𝑒^{−𝜇ℎ𝜏1}𝑒^−𝜆𝜏1+𝜇_ℎ+𝜃𝐴₂+𝜇_𝑣+𝑏𝛽_𝑣𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2}𝑒^−𝜆𝜏2×𝜇_ℎ+𝜌𝜇_𝑣+𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣+𝜇_ℎ+𝜌+𝛾_ℎ+𝜇_ℎ+𝜃𝑏𝛽_𝑣𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2}𝑒^−𝜆𝜏2+𝛾_ℎ+𝜇_ℎ+𝜌+𝜇_𝑣+𝑏𝛽_𝑣𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2}𝑒^−𝜆𝜏2×𝑚𝑏𝛽_ℎ𝑖_𝑣̂𝑒^{−𝜇ℎ𝜏1}𝑒^−𝜆𝜏1−𝜇_𝑣𝛾_ℎ+𝜇_ℎ+𝜃𝑒^{−(𝜏1+𝜏2)𝜆}+𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣𝑒^{−(𝜏1+𝜏2)𝜆}>0. (21)

Case 1. When 𝜏₁=𝜏₂=0, as a result of Proposition 5 and Hurwitz criterion, all roots of the characteristic equation (18) have negative real parts and the endemic equilibrium 𝑃^∗ of (3) is stable when 𝜏₁=𝜏₂=0.

Case 2. When 𝜏₁>0,  𝜏₂=0, the characteristic equation (18) becomes 𝜆³+𝐴₁₁𝜆²+𝐴₂₁𝜆+𝐴₃₁=𝑒^−𝜆𝜏1𝑇₁₁𝜆²+𝑇₂₁𝜆+𝑇₃₁, (22) where 𝐴₁₁=2𝜇_ℎ+𝛾_ℎ+𝜌+𝜃+𝜇_𝑣+𝑏𝛽_𝑣𝑖_ℎ̂,𝐴₂₁=𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃+𝜇_ℎ+𝜌𝜇_𝑣+𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣+𝜇_ℎ+𝜌+𝛾_ℎ+𝜇_ℎ+𝜃𝑏𝛽_𝑣𝑖_ℎ̂,𝐴₃₁=𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣+𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃𝑏𝛽_𝑣𝑖_ℎ̂,𝑇₁₁=−𝑚𝑏𝛽_ℎ𝑖_𝑣̂𝑒^{−𝜇ℎ𝜏1},𝑇₂₁=𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣−𝛾_ℎ+𝜇_ℎ+𝜌+𝜇_𝑣+𝑏𝛽_𝑣𝑖_ℎ̂𝑚𝑏𝛽_ℎ𝑖_𝑣̂𝑒^{−𝜇ℎ𝜏1},𝑇₃₁=𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣−𝛾_ℎ+𝜇_ℎ+𝜌𝜇_𝑣+𝑏𝛽_𝑣𝑖_ℎ̂𝑚𝑏𝛽_ℎ𝑖_𝑣̂𝑒^{−𝜇ℎ𝜏1}. (23) By the implicit function theorem and the continuity of the left-hand side function of (18), all roots of (22) have negative real parts for small 𝜏₁. Notice that the condition 𝑅₀>1 is equivalent to 𝜏₁<𝜏^∗₁=1𝜇_ℎln𝑚𝑏²𝛽_ℎ𝛽_𝑣𝜇_𝑣𝛾_ℎ+𝜇_ℎ.+𝜃 (24) Furthermore, we claim that (22) does not have any nonnegative real roots for any 𝜏₁>0. Rewrite (22) by moving the positive terms from the right-hand side to the left-hand side. The rewritten (22) takes the form 𝜆³+𝐴₁₁𝜆²+𝐴^∼₂₁𝜆+𝐴^∼₃₁=𝑒^−𝜆𝜏1𝑇₁₁𝜆²+𝑇^∼₂₁𝜆+𝑇^∼₃₁, (25) where 𝐴^∼₂₁=𝐴₂₁−(𝛾_ℎ+𝜇_ℎ+𝜃)𝜇_𝑣𝑒^−𝜆𝜏1, 𝐴^∼₃₁=𝐴₃₁−(𝜇_ℎ+𝜌)(𝛾_ℎ+𝜇_ℎ+𝜃)𝜇_𝑣𝑒^−𝜆𝜏1. It is easy to see 𝐴^∼₂₁>0 and 𝐴^∼₃₁>0 for all 𝜆≥0 and 𝜏₁∈(0,𝜏^∗₁). Consequently, the left-hand side in (25) is positive for all 𝜆≥0 while the right-hand side is negative for all 𝜆≥0 and the two cannot be equal for any 𝜆≥0. Therefore, (22) does not have any nonnegative real roots for any 𝜏₁∈(0,𝜏^∗₁). Now we want to show that all roots of (22) have negative real parts for all 𝜏₁∈(0,𝜏^∗₁). To do so, we show that (22) does not have any purely imaginary roots for all 𝜏₁∈(0,𝜏^∗₁).
We assume that 𝜆=𝑖𝜔 with 𝜔>0 being a root of (22). Then 𝜔 must satisfy the following system: 𝐴₂₁𝜔−𝜔³=𝑇₂₁𝜔cos𝜔𝜏₁−𝑇₃₁−𝑇₁₁𝜔²sin𝜔𝜏₁,𝐴₃₁−𝐴₁₁𝜔²=𝑇₃₁−𝑇₁₁𝜔²cos𝜔𝜏₁+𝑇₂₁𝜔sin𝜔𝜏₁. (26) Thus, 𝜔 must be a positive root of 𝜔⁶+𝐵₁𝜔⁴+𝐵₂𝜔²+𝐵₃=0, (27) where 𝐵₁=𝐴²₁₁−2𝐴₂₁−𝑇²₁₁,𝐵₂=𝐴₂₁−2𝐴₁₁𝐴₃₁+2𝑇₁₁𝑇₃₁−𝑇²₂₁,𝐵₃=𝐴²₃₁−𝑇²₃₁. (28) Let 𝑧=𝜔², then (27) becomes 𝑧³+𝐵₁𝑧²+𝐵₂𝑧+𝐵₃=0. (29)
Clearly, If 𝐵₁≥0, 𝐵₂≥0 and 𝐵₃≥0, then (29) has no positive real roots. Therefore, (22) does not have any purely imaginary roots for all 𝜏₁∈(0,𝜏^∗₁) so that all roots of the characteristic equation (22) have negative real parts and the endemic equilibrium 𝑃^∗ of (3) is stable.

Case 3. When 𝜏₂>0,  𝜏₁=0, the characteristic equation (18) becomes 𝜆³+𝐴₁₂𝜆²+𝐴₂₂𝜆+𝐴₃₂=𝑒^−𝜆𝜏2𝑇₁₂𝜆²+𝑇₂₂𝜆+𝑇₃₂, (30) where 𝐴₁₂=2𝜇_ℎ+𝛾_ℎ+𝜌+𝜃+𝜇_𝑣+𝑚𝑏𝛽_ℎ𝑖_𝑣̂,𝐴₂₂=𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃+𝜇_ℎ+𝜌𝜇_𝑣+𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣+𝜇_ℎ+𝜌+𝛾_ℎ+𝜇_𝑣𝑚𝑏𝛽_ℎ𝑖_𝑣̂,𝐴₃₂=𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣+𝜇_𝑣𝛾_ℎ+𝜇_ℎ+𝜌𝑚𝑏𝛽_ℎ𝑖_𝑣̂,𝑇₁₂=−𝑏𝛽_𝑣𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2},𝑇₂₂=𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣−2𝜇_ℎ+𝛾_ℎ+𝜌+𝜃+𝑚𝑏𝛽_ℎ𝑖_𝑣̂𝑏𝛽_𝑣𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2},𝑇₃₂=𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃𝜇_𝑣−𝛾_ℎ+𝜇_ℎ+𝜌𝑚𝑏²𝛽_ℎ𝛽_𝑣𝑖_ℎ̂𝑖_𝑣̂𝑒^{−𝜇𝑣𝜏2}−𝜇_ℎ+𝜌𝛾_ℎ+𝜇_ℎ+𝜃𝑏𝛽_𝑣𝑖_ℎ̂𝑒^{−𝜇𝑣𝜏2}. (31) The condition 𝑅₀>1 is equivalent to 𝜏₂<𝜏^∗₂=1𝜇_𝑣ln𝑚𝑏²𝛽_ℎ𝛽_𝑣𝜇_𝑣𝛾_ℎ+𝜇_ℎ.+𝜃 (32) Using a similar argument as in Case 2, we know that all roots of (30) have negative real parts for 𝜏₂∈(0,𝜏^∗₂) when 𝐶₁≥0,  𝐶₂≥0, and 𝐶₃≥0, where 𝐶₁=𝐴²₁₂−2𝐴₂₂−𝑇²₁₂,𝐶₂=𝐴₂₂−2𝐴₁₂𝐴₃₂+2𝑇₁₂𝑇₃₂−𝑇²₂₂,𝐶₃=𝐴²₃₂−𝑇²₃₂. (33)

Case 4. When 𝜏₂>0,  𝜏₁>0, the condition 𝑅₀>1 is equivalent to 𝜏₂<𝜏^∗₂𝜏₁=1𝜇_𝑣ln𝑚𝑏²𝛽_ℎ𝛽_𝑣𝑒^{−𝜇ℎ𝜏1}𝜇_𝑣𝛾_ℎ+𝜇_ℎ.+𝜃 (34) From Cases 1, and 2 the roots of (18) only have negative real parts for 𝜏₁∈[0,𝜏^∗₁) and 𝜏₂=0 under the condition of 𝐵_𝑖≥0, 𝑖=1,2,3. By the implicit function theorem and the continuity of the left-hand side function of (18), there is a 𝜏₂(𝜏₁) satisfying 0<𝜏₂(𝜏₁)≤𝜏^∗₂(𝜏₁), such that all roots of (18) have negative real parts for 0<𝜏₂<𝜏₂(𝜏₁). We show that 𝜏₂(𝜏₁)=𝜏^∗₂(𝜏₁) when 𝐵_𝑖≥0 and 𝐶_𝑖≥0, 𝑖=1,2,3. Suppose 𝜏₂(𝜏₁)<𝜏^∗₂(𝜏₁) for 𝜏₁∈[0,𝜏^∗₁), then there must be a 𝜏^∼₂(𝜏₁), 𝜏₂(𝜏₁)<𝜏^∼₂(𝜏₁)<𝜏^∗₂(𝜏₁), such that one root of (18) has nonnegative real part for 𝜏₂=𝜏^∼₂(𝜏₁). As a result of the continuity of 𝜏₂ in 𝜏₁, we have 𝜏^∼₂(0)<𝜏^∗₂(0)=𝜏^∗₂. However, from the argument in Case 3, we know that all the roots of (18) have negative real parts for 𝜏₁=0 and 𝜏₂∈(0,𝜏^∗₂). Contradict. Thus, 𝜏₂(𝜏₁)=𝜏^∗₂(𝜏₁) which implies that the endemic equilibrium 𝑃^∗ is stable when 𝜏₁∈[0,𝜏^∗₁), 𝜏₂∈[0,𝜏^∗₂(𝜏₁)), 𝐵_𝑖≥0, and 𝐶_𝑖≥0, 𝑖=1,2,3. The above analysis can be summarized into the following theorem.

Theorem 6. If 𝑅₀>1, 𝐵_𝑖≥0 and 𝐶_𝑖≥0, 𝑖=1,2,3, the unique endemic equilibrium 𝑃^∗ of system (3) is stable.

6. Discussions

In this paper, in order to study the impact of the incubation periods in human and mosquitoes, we formulate a model with two delays for the transmission dynamics of malaria.

Existence of equilibria is obtained under different conditions and their stabilities are analyzed too. We have also identified the basic reproduction number 𝑅₀, which gives the expected number of new infections (in mosquitoes or humans) from one infectious individual (human or mosquito) over the duration of the infectious period, given that all other members of the population are susceptible, in terms of the model parameters. Especially, we have proved mathematically backward bifurcation may occur for 𝑅₀<1 which implies that bringing the basic reproduction number below 1 is not enough to eradicate malaria.

The parameters 𝜏₁ and 𝜏₂ are the incubation periods in human and mosquitoes, which are temperature dependent. With the increasing of temperature in a range, they become shorter. Since the basic reproduction number is a monotone decreasing function of both time delays, global warming will then exacerbate the transmission of malaria.

The reason inducing backward bifurcation in [6] is the mosquito biting rate. While, the reasons inducing the backward bifurcation in our paper are the standard incidence rate and the disease-induced death rate which is big enough. By our results, the disease-free equilibrium is globally stable when 𝑅₀<1 if the disease-induced death is omitted. However, the result about backward bifurcation in this paper has important implications for malaria control, which implies that bringing the basic reproduction number below 1 is not enough to eradicate malaria.

Malaria transmission can be affected by a lot of aspects. In this paper, we are trying to model the impact of temperature increasing on malaria transmission. In some extent, our results can provide a theoretical principle for allocating and using medical health resource reasonably, which can be applied in the practice of malaria prevention and control.

Acknowledgments

H. Wan is supported by NSFC (no. 11201236 and no. 11271196) and the NSF of the Jiangsu Higher Education Committee of China (no. 11KJA110001 and no. 12KJB110012). J.-a. Cui is supported by NSFC (no. 11071011) and Funding Project for Academic Human Resources Development in Institutions of Higher Learning Under the Jurisdiction of Beijing Municipality (no. PHR201107123).