A mathematical model of malaria transmission in a periodic environment

ABSTRACT In this paper, we present a mathematical model of malaria transmission dynamics with age structure for the vector population and a periodic biting rate of female anopheles mosquitoes. The human population is divided into two major categories: the most vulnerable called non-immune and the least vulnerable called semi-immune. By applying the theory of uniform persistence and the Floquet theory with comparison principle, we analyse the stability of the disease-free equilibrium and the behaviour of the model when the basic reproduction ratio is greater than one or less than one. At last, numerical simulations are carried out to illustrate our mathematical results.


Introduction
Malaria is a common and life-threatening infectious disease in many tropical and subtropical areas. It is caused by the Plasmodium parasite which is transmitted by female anopheles mosquitoes while they bite humans for a blood meal for the development of their eggs. During the blood meal, the mosquito injects sporozoites into the blood stream. In few minutes, the sporozoites enter the liver cells where each sporozoite develops into a tissue schizont that contains 10,000 to 30,000 merozoites. After 1-2 weeks, the schizont ruptures and releases the merozoites into the blood stream which then invade the red blood cells. The clinical symptoms of malaria are due to the rupture of the red blood cells and release of the parasites waste and cells debris into the blood stream. Note that human malaria is caused by five different species of Plasmodium: Plasmodium falciparum, Plasmodium malariae, Plasmodium ovale, Plasmodium vivax and Plasmodium knowlesi. However, P. falciparum is the most prevalent in Africa and it causes the highest mortality rate induced by the disease [27]. The biology of the five species of Plasmodium is generally similar and consists of two distinct phases: a sexual stage at the mosquito host and an asexual stage at the human host.
The World Health Organization (WHO) estimated that there were 214 million malaria cases in 2015, resulting in about 438,000 deaths [39]. Moreover, in endemic regions, children under 5, pregnant women and non-immune adults are most at risk of malaria mortality [12]. Indeed, there are currently over 100 countries where there is a risk of malaria transmission, and these are visited by more than 125 million international travellers every year. International travellers to countries with ongoing local malaria transmission arriving from countries with no transmission are at high risk of malaria infection and its consequences because they lack immunity. Migrants from countries with malaria transmission living in malaria-free countries and returning to their home countries to visit friends and relatives are similarly at risk because of waning or lack of immunity. Despite extensive efforts to eradicate it, malaria caused by P. falciparum remains a significant problem. A characteristic of falciparum malaria disease that complicates control efforts is clinical immunity: an immune response that develops with exposure to parasites and provides protection against the clinical symptoms of malaria, despite the presence of parasites [32]. This immunity is not complete and one can lose it and become susceptible after interruption of exposure. Those who have acquired immunity can host and tolerate malaria parasites without developing any clinical symptoms. They may become asymptomatic carriers of parasites and may transmit slightly the parasites to mosquitoes [15,19,34].
In order to reduce the spread of infectious diseases, mathematical models have been proposed to study their dynamics [5,25,40]. Models can provide estimates of underlying parameters of a real-world problem which are difficult or expensive to obtain through experiment or otherwise [28,29]. They can predict whether the associated disease will spread through the population or die out [3,8,18]. It can also estimate the impact of a control measure and provide useful guidelines to public health for further efforts required for disease elimination.
Concerning the mathematical modelling of malaria, significant breakthroughs have been made in the recent years since the first model introduced by Ronald Ross [31]. According to Ross, if the mosquito population can be reduced to below a certain threshold, then malaria can be eradicated. Some years later, Macdonald [22] improved the model of Ross. He showed that reducing the number of mosquitoes has little effect on epidemiology of malaria in areas of intense transmission. Furthermore, Aron and May [1,2] added various characteristics of malaria to the model of Macdonald, such as an incubation period in the mosquito, superinfection and a period of immunity in humans. An important addition to the malaria models was the inclusion of acquired immunity proposed by Dietz et al. [9].
Other reviews on mathematical modelling in malaria include Ngwa et al. [26] and Chitnis et al. [6]. Indeed, in the Ngwa and Shu model, humans follow an SEIRS-like pattern and mosquitoes follow an SEI pattern, similar to that described by Yang [42], but with only one immune class for humans. Humans move from the susceptible to the exposed class at some probability when they come into contact with an infectious mosquito, and then to the infectious class, as in conventional SEIRS models. However, infectious people can then recover with, or without, a gain in immunity, and either return to the susceptible class or move to the recovered class. Moreover, Chitnis et al. extended the model of Ngwa and Shu by assuming that although individuals in the recovered class are immune, in the sense that they do not suffer from serious illness and do not contract clinical malaria, they still have low levels of Plasmodium in their blood stream and can pass the infection to susceptible mosquitoes.
Recent works have shown that the age structure of vector population and the climate effects are very important factors on the dynamics of malaria transmission [11,20,37] because the dynamics of vector population and the biting rate from mosquitoes to humans are greatly influenced by environmental and climatic factors. In addition, based on the susceptibility, the exposedness and the infectivity of human host, Ducrot et al. [10] have developed a deterministic mathematical model for the transmission of malaria with two host types in the human population. The first type called non-immune is supposed to be very vulnerable to malaria because it has never acquired immunity against the disease. The second type called semi-immune is supposed to be non-vulnerable because it has at least once acquired immunity in his life.
This article is an extension of the model studied in [10] in the sense that we consider the life cycle of vector population and the climatic factors on the biting rate of female anopheles mosquitoes [34]. Thus the model is formulated as a non-autonomous system of differential equations. Through rigorous analysis via theories and methods of dynamical systems [16,33,43], we derive the epidemic threshold parameter R 0 , for predicting disease persistence or extinction in periodic environments and we explore the global stability of the equilibrium state and the existence of positive periodic solution under certain conditions. This paper is organized as follows. In Section 2, we formulate the mathematical model. In Section 3, we show that the dynamical properties of the model are completely determined by the basic reproduction ratio R 0 . Computational simulations are provided in Section 4 in order to illustrate our mathematical results. Finally, we conclude in Section 5.

Mathematical model formulation
In the study of malaria transmission, it has been shown that the susceptibility and the infectivity of the human host depend on certain genetic factors. These two clinical states depend on whether the host has lost his immunity or if he has not yet acquired it. Indeed, acquired immunity is a very important factor in the dynamics of malaria transmission. Several studies have shown that humans who have never acquired immunity during their lifetime have a very high risk of succumbing to malaria compared to those who have already acquired it at least once in their lifetime [15,21,23]. This is the main reason why many children under 5 old years die of malaria in endemic areas. In fact, newborns (from an immunized mother) are protected because of the passive transfer of maternal antibodies through the placenta in the first 3-6 months of life. After these first few months, they are vulnerable to clinical episodes of malaria until they have built their own immunity. Thus, to better understand the dynamics of malaria transmission, it is important to divide human hosts according to their immune status. This will help to develop better strategies against the disease.
Hence, using this assumption, we divide the human population into two major types: the non-immune, namely those who have not acquired immunity (the most vulnerable), and the semi-immune, those who have already acquired their immunity at least once in their life (least vulnerable).
The non-immune population is divided into three epidemiological categories representing the state variables: the susceptible class S e , the exposed class E e and the infectious class I e . Similarly, the semi-immune population is divided into four epidemiological categories representing the state variables: the susceptible class S a , the exposed class E a , the infectious class I a and the immune class R a . Humans in the class R a have some immunity to the disease and do not get clinically ill, but they still harbour low levels of parasite in their blood stream and can pass the infection to the susceptible mosquitoes [6].
Mosquito population is also divided into two major stages: the mature stage (those which can fly) and the aquatic stage (eggs, larvae and pupae). The mature stage is divided into three compartments: the susceptible class S v , the exposed class E v and the infectious class I v . The immature stage represented by the compartment L is constituted by eggs, larvae and pupae.
At any time t, the total size of humans, N h (t), and mature mosquitoes, N v (t), is respectively denoted by the following equations: Furthermore, we assume throughout this paper that: (A1) all vector population measures refer to densities of female anopheles mosquitoes, (A2) the mosquitoes bite only humans, with periodic biting rate β(t), (A3) there is no direct transmission (blood transfusion or from mother to baby) of malaria, (A4) all the new recruits are susceptible. There are no immigrations of infectious humans.
Note that in absence of disease, the vector population dynamics is described by the following diagram .
Hence, according to Figure 1 we obtain the following system: Moreover, female mosquitoes lay their eggs on the surface of water or near rivers. However, if there are too many eggs in the oviposition habitat or very few nutrients and water resources, then females choose another site or lay less eggs. In addition, to complete their development, larvae and pupae need water or nutrients [24]. Hence, using the logistic growth [41], the per capita oviposition rate is given by Finally, we get the following system: where the biological meaning of parameters b, r v , d v , d l and K is given in Table 1.

Interaction between humans and mosquitoes
When an infectious mosquito bites a susceptible non-immune, the parasite enters the body of the non-immune with probability c ve and he moves to the exposed class E e . Some time after, he moves from class E e to the infectious class I e with constant rate ν e . The infectious non-immune moves to the immune class R a at rate α e after acquisition of his immunity. The non-immune disappears from the population through natural death rate d h and malaria death rate γ e . In the same way, when an infectious mosquito bites a susceptible semi-immune, the parasite enters the body of the semi-immune with probability c va and he moves to the exposed class E a . Some time after, he moves from class E a to the infectious class I a with constant rate ν a . By immunology memory, immunity of infectious semi-immune might be rapidly restored at rate α a when they begin to be re-exposed to infection. Immune individuals in class R a can lose their immunity at rate β a and go back to susceptible class S a . The infectious semi-immune leaves the population through natural death rate d h and malaria death rate γ a (γ a ≤ γ e ). At any time, the non-immune population receives a new recruitment, p , and the semi-immune population receives a new recruitment, Similarly, when a susceptible mosquito bites an infectious non-immune (resp. semiimmune, immune), it enters the class E v with a probability c ev (resp. c av ,c av ). Some time after, it moves from class E v to the infectious class I v with rate ν v where it stays for life. Mature mosquitoes disappear from the population through natural mortality d v .
Using the standard incidence as in the model of Ngwa and Shu [26], we define respectively the infection incidence from mosquitoes to non-immune, k e (t), from mosquitoes to semi-immune, k a (t), and from humans to mosquitoes, k v (t).
According to the above assumptions, we get the following schematic diagram ( Figure 2) .

Mathematical model
Using the above assumptions and by making a balance of the movements in each class, we obtain the following system: , , with initial conditions: The growth of the whole human population and mature vector is respectively described by the following equations:Ṅ The evolution of the immature mosquitoes is described by the following equation: Hence, the system (2) can be written as follows: where We assume that (A5) the biting rate β(t) is a continuous ω-periodic positive function with ω = 12 months, (A6) all the parameters of the model are positive except the disease-induced death rates, γ e and γ a , which are assumed to be non-negative.

Positivity and boundedness of solutions
Proof: For any positive initial condition φ, the function g(t, φ) is continuous in (t, φ) and Lipschitzian in φ. So, thanks to Theorems 2.2.1 and 2.2.3 of Hale and Verduyn Lunel [13], the system (2) has a unique solution u(t, φ) on its maximal interval [0, t max ) of existence. Moreover, for all t ≥ 0, let Let us suppose that there exists t 1 > 0 such that e(t 1 ) / ∈ R * + and e(t) > 0 for all t ∈ [0, t 1 ). If e(t) = S e (t), then S e (t) > 0 and k e (t) > 0. Thus, from the first equation of system (2), we haveṠ It then follows that which leads to a contradiction.
If e(t) = E e (t), then S e (t) > 0 and k e (t) > 0. Thus, from the second equation of system (2), we haveĖ It then follows that which leads to a contradiction.
If e(t) = I e (t), then E e (t) > 0. Thus, from the third equation of system (2), we havė It then follows that Thus the system (2) is mathematically well-defined over the whole R 11 + . Nevertheless, the region of biological interest is defined by

Lemma 3.2:
The compact is a positively invariant set, which attracts all positive orbits in R 11 + . Moreover, all the solutions are bounded.
Proof: According to Equations (3), (4) and (5), we havė Thus if Moreover, let us consider the following ordinary differential equations: with respective general solutions: By applying the standard comparison theorem, we have for all t ≥ 0, Thus the compact set is positively invariant, and then the solutions are bounded.

Disease-free equilibria
Let us consider the following threshold parameter: This threshold plays an important role in the dynamics of malaria transmission because it allows to regulate the mosquitoes population dynamics.

Proposition 3.3:
The model (2) has a: • trivial disease-free equilibrium • non-trivial disease-free equilibrium Proof: Solving the system at the disease-free equilibrium, we get Replacing S * v by its expression in (6), we have Hence, it yields that Hence, we obtain E + .

Remark 3.1:
In practice, it is very difficult to find regions completely devoid of anopheles. Hence, we only consider the non-trivial disease-free equilibrium because it is more biologically realistic. So, in the rest of the paper we assume that κ > 1.

Threshold dynamics
Linearizing the system (2) at the equilibrium state E + , we obtain the following system (only the equations for the infectious classes): , This system can be rewritten as follows: where For all t ≥ s, let Y(t, s) be the evolution operator of the linear periodic systeṁ That is, for each s ∈ R, the 7 × 7 matrix Y(t, s) satisfies the equation: where I is the 7 × 7 identity matrix. Let C ω be the ordered Banach space of all ω-periodic functions from R to R 7 which is equipped with the maximum norm . and the positive cone Let us suppose φ(s) ∈ C ω is the initial distribution of infectious individuals in this periodic environment, then F(s)φ(s) is the rate of new infections produced by the infected individuals who were introduced at time s, and Y(t, s)F(s)φ(s) represents the distribution of those infected individuals who were newly infected at time s and remains in the infected compartments at time t for t ≥ s. Hence, the distribution of accumulative new infections at time t produced by all those infected individuals φ(s) introduced at the previous time is given by Let L : C ω −→ C ω be the linear operator defined by Then, L is the next infection operator, and the basic reproduction ratio is R 0 := ρ(L), the spectral radius of L.
In order to calculate R 0 , we consider the following linear ω-periodic system: Let W(t, s, λ), t ≥ s, s ∈ R, be the evolution operator of the system (7) on R 7 . It is clear that The following result will be used in our numerical calculation of the basic reproduction ratio. Proof: Ifβ(t) = εβ(t), then the linear system (7) becomeṡ andŴ(ω, 0, λ), the monodromy matrix of the following system: It is easy to remark thatŴ Thus it then follows that Hence,R 0 = εR 0 .
We have proved that the biting rate is an important parameter for the dynamics of malaria transmission. Thus it could be used as a very good strategy to control malaria. Remark 3.2: Letβ = (1/ω) ω 0 β(s) ds be the average number of bites. Then, the basic reproduction number,R 0 , of the associated autonomous system is given by the spectral radius of the matrixFV −1 [35]. Sō We observe that the basic reproduction number,R 0 , depends on the threshold κ. It could also be used to control the transmission of malaria.

Stability of disease-free equilibrium E +
In this part of the paper, we focus on the stability of the disease-free equilibrium E + . For that, we first study the global behaviour of the model (1). Let

Theorem 3.6: If κ > 1, then the equilibrium M 1 is locally asymptotically stable.
Proof: By linearizing the system (8) at the steady state M 1 , we obtain the following system: Let us consider the following matrix: The matrix B has two eigenvalues: It is clear that the eigenvalue λ 1 always is negative. Now, we aim to show that the eigenvalue λ 2 is also negative if κ > 1. Indeed, we have Thus if κ > 1 both λ 1 and λ 2 are negative. Then, from Poincaré-Lyapunov theorem [24], we conclude that the equilibrium M 1 is locally asymptotically stable if κ > 1.

Theorem 3.7: If κ > 1, then the equilibrium M 1 is globally asymptotically stable in , where
Proof: To prove the global stability of M 1 , we consider the following Lyapunov function: where δ 1 and δ 2 are positive constants. Note that • otherwise, calculating the derivative of the function V, we geṫ Denote ., . the scalar product in R 2 . Then, we havė Now, we construct the symmetric matrix S = −M + 1 2 (N + N T ). Let Then, we have Thus the matrix S has two eigenvalues μ 1 and μ 2 defined as follows: Taking then we obtain It then follows that AX,X = SX,X ≤ 0 for allX ∈ R 2 .
is a consequence of the conclusions (i) and (ii) above.

(t) such that e rt v(t) is a solution ofż(t) = N (t)z(t).
Theorem 3.10: The disease-free equilibrium E + is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.

Proof:
Let A(t) be the Jacobian matrix of (2) evaluated at E + . Then, we have E + is locally asymptotically stable if ρ( A 22 (ω)) < 1 and ρ( F−V (ω)) < 1 [33]. We note that A 22 is a constant matrix and its eigenvalues are given by Since κ > 1, then all the eigenvalues λ 1 , λ 2 , λ 3 and λ 4 are negative. It then follows that Thus, thanks to Theorem 3.8, E + is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.

Proof:
If γ e = γ a = 0, we can rewrite (3) and (4) as follows: Thus there exists a period ω such that for all t > ω , It then follows that . Thus from the system (2), we havė , Let us consider the following auxiliary system:

(E e (t), I e (t), E a (t), I a (t), R a (t), E v (t), I v (t)) =
Hence, from the first and fourth equations of system (2), we get It then follows that the equilibrium E + is globally attractive.

Persistence of malaria
Let us consider the following sets: Let u(t, φ) be the unique solution of (2) with initial condition φ, (t) the periodic semiflow generated by periodic system (2) and P : X −→ X the Poincaré map associated with system (2), namely: Proof: For any initial condition ψ ∈ X 0 , by solving the equations of system (2) we derive that where Thus X 0 is positively invariant. Since X is positively invariant and ∂X 0 is relatively closed in X, it yields that ∂X 0 is positively invariant.
Note that from Lemma 3.2, is a compact set which attracts all positive orbits in X that implies that the discrete-time system P : X −→ X is point dissipative. Moreover, ∀n 0 ≥ 1, P n 0 is compact, it then follows that P admits a global attractor in X.
Theorem 3.14: and the system (2) has at least one positive periodic solution.

Numerical simulations and sensitivity analysis
In this section, we simulate the model in order to illustrate our mathematical results. The numerical values of the model are given in Table 1.
Using the periodicity of the function β(t), we can write it in the following form [4,30]: where a 1 and b 1 are constants chosen such that the function β(t) is positive. Figures 3-5 show that malaria is persistent and the system admits at least one positive periodic solution that illustrates our mathematical result of Theorem 3.14. Furthermore, we remark that the number of infected humans is very high (see Figures 3 and 4) that means   that we are in the endemic region. In addition, we notice that the density of infectious nonimmune is critical (see Figure 3b). As we mentioned it above, the non-immune humans are very vulnerable to malaria virus. Hence, it is very important to develop control measures in order to reduce their infectious number. These control measures can consist in fighting against immature mosquitoes proliferation (by reducing κ), avoiding contact between humans and mosquitoes (by reducing the biting rate β(t)), by protecting a proportion of non-immune through vaccination or by reducing the number of female anopheles already present in the area.
Thus, let us assume that after some years, humans become more conscious about malaria mortality and then, develop some personal methods to avoid mosquitoes bites and also  to reduce the mosquito population. By applying these measures, we obtain the following figures (Figures 6-8). Figures 6-8 show that the solution of the system approaches the disease-free equilibrium E + which is globally asymptotically stable; that illustrates the result of our Theorem 3.11. Moreover, Figure 6 shows that the number of infected non-immune is rising in the short time but the disease cannot persist due to the current prevention and control measures. However, the application of these control measures must be intensified over the first 60 months of infection in order to reduce the malaria death rate for the non-immune.  Next, we perform some sensitivity analysis to determine the influence of parameters κ and β(t) on the dynamics of malaria transmission. First, we analyse the impact of the parameter κ.
If we take the precautionary measures in the aquatic stage by reducing the parameter κ, the number of infectious humans and mosquitoes decreases rapidly as Figures 9 and 10 show it. Thus the parameter κ is very important in the dynamics of malaria transmission. Now, let us suppose that the only control measure applied by the humans consists in avoiding the bites of the mosquitoes. Let θ be the efficiency of intervention measures. By  considering this prevention, we useβ(t) = (1 − θ)β(t) to replace β(t) in the model. If θ = 0.85, then we obtain the following results.
By applying this prevention measure, we remark that malaria disappears in the populations as shown in Figures 11 and 12. It must also be noticed that the basic reproduction decreases with this effort. In fact for the biting rateβ(t) = (1 − θ)β(t), we have obtainedR 0 = (1 − θ)R 0 . This numerical result is consistent with our theoretical result of Proposition 3.5.

Conclusion
Based on the model in [10], we have formulated a mathematical model of malaria transmission with periodic biting rate and mosquitoes population structured in the heterogeneous humans host population. The basic reproduction ratios associated to the model has been determined and a new other threshold dynamics, κ, has been determined too. We have shown that κ is an important parameter which allows to control the dynamics of mosquito populations. Thus for the autonomous model, we clearly observe that basic reproduction number increases with κ. It also emerged from our study that the basic reproduction ratio is highly influenced by the biting rate. Then, these parameters could be used to fight against the disease in the populations.
Furthermore, we have determined the biological realistic disease-free equilibrium, E + , of the model and then we have completely investigated its stability by using the Floquet theory. Indeed, we have proved that if the basic reproduction ratio, R 0 , is less than one, then E + is globally asymptotically stable and malaria dies out of the populations. Otherwise, if R 0 is greater than one, the system admits at least one positive periodic solutions and malaria persists in the populations.
However, it must be noticed that our model is limited due to the following reasons: (i) seasonality is not incorporated in the evolution of immature mosquitoes, that is a mathematical simplification, (ii) the immature class is not differentiated (eggs, larvae and pupae).
One can develop a more realistic model by incorporating these above assumptions and also include the investigation of the stability of positive periodic solution. Otherwise, this study could be applied to some regions with realistic data on malaria transmission that would allow to show the impact of climate on the transmission.