A global mathematical model of malaria transmission dynamics with structured mosquito population and temperature variations

Nonlinear Analysis: Real World Applications 53 (2020) 103081

Contents lists available at ScienceDirect

Nonlinear Analysis: Real World Applications www.elsevier.com/locate/nonrwa

A global mathematical model of malaria transmission dynamics with structured mosquito population and temperature variations Bakary Traoré, Ousmane Koutou, Boureima Sangaré ∗ Laboratoire de Mathématiques Informatique et Applications (La.M.I.A), UFR Sciences et Techniques, Université Nazi BONI, 01 BP 1091 Bobo Dsso 01, Burkina Faso

article

info

Article history: Received 16 August 2018 Received in revised form 15 December 2019 Accepted 16 December 2019 Available online xxxx Keywords: Malaria transmission Basic reproduction ratio Vector reproduction ratio Mosquito population Temperature variations Global stability

abstract In this paper, a mathematical model of malaria transmission which takes into account the four distinct mosquito metamorphic stages is presented. The model is formulated thanks to the coupling of two sub-models, namely the model of mosquito population and the model of malaria parasite transmission due to the interaction between mosquitoes and humans. Moreover, considering that climate factors have a great impact on the mosquito life cycle and parasite survival in mosquitoes, we incorporate seasonality in the model by considering some parameters which are periodic functions. Through a rigorous analysis via theories and methods of dynamical systems, we prove that the global behavior of the model depends strongly on two biological insightful quantities : the vector reproduction ratio Rv and the basic reproduction ratio R0 . Indeed, if Rv < 1, mosquitoes and malaria die out; if Rv > 1 and R0 < 1, the disease-free periodic equilibrium is globally attractive; and if Rv > 1 and R0 > 1 the disease remains persistent. Finally, using the reported monthly mean temperature of Burkina Faso, we perform some numerical simulations to illustrate our theoretical results. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Malaria is a protozoan infection of red blood cells caused in the human body by five species of the genus Plasmodium: P. falciparum, P. vivax, P. ovale, P. malariae and P. knowlesi. The parasites are generally transmitted to the human host through the bite of an infectious female anopheline mosquito. The incidence of malaria is approximately 300–500 million clinical cases, resulting in 1 million deaths each year, where children under 5 are the most chiefly burdened. Of all infectious diseases, malaria continues to be one of the biggest contributors to the global disease burdens in terms of suffering and death and keeps on receiving worldwide attention [1,2]. In order to stop the spread of infectious diseases, Epidemic Dynamics formulates mathematical models, based on the occurrence and progressions of diseases and the surroundings, ∗ Corresponding author. E-mail addresses: [email protected] (B. Traoré), [email protected] (O. Koutou), [email protected] (B. Sangaré).

https://doi.org/10.1016/j.nonrwa.2019.103081 1468-1218/© 2019 Elsevier Ltd. All rights reserved.

2

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

to characterize the infectious agents, to describe the transmission processes, to analyze the origins of the diseases and the factors involved in the transmissions, and to predict the prevalence of the diseases and their patterns. Thus, malaria is one of the first human diseases to be subject to mathematical inquiry [3–5]. Indeed, the mathematical modeling of malaria transmission began in the 20th century (1911) with the first model of Ronald Ross [6]. Ross proposed a mechanistic transmission model including the human host and the mosquito. In the early 1950s, Georges Macdonald [7] improved Ross’s model and introduced an expression for the basic reproduction number. Since then, there has been a great deal of work about using mathematical models to study malaria transmission [8–10]. The standard technique for developing mathematical descriptions of mosquito-plasmodium interactions is to model the system as a set of autonomous ordinary differential equations. This is an immensely powerful approach which has led to many insights into the factors that affect malaria prevalence and control. However, in the mathematical modeling of malaria, two major points are often less investigated by authors. The first point is the stage structure of the vector. While mathematical models for malaria transmission growth abound in literature, the mosquito population has been assumed to be homogeneous, or its stage structure has been simplified in most of the models. In fact, immature mosquitoes are not considered because they do not fly and do not bite humans, so they do not participate in the infection cycle. However, since the abundance and dispersal of adult females are the key factors for the disease transmission, then the dispersal pattern of adult females needs to be studied and modeled. Thus, in order to have more realistic models of mosquitoes, we need to include the immature stages. In particular, the different stages have different responses to the environment and the regulating factors to the population. Moreover, the immature control is a hotly recommended strategy in the fight against malaria, so, it is important to incorporate the immature stages in the malaria models. Several recent works taking into account the stage structure of mosquito have proved that including immature mosquitoes in the model has a profound impact on malaria transmission [11–13]. The second point is the climate effects on the dynamics of vector population. Since mosquito biology and disease ecology are strongly linked to environmental conditions, then the impacts of mosquito-borne diseases increase with warm temperatures and extreme precipitations. There is an evidence to show that climate change and increasing precipitation of the earths temperature will hasten maturation of the parasite in mosquitoes and so increase the biting frequency and create favorable conditions for mosquito breeding. In fact, increasing temperature is known to decrease the length of time spent in each of the immature stages and to decrease the lifespan of adults. In addition, the death rate of the immature stages is strongly linked to temperature. Clearly, temperature exhibits complex and opposing effects on different parts of the mosquito life and transmission cycles. As such, mathematical models are needed which explicitly incorporate the impact of temperature on each life stage if we wish to understand its effect on seasonal abundance [14–17]. Most recently, Wang et al. [18] have formulated a malaria model which includes immature mosquitoes and climate effects. However, their model is limited due to the following aspects: the immature mosquitoes have been grouped in a single compartment, the extrinsic incubation period in mosquitoes and human hosts have been neglected and the oviposition rate of juvenile mosquito is assumed to be linear. Further, in a study, Okuneye et al. [16] developed and analyzed several complex weather-driven mechanistic models that extend the prior studies by incorporating a broader array of biological, ecological and epidemiological factors such as the dynamics of immature mosquitoes, host age-structure, host immunity-boosting due to repeated exposure to malaria infection and the effects of climate on the vector dynamics. Moreover, through sensitivity analysis of basic reproduction ratio, they found some important parameters which influence malaria transmission in the periodic environments. Despite its virtues, the model of Okuneye et al. does not include all the immature stages of mosquito life cycle. Motivated by the works done in [13,16,19], we formulate a mosquito-stage-structured malaria model by taking into account the four distinct metamorphic stages of mosquitoes which is subject to temperature

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

3

Fig. 1. Different stages of biological evolution of mosquito population.

variations. Thus, the model is a non-autonomous system of ordinary differential equations constructed thanks to the coupling of two subsystems. The first subsystem describes the dynamics of mosquito population and the second subsystem describes the dynamics of malaria transmission. In order to get qualitative dynamics, we use the Floquet theory of periodical system [20], the theory of uniform persistence, the limiting system and the theory of chain transitive set [21,22] to analyze the global dynamics of the model. This paper is organized as follows. In Section 2, we formulate the model. In Section 3, we introduce some basic results and then analyze the global behavior of the models. In Section 4, we perform some numerical analysis in order to support our main results. In addition, we examine the impact of some control measures to detect the most optimal for the country. A brief discussion is given in Section 5 to conclude this work and some possible future works are presented. 2. Model formulation 2.1. Formulation of vector population growth dynamics The complete metamorphosis of the mosquito entails goes through four distinct stages of development, namely egg, larva, pupa, and adult mosquito stages (see Fig. 1). The first three stages are aquatic and the last stage is aerial [23,24]. Malaria is transmitted to human through an effective bite from an infectious female mosquito. Female mosquitoes take blood meals to carry out their egg production, and such blood meals are the link between the human and the mosquito. Adult females lay hundreds of eggs per oviposition. The eggs are laid singly directly on the surface of the water. Hatching of eggs into larvae may occur either within a few days or may be delayed for several months depending on the period (temperature) of the year during which the eggs are laid. Anopheles mosquitoes develop through four larval sizes before pupating. Larvae are very small in the first size and increase in size until reaching 5–6 mm by the completion of the fourth size. They feed on organic matters and algae. Moreover, it has been observed that larvae are cannibal because they are able to eat earlier-stage larvae under certain conditions [18,25]. Thus, after the larvae have completed moulting, depending on the temperature, they become pupae. A pupa is a resting, non-feeding development stage. This is the time the pupa changes into an adult. When the development is complete, depending on the temperature the pupal skin splits and the adult mosquito emerges [26,27]. Therefore at each stage, mosquitoes leave the population through natural death rates which depend on the climatic profile. Hence, in the model we divided the mosquito population into two groups: adult mosquitoes denoted by Nv and immature mosquitoes. Immature mosquitoes are

4

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

Fig. 2. Transfer diagram of mosquito population growth: the solid arrows represent the transition from one class to another and the dashed arrow represents the eggs laying of female mosquitoes.

also divided into three compartments: E, L and P which represent the number of eggs, larvae and pupae respectively [28]. From Fig. 2, we obtain the following system: ⎧ ⎪ ⎪ ⎪ ⎨

˙ E(t) = b(t)Nv (t) − (α1 (t) + γ1 (t))E(t), ˙ L(t) = α1 (t)E(t) − (α2 (t) + γ2 (t))L(t),

⎪ P˙ (t) = α2 (t)L(t) − (α3 (t) + γ3 (t))P (t), ⎪ ⎪ ⎩ N˙ v (t) = α3 (t)P (t) − γ4 (t)Nv (t).

(1)

Mosquitoes develop various approaches to avoid predation. The oviposition habitat selection is made to make sure that the mosquito larvae can develop relatively unmolested. These selected habitats are largely devoid of predators, thus offering a relatively secure refuge for larval development [29–32]. Indeed, if there are too many eggs in the oviposition habitat or too few nutrients and water resources, then females laid less eggs or choose another site. It seems reasonable to express this biological phenomenon with a mathematical model which explicitly incorporates the idea of limited carrying capacity resources. This model should take into account the availability of nutrients and the occupation by eggs or larvae of the available breeder sites [11,25]. Thus, the temperature-dependent per capita oviposition rate is given by ( ) E(t) b(t) 1 − Nv (t), KE the number of eggs that survive and become larvae is given by ) ( L(t) E(t), α1 (t) 1 − KL and the number of larvae that survive and become pupae is given by ( ) P (t) α2 (t) 1 − L(t). KP Finally, the number of mosquitoes that survive from recruitment into one class, to the next is given by the following sequence of ordinary differential equations: ( ) ⎧ E(t) ⎪ ˙ ⎪ E(t) = b(t) 1 − Nv (t) − (α1 (t) + γ1 (t))E(t), ⎪ ⎪ KE ⎪ ⎪ ( ) ⎪ ⎪ ⎪ L(t) ⎨ L(t) ˙ = α1 (t) 1 − E(t) − (α2 (t) + γ2 (t))L(t), KL (2) ( ) ⎪ ⎪ ⎪ P (t) ⎪ ⎪ L(t) − (α3 (t) + γ3 (t))P (t), ⎪ P˙ (t) = α2 (t) 1 − K ⎪ ⎪ P ⎪ ⎩ ˙ Nv (t) = α3 (t)P (t) − γ4 (t)Nv (t).

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

5

Table 1 Biological description of model (2) parameters. Parameters

Biological descriptions

KE KL KP α1 (t) α2 (t) α3 (t) b(t) γ1 (t) γ2 (t) γ3 (t) γ4 (t)

Carrying capacity of eggs Carrying capacity of larvae Carrying capacity of pupae Temperature-dependent transfer rate from eggs to larvae Temperature-dependent the transfer rate from larvae to pupae Temperature-dependent transfer rate from pupae to adult Temperature-dependent eggs laying rate per female Temperature-dependent death rate of eggs Temperature-dependent death rate of larvae Temperature-dependent death rate of pupae Temperature-dependent death rate of adult mosquitoes

Table 2 Parameters of malaria transmission model. Parameters

Biological descriptions

β(t)

Temperature-dependent biting rate of mosquitoes Constant recruitment for humans Human natural death rate Transfer rate of humans from exposed class to infectious class Disease-induced death rate for humans Recovery rate of humans Per capita rate of loss of immunity for humans Transfer rate of mosquitoes from exposed class to infectious class Probability of transmission of infection from Ih to Sh Probability of transmission of infection from Ih to Sv Probability of transmission of infection from Rh to Sv

Λh

dh α κ rh γ νv cvh chv c ¯hv

2.2. Formulation of malaria transmission model In the presence of the disease, the human population is divided into four epidemiological categories representing the state variables: the susceptible class Sh , the exposed class Eh , the infectious class Ih and the immune (asymptomatic, but slightly infectious humans) class Rh . Indeed, individuals in the class Rh are immune in the sense that they do not suffer from serious illness and do not contract clinical malaria, but they still have low levels of Plasmodium in their blood stream and can transmit the infection to susceptible mosquitoes [17,33,34]. Similarly, adult mosquitoes are divided into three compartments Sv , Ev and Iv that represent the susceptible mosquitoes, the exposed mosquitoes and the infectious mosquitoes, respectively. When males and females finish to mate, females will take blood meals for their eggs development. Thus, if a susceptible mosquito bites an infectious human, the mosquito becomes infected and it moves into the class Ev . Some time after, it becomes infectious and enters into the class Iv . Hence, this infectious mosquito can infect a susceptible human through another bite. Some time after infection, the human becomes infectious. If he is continually exposed to mosquito bites, he becomes immune. However, he can lose this immunity and become susceptible if the exposure to infection is stopped. Humans leave the total population through natural death rate and malaria death rate; and mosquitoes leave the population through natural death. At any time t ≥ 0, the total size of humans and adults mosquitoes is given respectively by Nh (t) = Sh (t) + Eh (t) + Ih (t) + Rh (t),

(3)

Nv (t) = Sv (t) + Ev (t) + Iv (t).

(4)

The parameters of the model are given in Table 1 (see Table 2).

6

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

It is assumed throughout this paper that: (H1): all vector population measures refer to densities of female mosquitoes, (H2): the mosquitoes bite only humans, (H3): there is no direct transmission (blood transfusion or from mother to baby) of malaria, (H4): all the new recruits are susceptible, (H5): there is no immigration of infectious humans. Using the standard incidence, we define respectively the infection incidence from mosquitoes to humans, kh (t) and from humans to mosquitoes, kv (t): kh (t) = cvh β(t)

Iv (t) Nh (t)

and kv (t) = chv β(t)

Ih (t) Rh (t) + c¯hv β(t) . Nh (t) Nh (t)

where β(t) is the average number of bites per mosquito per unit time t. According to the above assumptions, our model of malaria transmission is represented by Fig. 3. 2.3. Main model Following the above assumptions, the number of individuals which survive from recruitment into one class, to the next (see Fig. 3) is defined by the following sequence of ordinary differential equations: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

( ) E(t) ˙ Nv (t) − (α1 (t) + γ1 (t))E(t), E(t) = b(t) 1 − KE ( ) L(t) ˙ L(t) = α1 (t) 1 − E(t) − (α2 (t) + γ2 (t))L(t), KL ( ) P (t) P˙ (t) = α2 (t) 1 − L(t) − (α3 (t) + γ3 (t))P (t), KP S˙ v (t) = α3 (t)P (t) − (γ4 (t) + kv (t))Sv (t), E˙ v (t) = kv (t)Sv (t) − (νv + γ4 (t))Ev (t), I˙v (t) = νv Ev (t) − γ4 (t)Iv (t), S˙ h (t) = Λh + γRh (t) − (dh + kh (t))Sh (t), E˙ h (t) = kh (t)Sh (t) − (dh + α)Eh (t), I˙h (t) = αEh (t) − (dh + κ + rh )Ih (t), R˙ h (t) = rh Ih (t) − (dh + γ)Rh (t).

Since, Sv (t) = Nv (t) − Ev (t) − Iv (t), ∀t ≥ 0,

(5)

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

7

Fig. 3. Transfer diagram of malaria transmission to humans: the black dashed arrows indicate the direction of the infection and the solid arrows represent the transition from one class to another.

then, system (5) reads as follows: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

( ) E(t) ˙ E(t) = b(t) 1 − Nv (t) − (α1 (t) + γ1 (t))E(t), KE ( ) L(t) ˙ L(t) = α1 (t) 1 − E(t) − (α2 (t) + γ2 (t))L(t), KL ( ) P (t) P˙ (t) = α2 (t) 1 − L(t) − (α3 (t) + γ3 (t))P (t), KP N˙ v (t) = α3 (t)P (t) − γ4 (t)Nv (t), E˙ v (t) = kv (t)Nv (t) − kv (t)Iv (t) − (νv + γ4 + kv (t))Ev (t),

(6)

I˙v (t) = νv Ev (t) − γ4 (t)Iv (t), S˙ h (t) = Λh + γRh (t) − (dh + kh (t))Sh (t), E˙ h (t) = kh (t)Sh (t) − (dh + α)Eh (t), I˙h (t) = αEh (t) − (dh + κ + rh )Ih (t), R˙ h (t) = rh Ih (t) − (dh + γ)Rh (t),

with initial conditions: E(0) > 0, L(0) > 0, P (0) > 0, Sh (0) > 0, Eh (0) > 0, Ih (0) > 0, Rh (0) > 0, Nv (0) > 0, Ev (0) > 0, Iv (0) > 0. At any time t ≥ 0, we have: N˙ h (t) = Λh − dh Nh (t) − κIh (t), N˙ v (t) = α3 (t)P (t) − γ4 (t)Nv (t).

(7) (8)

8

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

3. Preliminaries and main results In this section, we focus on the global behavior of model (6). Before starting our study, we enumerate some hypotheses and recall some fundamental results which will be used in our further analysis. (H6): Constant parameters of the model are positive, (H7): γ1 (t), γ2 (t), γ3 (t), γ4 (t), b(t) and β(t) are positive, continuous and periodic functions with the same period ω. 3.1. Preliminaries Let (Rn+ , Rn+ ) be the standard ordered n-dimensional Euclidean space with a norm ∥.∥. For u, v ∈ Rn , we denote u ≥ v, if u − v ∈ Rn+ , u > v, if u − v ∈ Rn+ \ {0}, u ≫ v, if u − v ∈ int(Rn+ ). For a continuous positive ω-periodic function W, we define ˆ = sup W(t) and W ¯ = inf W(t). W t∈[0,ω]

t∈[0,ω]

Let us consider the following linear ordinary differential system: z(t) ˙ = N (t)z(t),

(9)

where N (t) is a continuous, cooperative, irreducible and ω-periodic n × n matrix function. Let ΦN (t) be ( ) the fundamental solution matrix of (9) and ρ ΦN (ω) be the spectral radius of the matrix ΦN (ω). By the ( ) Perron–Frobenius theorem (see [35], Theorem A.3), ρ ΦN (ω) is the principal eigenvalue of ΦN (ω) which is associated to the eigenvector v ∗ ≫ 0. The following result is useful for our subsequent comparison arguments. Lemma 3.1 ([36]). Let r = v(t)ert is a solution of (9).

1 ω

ln ρ(ΦN (ω)), then there exists a positive ω-periodic function v(t) such that

Definition 3.1 ([1]). Let X be a metric space with metric d and g : X → X a continuous map. For any x ∈ X , we denote g n (x) = g(g n−1 (x)) for any integer n > 1 and g 1 (x) = g(x). g is said to be compact in X if for any bounded set H ∈ X , the set g(H) = {g(x) : x ∈ H} is precompact in X . g is said to be point dissipative if there is a bounded set B0 ∈ X such that B0 attracts each point in X . ( ) lim d g n (x), B0 = 0 for any x0 ∈ X . n→∞

The positive semiorbit through x0 is defined by γ + (x0 ) = {xn = g n (x0 ), n = 1, 2, . . .}

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

9

and the negative semiorbit through x0 is defined as a sequence γ − (x0 ) = {xk } satisfying g(xk−1 ) = xk for integers k ≤ 0. The omega limit set of γ + (x0 ) is defined by ω ′ (x0 ) = {c ∈ X : there is a sequence nk → ∞ such that

lim xnk = c} k→∞

and the alpha limit set of γ − (x0 ) is defined by α′ (x0 ) = {c ∈ X : there is a sequence nk → −∞ such that

lim xnk = c}. k→∞

A nonempty set C ⊂ X is said to be invariant if g(C) ⊆ C. A nonempty invariant set M of X is called to be isolated in X if it is the maximal invariant set in the neighborhood of itself. For a nonempty set M of X , the set W s (M ) := {x ∈ X : lim d(g n (x), M ) = 0} n→∞

is called the stable set of M . Let A1 and A2 be two isolated invariant sets. The set A1 is said to be chained to set A2 , written A1 → A2 if there exists a full orbit through some x ∈ / A1 ∪ A2 such that ω ′ (x) ⊂ A2 and α′ (x) ⊂ A1 . A finite sequence {M1 , . . . , Mk } of isolated invariant sets is called a chain if M1 → M2 → M3 → ... → Mk , and if Mk = M1 then the chain is called a cycle. Let X0 be a nonempty open set of X . We denote ∂X0 = X \ X0 , and M∂ := {x ∈ ∂X0 : g n (x) ∈ ∂X0 , ∀ n ≥ 0}. Lemma 3.2 ([22]). Let g : X → X a continuous map. Assume that the following conditions hold: (i) g is compact and point dissipative, and g(∂X0 ) ⊆ ∂X0 , (ii) there exists a finite sequence M = {M1 , . . . , Mk } of compact and isolated invariant sets such that (a) (b) (c) (d)

Mi ∩ Mj = ∅ for any i, j = 1, 2, . . . , k and i ̸= j; ⋃ ⋃k ′ x∈M∂ ω (x) ⊂ i=1 Mi ; no subset of M forms a cycle in ∂X0 ; W s (Mi ) ∩ X0 = ∅ for each 1 ≤ i ≤ k.

Then, g is uniformly persistent with respect to (X0 , ∂X0 ), that is, there exists a positive constant η such that lim inf n→∞ d(g n (x), ∂X0 ) ≥ η. Theorem 3.1 ([37,38]). Let G : Rn −→ R be a differential function and let f (x) be a vector fields. Let λ ∈ R { } { } and let Γ = x ∈ Rn : G(x) ≤ λ be such that ∇G(x) ̸= 0 for all x ∈ x ∈ Rn : G(x) = λ = G−1 (λ). If ⟨ ⟩ f (x), ∇G(x) ≤ 0 for all x ∈ G−1 (λ), then the set Γ is positively invariant. 3.2. Dynamics analysis of vector population model Let us consider the following system: ) ( ⎧ E(t) ⎪ ˙ ⎪ E(t) = b(t) 1 − Nv (t) − (α1 (t) + γ1 (t))E(t), ⎪ ⎪ KE ⎪ ⎪ ( ) ⎪ ⎪ ⎪ L(t) ⎨ L(t) ˙ = α1 (t) 1 − E(t) − (α2 (t) + γ2 (t))L(t), KL ( ) ⎪ ⎪ ⎪ P (t) ⎪ ⎪ P˙ (t) = α2 (t) 1 − L(t) − (α3 (t) + γ3 (t))P (t), ⎪ ⎪ KP ⎪ ⎪ ⎩ ˙ Nv (t) = α3 (t)P (t) − γ4 (t)Nv (t).

(10)

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

10

System (10) can be written as follows: {

( ) Z˙ 1 (t) = f1 t, Z1 (t) , Z1 (0) > 0,

(11)

where Z1 (t) = (E(t), L(t), P (t), Nv (t))T . 3.2.1. Positivity and boundedness of solutions Lemma 3.3. For any positive initial condition φ1 = (E(0), L(0), P (0), Nv (0)), system (10) admits a unique positive solution u1 (t, φ1 ) for all t ≥ 0. Moreover, the closed set } { α ˆ 3 KP ∆1 = (E, L, P, Nv ) ∈ R4+ : E ≤ KE , L ≤ KL , P ≤ KP , Nv ≤ γ¯4 is positively invariant and attracting under the flow described by system (10). Proof . For any positive initial condition φ1 , the function f1 (t, φ1 ) is continuous in (t, φ1 ) and Lipschitzian in φ1 . So, thanks to Theorem 2.2.1 and 2.2.3 of Hale and Verduyn Lunel [21], system (10) has a unique solution u1 (t, φ1 ) on its maximal interval [0, tmax ) of existence. In addition, for all t ≥ 0, let { } e(t) = min E(t), L(t), P (t), Nv (t) . Suppose that there exists t1 > 0 such that e(t1 ) ∈ / R∗+ and e(t) > 0 for all t ∈ [0, t1 ). If e(t) = Nv (t) then, from the fourth equation of system (10), we have N˙ v (t) > −γ4 (t)Nv (t). It then follows that

( ∫ Nv (t1 ) > Nv (0) exp −

t1

) γ4 (t)dt > 0,

0

which leads to a contradiction. If e(t) = E(t) then, from the first equation of system (10), we have ( ) ˙ E(t) > − α1 (t) + γ1 (t) E(t). It then follows that

( ∫ E(t1 ) > E(0) exp −

t1 (

) α1 (t) + γ1 (t) dt > 0, )

0

which leads to a contradiction. Similar contradictions can be obtained if e(t) = L(t) and e(t) = P (t). Hence, the solution u1 (t, φ1 ) is positive. Moreover, let { } Q1 = (E, L, P, Nv ) ∈ R4+ : E ≤ KE , { } Q3 = (E, L, P, Nv ) ∈ R4+ : P ≤ KP ,

{ } Q2 = (E, L, P, Nv ) ∈ R4+ : L ≤ KL , { } α ˆ3 KP . Q4 = (E, L, P, Nv ) ∈ R4+ : Nv ≤ γ¯4

Let us show that Q1 , Q2 , Q3 and Q4 are positively invariant. Let G1 (E, L, P, Nv ) = E, G2 (E, L, P, Nv ) = L, G3 (E, L, P, Nv ) = P and G4 (E, L, P, Nv ) = Nv .

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

11

It is clear that Gi (E, L, P, Nv ) is differentiable for i = 1, 2, 3, 4. It then follows that: ∇G1 (E, L, P, Nv ) = (1, 0, 0, 0),

∇G2 (E, L, P, Nv ) = (0, 1, 0, 0),

∇G3 (E, L, P, Nv ) = (0, 0, 1, 0),

∇G4 (E, L, P, Nv ) = (0, 0, 0, 1).

Moreover, let { } Γ1 = (E, L, P, Nv ) ∈ R4+ : E = KE , { } Γ3 = (E, L, P, Nv ) ∈ R4+ : P = KP ,

{ } Γ2 = (E, L, P, Nv ) ∈ R4+ : L = KL , } { α ˆ3 KP . Γ4 = (E, L, P, Nv ) ∈ R4+ : Nv = γ¯4

Thus, if: • • • •

x ∈ Γ1 , x ∈ Γ2 , x ∈ Γ3 , x ∈ Γ4 ,

then then then then

⟨ ⟩ ( ) f1 (t, x), ∇G1 (x) = − α1 (t) + γ1 (t) < 0. ⟨ ⟩ ( ) f1 (t, x), ∇G2 (x) = − α2 (t) + γ2 (t) < 0. ⟨ ⟩ ( ) f1 (t, x), ∇G3 (x) = − α3 (t) + γ3 (t) < 0. ⟨ ⟩ f1 (t, x), ∇G4 (x) = α3 (t)P (t) − γ4 (t) αγˆ¯43 KP ≤ α ˆ 3 (P (t) − KP ) ≤ 0.

So, thanks to Theorem 3.1, we conclude that the compact set ∆1 is positively invariant and then the solution u1 (t, φ1 ) is bounded. It then follows that tmax = +∞. □ 3.2.2. Vector reproduction ratio We introduce the threshold dynamics for system (10) according to the theory developed by Wang and Zhao [39] which is a generalization of Van den Driessche and James Watmough method, [40]. Linearizing system (10) at the unique mosquito-free equilibrium H 0 = (0, 0, 0, 0), we obtain the following system: ⎧ ˙ E(t) = b(t)Nv (t) − (α1 (t) + γ1 (t))E(t), ⎪ ⎪ ⎪ ⎨ L(t) ˙ = α1 (t)E(t) − (α2 (t) + γ2 (t))L(t),

(12)

⎪ P˙ (t) = α2 (t)L(t) − (α3 (t) + γ3 (t))P (t), ⎪ ⎪ ⎩ N˙ v (t) = α3 (t)P (t) − γ4 (t)Nv (t), which can be written as follows: U˙ 1 (t) =

( ) F1 (t) − V1 (t) U1 (t)

with U1 (t) = (E(t), L(t), P (t), Nv (t))T and ⎛ ⎜ ⎜ V1 (t) = ⎜ ⎝

α1 (t) + γ1 (t)

0

0

0

⎞

−α1 (t)

α2 (t) + γ2 (t)

0

0

0

−α2 (t)

α3 (t) + γ3 (t)

0

⎟ ⎟ ⎟, ⎠

0

0

−α3 (t) ⎞

γ4 (t)

⎛

0 0 0 ⎜ 0 0 0 F1 (t) = ⎜ ⎝ 0 0 0 0 0 0

b(t) 0 ⎟ ⎟. 0 ⎠ 0

For all t ≥ s, let Y˜ (t, s) be the evolution operator of the linear periodic system y˙ = −V1 (t)y.

12

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

For each s ∈ R, the 4 × 4 matrix Y˜ (t, s) satisfies Y˙ (t, s) = −V1 (t)Y (t, s), ∀t ≥ s, Y (s, s) = I,

(13)

where I is the 4 × 4 identity matrix. Let Cω be the ordered Banach space of all ω-periodic functions from R to R4 which is equipped with the maximum norm ∥.∥ and the positive cone { } C+ ω := ϕ ∈ Cω : ϕ(t) ≥ 0, ∀t ∈ R . Suppose ϕ(s) ∈ Cω is the initial distribution of juvenile and adults mosquitoes. Then, F1 (s)ϕ(s) ∈ Cω is the distribution of new juvenile mosquitoes in the breeding habits produced by the adult ones which were introduced at time s. Hence, for any t ≥ s, Y (t, s)F1 (s)ϕ(s) is the distribution of those mosquitoes which were newly born into the juvenile mosquito compartment at time s and remain alive. Thus, ∫ t ∫ ∞ ˜ = ψ(t) Y (t, s)F1 (s)ϕ(s)ds = Y (t, t − a)F1 (t − a)ϕ(t − a)da, a ∈ [0, ∞) −∞

0

is the distribution of new eggs at time t, hatched by all female adult mosquitoes ϕ(s) introduced at the previous time s. Let L : Cω −→ Cω be the linear operator defined by ∫ ∞ (Lϕ)(t) = Y (t, t − a)F1 (t − a)ϕ(t − a)da, ∀t ∈ R, ϕ ∈ Cω . 0

Then, the vector reproduction ratio is Rv := ρ(L), the spectral radius of L. 3.2.3. Global stability of mosquito-free equilibrium Lemma 3.4 ([39]). (i) Rv = 1 if and only if ρ(ΦF1 −V1 (ω)) = 1. (ii) Rv < 1 if and only if ρ(ΦF1 −V1 (ω)) < 1. (iii) Rv > 1 if and only if ρ(ΦF1 −V1 (ω)) > 1. Hence, if Rv < 1 then the mosquito-free equilibrium is locally asymptotically stable. Theorem 3.2.

The mosquito-free equilibrium H 0 = (0, 0, 0, 0) is globally attractive if Rv < 1.

Proof . For all t ≥ 0, since E(t) ≤ KE , L(t) ≤ KL and P (t) ≤ KP , we have: 1−

L(t) P (t) E(t) ≤ 1, 1 − ≤ 1 and 1 − ≤ 1. KE KL KP

Thus, system (10) becomes: ⎧ ˙ E(t) ≤ b(t)Nv (t) − (α1 (t) + γ1 (t))E(t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ˙ ⎪ ≤ α1 (t)E(t) − (α2 (t) + γ2 (t))L(t), ⎨ L(t) ⎪ ⎪ P˙ (t) ≤ α2 (t)L(t) − (α3 (t) + γ3 (t))P (t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ˙ Nv (t) = α3 (t)P (t) − γ4 (t)Nv (t).

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

13

Let us consider the following auxiliary system: ⎧ ˙ ˜ = b(t)N ˜v (t) − (α1 (t) + γ1 (t))E(t), ˜ E(t) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ˜˙ ˜ − (α2 (t) + γ2 (t))L(t), ˜ ⎪ = α1 (t)E(t) ⎨ L(t) ⎪ ⎪ ˜ − (α3 (t) + γ3 (t))P˜ (t), ⎪ P˜˙ (t) = α2 (t)L(t) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ˜˙ ˜v (t). Nv (t) = α3 (t)P˜ (t) − γ4 (t)N For the convenience, we will rewrite it as follows ( ) ˜˙ = F1 (t) − V1 (t) J(t), ˜ J(t)

(14)

with ˜ = (E(t), ˜ ˜ ˜v (t))T . J(t) L(t), P˜ (t), N If Rv < 1, then from Lemma 3.4, ρ(ΦF1 −V1 (ω)) < 1. By Lemma 3.1, there exists a positive ω-periodic ˜ = v(t)ert with r = 1 ln ρ(ΦF −V (ω)). Since the function v(t) is bounded and function v(t) such that J(t) 1 1 ω ˜ → 0 as t → ∞. By applying the comparison theorem [41] ρ(ΦF1 −V1 (ω)) < 1 then, we have r < 0 and J(t) on system (14), we get lim (E(t), L(t), P (t), Nv (t)) = (0, 0, 0, 0). t→+∞

It then follows that the mosquito-free equilibrium H 0 is globally attractive which completes the proof.

□

3.2.4. Persistence of mosquitoes Let us define the following sets: X := R4+ , X0 := {(E, L, P, Nv ) ∈ X : E > 0, L > 0, P > 0, Nv > 0}, ∂X0 := {(E, L, P, Nv ) ∈ X : ELP Nv = 0}. Let Φ(t) be the periodic semiflow generated by the periodic system (10) and S1 : X −→ X the Poincar´e map associated with system (10), namely: S1 (φ1 ) = Φ(ω)φ1 = u1 (ω, φ1 ), ∀φ1 ∈ X. S1m (φ1 ) = Φ(mω)φ1 = u1 (mω, φ1 ), ∀m ≥ 0. { } It then follows that the study of the dynamics (10) reduces to that of the discrete dynamical system S1m on X. Otherwise, note that for all φ1 ∈ X0 , Φ(t)(φ1 ) = u1 (t, φ1 ) ∈ X0 . Thus, Φ(t)(X0 ) ⊂ X0 , for all t ≥ 0. So, X0 and ∂X0 are positively invariant. Therefore, Lemma 3.1 implies that discrete-time system S1 is point dissipative. Lemma 3.5.

If Rv > 1, there exists δ > 0 such that when   φ1 − H 0  ≤ δ, for any φ1 ∈ X0 ,

one has ( ) lim sup d S1k (φ1 ), H 0 ≥ δ.

k→∞

14

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

Proof . Suppose by contradiction that ( ) lim sup d S1k (φ1 ), H 0 < δ for some

φ1 ∈ X0 .

(15)

k→∞

( ) Then, there exists an integer k2 such that for all k ≥ k2 , d S1k (φ1 ), H 0 < δ. By the continuity of the solution u1 (t, φ1 ), we have  ( k  ) u1 t, S1 (φ1 ) − u1 (t, H 0 ) ≤ α, ∀t ≥ 0 and α > 0. [ ] [ ] Let t = kω + t1 , where t1 ∈ [0, ω] and k = ωt . ωt is the greatest integer less than or egal to ωt . If   φ1 − H 0  ≤ δ, then we get  (   (  ) ) u1 t, φ1 − u1 (t, H 0 ) = u1 t1 + kω, φ1 − u1 (t1 + kω, H 0 )  (  ) = u1 t1 , S1k φ1 − u1 (t1 , H 0 ) ≤ α, for any t ≥ 0. It then follows that for all t ≥ 0, 0 ≤ E(t) ≤ α, 0 ≤ L(t) ≤ α and 0 ≤ P (t) ≤ α. } { Thus, there exists α∗ = max

1−

α α α KE , KL , KP

such that

E(t) L(t) P (t) ≥ 1 − α∗ , 1 − ≥ 1 − α∗ and 1 − ≥ 1 − α∗ , for all t ≥ 0. KE KL KP

Furthermore, system (10) can ⎧ ⎪ ⎪ ⎪ ⎨

be rewritten as follows: ˙ E(t) ≥ b(t)(1 − α∗ )Nv (t) − (α1 (t) + γ1 (t))E(t), ˙ L(t) ≥ α1 (t)(1 − α∗ )E(t) − (α2 (t) + γ2 (t))L(t),

(16)

⎪ P˙ (t) ≥ α2 (t)(1 − α∗ )L(t) − (α3 (t) + γ3 (t))P (t), ⎪ ⎪ ⎩ N˙ v (t) = α3 (t)P (t) − γ4 (t)Nv (t). Let us consider the following auxiliary linear system ˙ R(t) = Mα∗ (t)R(t),

(17)

where R(t) = (E(t), L(t), P (t), Nv (t))T and

⎛ ⎜ ⎜ Mα∗ (t) = ⎜ ⎝

( ) − α1 (t) + γ1 (t) α1 (t)(1 − α∗ )

0 ( ) − α2 + γ2 (t)

0

0

α2 (t)(1 − α∗ )

− α3 (t) + γ3 (t)

0

0

α3 (t)

0 (

)

b(t)(1 − α∗ )

⎞

0

⎟ ⎟ ⎟. ⎠

0 −γ4 (t)

Since limα∗ →0+ ΦMα∗ (ω) = ΦF1 −V1 (ω) and by the continuity of the spectral radius, we have ( ) limα∗ →0+ ρ ΦMα∗ (ω) = ρ(ΦF1 −V1 (ω)). From Lemma 3.4 if Rv > 1, then ρ(ΦF1 −V1 (ω)) > 1 and ( ) ( ) limα∗ →0+ ρ ΦMα∗ (ω) > 1. It then follows that there exists α1 > 0 such that ρ ΦMα∗ (ω) > 1, for all α ∈ [0, α1 ]. From Lemma 3.1, there exists a positive ω-periodic function v(t) such that R(t) = v(t)ert with ( ) ( ) r = ω1 ln ρ ΦMα∗ (ω) . Since ρ ΦMα∗ (ω) > 1 and function v(t) is bounded, then r > 0 and R(t) → ∞ as t → ∞. By applying the comparison theorem on system (17), we get lim |(E(t), L(t), P (t), Nv (t))| = ∞,

t→∞

which is in contradiction with (15).

□

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

15

Theorem 3.3. If Rv > 1, there exists η > 0 such that any solution (E(t), L(t), P (t), Nv (t)) with initial condition φ1 ∈ X0 satisfies lim inf E(t) ≥ η, lim inf L(t) ≥ η, lim inf P (t) ≥ η, lim inf Nv (t) ≥ η t→∞

t→∞

t→∞

t→∞

and system (10) has at least one positive periodic solution. Proof . Denote { } M∂ = φ1 ∈ ∂X0 : S1k (φ1 ) ∈ ∂X0 , k ≥ 0 . At first, we prove that M∂ = {(0, 0, 0, 0)}. It is clear that {(0, 0, 0, 0)} ⊂ M∂ . So, we only need to show that M∂ ⊂ {(0, 0, 0, 0)}, which means that for any initial condition φ1 ∈ ∂X0 , E(kω) = 0, or L(kω) = 0 or P (kω) = 0 or Nv (kω) = 0, ∀k ≥ 0. Let φ1 ∈ ∂X0 . Suppose by contradiction that there exists an integer k1 ≥ 0 such that (E(k1 T ), L(k1 ω), P (k1 ω), Nv (k1 ω))T > 0. Using the constant variation formula, we derive that ( ∫ t )[ (∫ s ) ] ∫ t E(t) = exp − h1 (s)ds E(0) + b(s)Nv (s) exp h1 (τ )dτ ds , (18) 0 0 0 ( ∫ t )[ (∫ s ) ] ∫ t L(t) = exp − h2 (s)ds L(0) + α1 (s)E(s) exp h2 (τ )dτ ds , (19) 0 0 ( ∫0 t )[ (∫ ) ] ∫ t s P (t) = exp − h3 (s)ds P (0) + α2 (s)L(s) exp h3 (τ )dτ ds , (20) 0 0 ( ∫0 t )[ ( ) ] ∫ t ∫ s Nv (t) = exp − γ4 (s)ds Nv (0) + α3 (s)P (s) exp − γ4 (τ )dτ ds , (21) 0

0

0

where b(t) Nv (t) + α1 (t) + γ1 (t), KE α1 (t) E(t) + α2 (t) + γ2 (t), h2 (t) = KL α2 (t) h3 (t) = L(t) + α3 (t) + γ3 (t). KP Thus, taking k1 ω as the initial time in (18)–(21), we obtain that E(t) > 0, L(t) > 0, P (t) > 0 and Nv (t) > 0, for all t > k1 ω, which contradicts the fact that ∂X0 is positively invariant so M∂ ⊂ {(0, 0, 0, 0)}. Hence, it then follows that M∂ = {(0, 0, 0, 0)}. Equality M∂ = {(0, 0, 0, 0)} implies that H 0 is a fixed point of S1 and acyclic in M∂ . Thus, every solution in M∂ approaches to H 0 . Moreover, Lemma 3.5 implies that H 0 is an isolate invariant set in X and W s (H 0 )∩X0 = ∅. Therefore, by Theorem 3.1.1 in [22], we finally obtain that S1 is uniformly persistent with respect to (X0 , ∂X0 ). So, the periodic semiflow Φ(t) is also uniformly persistent. It then follows that there exists η > 0 such that any solution (E(t), L(t), P (t), Nv (t)) with initial condition (E(0), L(0), P (0), Nv (0)) ∈ X0 satisfies lim inf E(t) ≥ η, lim inf L(t) ≥ η, lim inf P (t) ≥ η, lim inf Nv (t) ≥ η. h1 (t) =

t→∞

t→∞

t→∞

t→∞

Besides, thanks to Theorem 1.3.6 in [22], system (10) has at least one periodic solution u∗1 (t, φ∗1 ) with φ∗1 = (E ∗ (0), L∗ (0), P ∗ (0), Nv∗ (0)) ∈ X0 . Since φ∗1 ∈ X0 , then from Eqs. (18)–(21) we obtain that E ∗ (t) > 0, L∗ (t) > 0, P ∗ (t) > 0 and Nv∗ (t) > 0 for all t ≥ 0. Thus, the ω-periodic solution u∗1 (t, φ∗1 ) is positive. □ Lemma 3.6 ([42]). If Rv > 1, then solution Nv∗ (t) of Eq. (8) with initial value Nv∗ (0) is bounded and globally uniformly attractive on [0, ∞).

16

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

3.3. Dynamics analysis of malaria transmission model The previous results indicate that the mosquito population will die out if Rv < 1, and if Rv > 1, the mosquito population persists at the ω-periodic steady state (E ∗ (t), L∗ (t), P ∗ (t), Nv∗ (t)). Consequently, when the parasite is introduced in human population or mosquito population, the dynamics of malaria transmission model due to the interaction between humans and mosquitoes is given by the following limiting system: ⎧ S˙ h (t) = Λh + γRh (t) − (dh + kh (t))Sh (t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ E˙ h (t) = kh (t)Sh (t) − (dh + α)Eh (t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ I˙h (t) = αEh (t) − (dh + κ + rh )Ih (t), (22) ⎪ ⎪ ˙ h (t) = rh Ih (t) − (dh + γ)Rh (t), R ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ E˙ v (t) = kv (t)Nv∗ (t) − kv (t)Iv (t) − (νv + γ4 (t) + kv (t))Ev (t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ˙ Iv (t) = νv Ev (t) − γ4 (t)Iv (t). System (22) can be written as follows: {

( ) Z˙ 2 (t) = f2 t, Z2 (t) , Z2 (0) > 0,

(23)

where Z2 (t) = (Sh (t), Eh (t), Ih (t), Rh (t), Ev (t), Iv (t))T . 3.3.1. Positivity and boundedness of solutions Lemma 3.7. For any positive initial condition φ2 = (Sh (0), Eh (0), Ih (0), Rh (0), Ev (0), Iv (0)), system (22) admits a unique positive solution u2 (t, φ2 ) for all t ≥ 0. Moreover, the closed set } { α ˆ3 Λh 6 , Ev + Iv ≤ KP . ∆2 = (Sh , Eh , Ih , Rh , Ev , Iv ) ∈ R+ : Nh ≤ dh γ¯4 is positively invariant and attracting under the flow described by system (22). Proof . For any positive initial condition φ2 , the function f2 (t, φ2 ) is continuous in (t, φ2 ) and Lipschitzian in φ2 . So, thanks to Theorems 2.2.1 and 2.2.3 of Hale and Verduyn Lunel [21], system (22) has a unique solution u(t, φ2 ) on its maximal interval [0, tmax ) of existence. In addition, for all t ≥ 0, let { } e′ (t) = min Sh (t), Eh (t), Ih (t), Rh (t), Ev (t), Iv (t) . Suppose that there exists t1 > 0 such that e′ (t1 ) ∈ / R∗+ and e′ (t) > 0 for all t ∈ [0, t1 ). ′ If e (t) = Sh (t) then, from the second equation of system (22), we have E˙ h (t) > −(dh + α)Eh (t). It then follows that [ ] Eh (t1 ) > Eh (0) exp −(dh + α)t1 > 0, which leads to a contradiction.

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

17

If e′ (t) = Eh (t) then, from the third equation of system (22), we have I˙h (t) > −(dh + κ + rh )Eh (t). It then follows that [ ] Ih (t1 ) > Ih (0) exp −(dh + κ + rh )t1 > 0, which leads to a contradiction. If e′ (t) = Ih (t) then, from the fourth equation of system (22), we have R˙ h (t) > −(dh + γ)Rh (t). It then follows that [ ] Rh (t1 ) > Rh (0) exp −(dh + γ)t1 > 0, which leads to a contradiction. Similar contradictions can be obtained if e′ (t) = Iv (t) and e′ (t) = Ev (t). Hence, the solution u2 (t, φ2 ) is positive. Moreover, let } { Λh 6 , Q5 = (Sh , Eh , Ih , Rh , Ev , Iv ) ∈ R : Sh + Eh + Ih + Rh ≤ dh { } α ˆ3 6 Q6 = (Sh , Eh , Ih , Rh , Ev , Iv ) ∈ R : Ev + Iv ≤ KP γ¯4 and G5 (Sh , Eh , Ih , Rh , Ev , Iv ) = Sh + Eh + Ih + Rh ,

G6 (Sh , Eh , Ih , Rh , Ev , Iv ) = Ev + Iv .

It is clear that G5 and G6 are differentiable and we have: ∇G5 (Sh , Eh , Ih , Rh , Ev , Iv ) = (1, 1, 1, 1, 0, 0) and ∇G6 (Sh , Eh , Ih , Rh , Ev , Iv ) = (0, 0, 0, 0, 1, 1). Moreover, let { Γ5 =

Λh (Sh , Eh , Ih , Rh , Ev , Iv ) ∈ R6 : Sh + Eh + Ih + Rh = dh

and

}

} α ˆ3 KP . (Sh , Eh , Ih , Rh , Ev , Iv ) ∈ R : Ev + Iv = γ¯4

{ Γ6 =

6

We have: ⟨ ⟩ f2 (t, x), ∇G5 (x) = Λh − dh Nh − γ3 Ih ⟨

≤ Λh − dh Nh . ⟩ f2 (t, x), ∇G6 (x) = kv (t)(Nv∗ (t) − Ev − Iv ) − dv(t)(Ev + Iv ) ( ) α ˆ3 ≤ kv (t) KP − (Ev + Iv ) − γ4 (t)(Ev + Iv ). γ¯4

Thus, if: ⟨ ⟩ • x ∈ Γ5 then, f2 (t, x), ∇G5 (x) ≤ 0, ⟨ ⟩ • x ∈ Γ6 then, f2 (t, x), ∇G6 (x) ≤ 0. So, from Theorem 3.1, we conclude that the compact set ∆2 is positively invariant and then all the solutions are nonnegative and bounded. □

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

18

3.3.2. Basic reproduction ratio Using the same method as mentioned above, we derive the threshold dynamics associated to system (22). Hence, by linearizing system (22) at the disease-free equilibrium W0 = (Sh∗ , 0, 0, 0, 0, 0), we obtain the following system: ⎧ E˙h (t) = cvh β(t)Iv (t) − (dh + α)Eh (t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ I˙h (t) = αEh (t) − (dh + κ + rh )Ih (t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ R˙h (t) = rh Ih (t) − (dh + γ)Rh (t), (24) ⎪ ⎪ ⎪ ⎪ ⎪ N ∗ (t) N ∗ (t) ⎪ E˙v (t) = chv β(t) Sv ∗ Ih (t) + c¯hv β(t) Sv ∗ Rh (t) − (νv + γ4 (t))Ev (t), ⎪ ⎪ ⎪ h h ⎪ ⎪ ⎪ ⎪ ⎩ ˙ Iv (t) = νv Ev (t) − γ4 (t)Iv (t). System (24) can be written as follows: U˙ 2 (t) =

(

) F2 (t) − V2 (t) U2 (t)

with U2 (t) = (Eh (t), Ih (t), Rh (t), Ev (t), Iv (t))T and ⎛

0

0

0

0 cvh β(t)

⎜ 0 0 0 0 ⎜ 0 ⎜ ⎜ 0 0 0 0 0 F2 (t) = ⎜ ∗ ∗ ⎜ N (t) N (t) ⎜ 0 chv β(t) v ∗ c¯hv β(t) Sv ∗ 0 0 Sh ⎝ h 0 0 0 0 0 ⎛ dh + α 0 0 0 ⎜ dh + κ + rh 0 0 ⎜ −α ⎜ ⎜ 0 −rh dh + γ 0 V2 (t) = ⎜ ⎜ 0 0 0 νv + γ4 (t) ⎝ 0

0

0

−νv

⎞ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎠

0

⎞

0

⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎠

0 0 γ4 (t)

Let us assume that Y˜ (t, s), t ≥ s is the matrix solution of the linear ω-periodic system { Y˜˙ (t, s) = −V2 (t)Y˜ (t, s), ∀t ≥ s, Y˜ (s, s) = I,

(25)

where I is the 5 × 5 identity matrix. Let Cω be the ordered Banach space of all ω-periodic functions from R to R5 which is equipped with the { } maximum norm ∥.∥ and the positive cone C+ ω := ϕ ∈ Cω : ϕ(t) ≥ 0, ∀t ∈ R . Then, we can define a linear operator L˜ : Cω −→ Cω by ∫ ∞ ˜ (Lϕ)(t) = Y˜ (t, t − a)F2 (t − a)ϕ(t − a)da, ∀t ∈ R, ϕ ∈ Cω . (26) 0

˜ the spectral radius of L. ˜ L˜ is the next infection operator, and the basic reproduction ratio is R0 = ρ(L),

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

19

3.3.3. Global stability of disease-free equilibrium Lemma 3.8 ([14]). (i) R0 = 1 if and only if ρ(ΦF2 −V2 (ω)) = 1. (ii) R0 < 1 if and only if ρ(ΦF2 −V2 (ω)) < 1. (iii) R0 > 1 if and only if ρ(ΦF2 −V2 (ω)) > 1. Hence, if R0 < 1 then the disease-free equilibrium, W0 is locally asymptotically stable. Theorem 3.4.

The disease-free equilibrium W0 is globally attractive if R0 < 1 and κ = 0.

Proof . If κ = 0, then from Eq. (7), there exists a period T (ϵ) > 0 such that Nh (t) ≥ Sh∗ − ϵ, ∀t ≥ T (ϵ). Thus, for all t ≥ T (ϵ), we have kv (t)Nv∗ (t) ≤ chv β(t)

Nv∗ (t) Nv∗ (t) I (t) + c ¯ β(t) Rh (t). h hv Sh∗ − ϵ Sh∗ − ϵ

Hence, from system (22) we obtain the following system: ⎧ E˙ h (t) ≤ cvh β(t)Iv (t) − (dh + α)Eh (t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ I˙h (t) = αEh (t) − (dh + rh )Ih (t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ R˙ h (t) = rh Ih (t) − (dh + γ)Rh (t), ⎪ ⎪ ⎪ ⎪ ⎪ N ∗ (t) N ∗ (t) ⎪ E˙ v (t) ≤ chv β(t) S ∗v −ϵ Ih (t) + c¯hv β(t) S ∗v −ϵ Rh (t) − (νv + γ4 (t))Ev (t), ⎪ ⎪ ⎪ h h ⎪ ⎪ ⎪ ⎪ ⎩ ˙ Iv (t) = νv Ev (t) − γ4 (t)Iv (t).

(27)

Now, let us consider the following auxiliary system, z˜˙ = Aϵ (t)˜ z (t)

(28)

with ˜h (t), I˜h (t), R ˜ h (t), E ˜v (t), I˜v (t))T z˜(t) = (E and

⎛ ⎜ ⎜ ⎜ ⎜ Aϵ (t) = ⎜ ⎜ ⎜ ⎝

−(dh + α)

0

0

0

cvh β(t)

⎞

α

−(dh + rh )

0

0

0

0

rh

−(dh + γ)

0

0

0

N ∗ (t) chv β(t) S ∗v −ϵ h

N ∗ (t) c¯hv β(t) S ∗v −ϵ h

−(νv + γ4 (t))

0

⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎠

−γ4 (t) ) ( ) Moreover, by continuity of spectral radius, we have limϵ→0+ ρ ΦAϵ (ω) = ρ ΦF2 −V2 (ω) . In addition, ( ) from Lemma 3.8, if R0 < 1, then ρ(ΦF2 −V2 (ω)) < 1, thus, limϵ→0+ ρ ΦAϵ (ω) < 1. It then follows that ( ) there exists ϵ1 > 0 such that ρ ΦAϵ (ω) < 1 for all ϵ ∈ [0, ϵ1 ]. From Lemma 3.1, there exists a positive ( ) ω-periodic function v(t) such that z˜(t) = v(t)ert is a solution of (28). Since ρ ΦAϵ (ω) < 1, so r < 0. 0

0

νv

0

(

20

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

Therefore, the ω-periodic function v(t) is bounded, hence z˜(t) → 0 as t → ∞. Applying the comparison theorem to system (27), we get lim (Eh (t), Ih (t), Rh (t), Ev (t), Iv (t)) = (0, 0, 0, 0, 0).

t→∞

Moreover, we have ( ) lim (Sh (t) − Sh∗ ) = lim Nh (t) − Eh (t) − Ih (t) − Rh (t) − Sh∗ = 0.

t→∞

t→∞

Hence, we have limt→∞ Sh (t) = Sh∗ . It then follows that the disease-free equilibrium W0 is globally attractive which completes the proof. □ 3.3.4. Persistence of malaria transmission Let us consider the following sets: Y := R6+ , { } Y0 := (Sh , Eh , Ih , Rh , Ev , Iv ) ∈ Y : Eh > 0, Ih > 0, Rh > 0, Ev > 0, Iv > 0 , { } ∂Y0 := (Sh , Eh , Ih , Rh , Ev , Iv ) ∈ Y : Eh Ih Rh Ev Iv = 0 . Let S2 : Y −→ Y be the Poincar´e map associated with system (22). Lemma 3.9. The sets Y0 and ∂Y0 are positively invariant with respect to system (22). Proof . For any initial condition φ2 ∈ Y0 , by solving system (22) we obtain, for all t > 0, ( ∫ t )[ (∫ τ ) ] ∫ t Sh (t) = exp − H2 (τ )dτ Sh (0) + H1 (τ ) exp H2 (s)ds dτ > 0, 0

0

(29)

0

[ ( ) ] ∫ t ( ) Eh (t) = exp −(dh + α)t Eh (0) + kh (τ )Sh (τ ) exp (dh + α)τ dτ > 0,

(30)

0

[ ( ) ] ∫ t ( ) Ih (t) = exp −(dh + κ + rh )t Ih (0) + αEh (τ ) exp (dh + κ + rh )τ dτ > 0

(31)

0

[ ( ) ] ∫ t ( ) Rh (t) = exp −(dh + γ)t Rh (0) + rh Ih (τ ) exp (dh + γ)τ dτ > 0,

(32)

0

( ∫ t )[ (∫ ∫ t Ev (t) = exp − H4 (τ )dτ Ev (0) + H3 (τ ) exp 0

0

τ

0

(33)

) ] γ4 (s)ds dτ > 0,

(34)

0

( ∫ t )[ (∫ ∫ t Iv (t) = exp − γ4 (τ )dτ Iv (0) + νv Ev (τ ) exp 0

) ] H4 (s)ds dτ > 0

τ

0

with H1 (t) = Λh + Ih (t) + γRh (t), H2 (t) = kh (t) + dh , ( ) H3 (t) = kv (t) Nv∗ (t) − Iv (t) , H4 (t) = kv (t) + γ4 (t) + νv . So, it then follows that Y0 is positively invariant. Therefore, ∂Y0 is relatively closed in Y , then, ∂Y0 is also positively invariant with respect to system (22). □

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

21

Moreover, the compact set ∆2 attracts all positive orbits in Y , which implies that the discrete-time system S2 : Y −→ Y is point dissipative. Lemma 3.10.

If R0 > 1, there exists ξ > 0 such that when ∥φ2 − W0 ∥ ≤ ξ, for any φ2 ∈ Y0 ,

one has ( ) lim sup d S2n (φ2 ), W0 ≥ ξ.

n→∞

Proof . Suppose by contradiction that ( ) lim sup d S2n (φ2 ), W0 < ξ for some φ2 ∈ Y0 .

n→∞

Then, it follows that there exists an integer N ≥ 1 such that for all n ≥ N we have ( ) d S2n (φ2 ), W0 < ξ. ( ) By the continuity of solution u2 t, φ2 , we have  ( n ( )) ( ) u2 t, S2 φ2 − u2 t, W0  ≤ ε, ∀t ≥ 0, ε > 0. [ ] [ ] Let t = nω + t′ , with t′ ∈ [0, ω] and n = ωt . ωt is the greatest integer less than or egal to have  ( ) ( )  ( ( )) ( ) u2 t, φ2 − u2 t, W0  = u2 t′ , S2n φ2 − u t′ , W0  < ε.

t ω.

Then, we

Hence, it follows that Sh∗ − ε ≤ Sh (t) ≤ Sh∗ + ε. Thus, there exists ζ(ε) > 0 such that Sh (t) ≥ 1 − ζ(ε) Nh (t)

and

Nv∗ (t) N ∗ (t) ≥ v∗ . Nh (t) Sh

From system (22) we have: ⎧ E˙h (t) ≥ cvh β(t)(1 − ζ(ε))Iv (t) − (dh + α)Eh (t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ I˙h (t) = αEh (t) − (κ + dh + rh )Ih (t), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ˙ Rh (t) = rh Ih (t) − (dh + γ)Rh (t), ⎪ ⎪ ( ) ⎪ ⎪ ⎪ Nv∗ (t) ⎪ ˙ Ev (t) ≥ β(t) S ∗ chv Ih (t) + c¯hv Rh (t) − (νv + γ4 (t))Ev (t), ⎪ ⎪ ⎪ h ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ˙ Iv (t) = νv Ev (t) − γ4 (t)Iv (t). Now, let us consider the following auxiliary linear system ˜˙ ˜ h(t) = Nζ(ε) (t)h(t) with ˜ = (E ˜h (t), I˜h (t), R ˜ h (t), E ˜v (t), I˜v (t))T h(t)

(35)

22

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

and ⎛ ⎜ ⎜ ⎜ ⎜ Nζ(ε) (t) = ⎜ ⎜ ⎜ ⎜ ⎝

−(dh + α)

0

0

0

cvh β(t)(1 − ζ(ε))

⎞

α

−(dh + rh )

0

0

0

⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎠

−(dh + γ) 0 0 Nv∗ (t) 0 c¯hv β(t) −(νv + γ4 (t)) 0 Sh∗ 0 0 νv −γ4 (t) ( ) Once again if R0 > 1, there exists σ > 0 such that ρ ΦNζ(ε) (ω) > 1 for all ζ(ε) ∈ [0, σ]. In this case ˜ r > 0 and from Lemma 3.1 we get h(t) → ∞ as t → ∞. Thus, by applying the theorem of comparison, we get lim |(Eh (t), Ih (t), Rh (t), Ev (t), Iv (t))| = ∞. rh Nv∗ (t) chv β(t) Sh∗ 0

0

t→∞

This leads to a contradiction. □ Theorem 3.5. If R0 > 1, then system (22) is uniformly persistent and there exists at least one positive periodic solution. Proof . Set { } M∂ = φ2 ∈ ∂Y0 : S2n (φ2 ) ∈ ∂Y0 , for any n ≥ 0 . Now, we aim to show that { } M∂ = (Sh , 0, 0, 0, 0, 0) ∈ Y : Sh ≥ 0 .

(36)

It is obvious that { } M∂ ⊃ (Sh , 0, 0, 0, 0, 0) ∈ Y : Sh ≥ 0 , so, we only need to show that { } M∂ ⊂ (Sh , 0, 0, 0, 0, 0) ∈ Y : Sh ≥ 0 , that means that for any initial condition φ2 ∈ ∂Y0 and n ≥ 0, we have Eh (nω) = 0 or Ih (nω) = 0 or Rh (nω) = 0 or Ev (nω) = 0 or Iv (nω) = 0. Let φ2 ∈ ∂Y0 . Suppose by contradiction that there exists an integer n1 ≥ 0 such that Eh (n1 ω) > 0, Ih (n1 ω) > 0, Rh (n1 ω) > 0, Ev (n1 ω) > 0 and Iv (n1 ω) > 0. Thus, by replacing the initial time t = 0 by t = n1 ω in (29)–(34) we get Sh (t) > 0, Eh (t) > 0, Ih (t) > 0, Rh (t) > 0, Ev (t) > 0, Iv (t) > 0. That { } contradicts the fact that ∂Y0 is positively invariant. Hence, M∂ ⊂ (Sh , 0, 0, 0, 0, 0) ∈ Y : Sh ≥ 0 and then the equality (36) is true. The equality (36) implies that W0 is a fixed point of S2 and acyclic in M∂ , every solution in M∂ approaches to W0 . Moreover, Lemma 3.10 implies that W0 is an isolate invariant set in Y and W s (W0 ) ∩ Y0 = ∅. By Theorem 3.1.1 in [22], it then follows that S2 is uniformly persistent with respect to (Y0 , ∂Y0 ). So, the solution of system (22) is uniformly persistent. It yields that system (22) has at least one periodic solution u∗2 (t, φ∗2 ) with φ∗2 ∈ Y0 . We further claim that u∗2 (t, φ∗2 ) is positive for all t ≥ 0. Suppose that u∗2 (0, φ∗2 ) = 0. Then, from Eqs. (29)–(34) we have u∗2 (t, φ∗2 ) > 0 for all t > 0. By the periodicity of the solution, we have u∗2 (0 + nω, φ∗2 ) = u∗2 (nω, φ∗2 ) = u∗2 (0, φ∗2 ) = 0 for all n ≥ 1. That leads to a contradiction. □ 3.4. Analysis of the main model Now, we use the theory of chain transitive set [22] to lift the threshold type result for system (22) to system (6).

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

23

3.4.1. Positivity and boundedness of solutions Lemma 3.11. For any positive initial condition φ = (φ1 , φ2 ), system (6) has a unique positive bounded solution u(t, φ) on [0, +∞) with u(0) = φ. Moreover, the closed set ∆ = ∆1 × ∆2 is positively invariant for system (6). ( )T Proof . For any φ, we have f (t, φ) = f1 (t, φ1 ), f2 (t, φ2 ) is continuous in (t, φ) and lipschitzian in φ in each compact of ∆. So, thanks to Theorems 2.2.1 and 2.2.3 in [21], system (6) has a unique solution v(t, φ) on its maximal interval [0, σφ ) of existence. Since ∆1 and ∆2 are both respectively positively invariant under (10) and (22), then ∆ is positively invariant for (6). Hence, the solution u(t, φ) is bounded. It then follows that σφ = +∞. □ 3.4.2. Stability of equilibria Theorem 3.6. for system (6).

( ) If Rv < 1, then the disease-free equilibrium 0, 0, 0, 0, Sh∗ , 0, 0, 0, 0, 0 is globally attractive

˜ Proof . Let Φ(t) the periodic semiflow generated by the periodic system (6) and S the Poincar´e map associated with system (6), namely: ˜ S(φ) = Φ(ω)φ = v(ω, φ), n ˜ S (φ) = Φ(nω)ϕ = v(nω, φ), ∀n ≥ 0. ˜ n (φ)}n≥0 . Then, from Lemma 1.2.1 in [22], O Let O be the omega limit set of the discrete-time orbit {Φ n ˜ on ∆. If Rv < 1, then thanks to Theorem 3.2, we have is an internally chain transitive set for Φ ( ) lim E(nω, φ), L(nω, φ), P (nω, φ), Nv (nω, φ), Ev (nω, φ), Iv (nω, φ) = (0, 0, 0, 0, 0, 0).

n→∞

Then, there exists O′ ⊂ R4 such that O = {(0, 0, 0, 0, 0, 0)} × O′ . For any x = (xi )1≤i≤n ∈ O, there exists a sequence nk → ∞ such that u(nk ω, φ) → x as k → ∞. Since, 0 ≤ Sh (nk ω, φ) ≤ Nh∗ ,

0 ≤ Eh (nk ω, φ) ≤ Nh∗ ,

0 ≤ Ih (nk ω, φ) ≤ Nh∗ ,

0 ≤ Rh (nk ω, φ) ≤ Nh∗ ,

then O′ ⊂ ∆2 . Moreover, it is obvious that ˜ n (0, 0, 0, 0, 0, 0, x5 , x6 , x7 , x8 ) = {(0, 0, 0, 0, 0, 0)} × Ψ n ′ (x5 , x6 , x7 , x8 ) Φ |O |O where {Ψ (t)}t≥0 is the semiflow associated to the following system: ⎧ ˙ Sh (t) = Λh + γRh (t) − dh Sh (t), ⎪ ⎪ ⎪ ⎨ E˙ (t) = −(d + α)E (t), h h h ˙ ⎪ I (t) = αE (t) − (d h h + κ + rh )Ih (t), ⎪ h ⎪ ⎩ ˙ Rh (t) = rh Ih (t) − (dh + γ)Rh (t).

(37)

˜ n , it then follows that O′ is internally chain transitive set for Since O is internally chain transitive set for Φ n ˜ Ψ . In addition, from the second equation of system (37), we have: lim Eh (t) = lim e−(dh +α)t = 0.

t→∞

t→∞

24

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

Thus, we obtain the following system: {

I˙h (t) ≤ αEh (t), R˙ h (t) ≤ rh Ih (t).

(38)

Using the comparison theorem, on system (38) we have: ( ) lim Ih (t), Rh (t) = (0, 0). t→∞

It yields from the first equation of system (37) that limt→∞ Sh (t) = Λd h = Sh∗ . Hence, (Sh∗ , 0, 0, 0, ) is h globally attractive for system (37). So, from Theorem 1.2.1 in [22], we have O′ = (Sh∗ , 0, 0, 0, ) and then O = (0, 0, 0, 0, Sh∗ , 0, 0, 0, 0, 0). □ (H8): Assume that the positive periodic solutions of systems (10) and (22) are globally attractive. Theorem 3.7. Let (H8) hold. If κ = 0, Rv > 1 and R0 < 1, then the disease-free periodic equilibrium ) ( ∗ E (t), L∗ (t), P ∗ (t), Nv∗ (t), Sh∗ , 0, 0, 0, 0, 0 is globally attractive for system (6). Proof . If Rv > 1, then thanks to Theorem 3.3, we have ( ) ( ) lim E(nω, φ), L(nω, φ), P (nω, φ), Nv (nω, φ) = E ∗ (0), L∗ (0), P ∗ (0), Nv∗ (0) . n→∞

) ( Then, there exists O′′ ⊂ R6 such that O = { E ∗ (0), L∗ (0), P ∗ (0), Nv∗ (0) } × O′′ . For any y = (yi )1≤i≤10 ∈ O, there exists a sequence nl → ∞ such that v(nl ω, φ) → y as l → ∞. Since, 0 ≤ Sh (nl ω, φ) ≤ Nh∗ ,

0 ≤ Eh (nl ω, φ) ≤ Nh∗ ,

0 ≤ Ih (nl ω, φ) ≤ Nh∗ ,

0 ≤ Rh (nl ω, φ) ≤ Nh∗ ,

0 ≤ Ev (nl ω, φ) ≤ Nv∗ (nl ω, φ),

0 ≤ Iv (nl ω, φ) ≤ Nv∗ (nl ω, φ),

then O′′ ⊂ ∆2 . Moreover, let y ′ = (y5 , y6 , y7 , y8 , y9 , y10 ). It is obvious that {( )} ˜ n (E ∗ (0), L∗ (0), P ∗ (0), Nv∗ (0), y ′ ) = E ∗ (0), L∗ (0), P ∗ (0), Nv∗ (0) × Ψ ˆ n ′′ (y ′ ), Φ |O |O ˆ (t)}t≥0 is the semiflow associated to system (22). Hence, from Theorem 3.4, if R0 < 1 and κ = 0 where {Ψ then, the equilibrium W0 is globally attractive. So, from Theorem 1.2.1 in [22], we have O′′ = (Sh∗ , 0, 0, 0, 0, 0) and then O = (E ∗ (0), L∗ (0), P ∗ (0), Nv∗ (0), Sh∗ , 0, 0, 0, 0, 0). □ Theorem 3.8. Let (H8) hold. If Rv > 1 and R0 > 1, then the endemic periodic equilibrium ( ∗ ) E (t), L∗ (t), P ∗ (t), Nv∗ (t), Sh∗ (t), Eh∗ (t), Ih∗ (t), Rh∗ (t), Ev∗ (t), Iv∗ (t) is globally attractive for system (6). Proof . If Rv > 1, then thanks to Theorem 3.3, we have ( ) ( ) lim E(nω, φ), L(nω, φ), P (nω, φ), Nv (nω, φ) = E ∗ (0), L∗ (0), P ∗ (0), Nv∗ (0) . n→∞

( ) Then, there exists T ⊂ R6 such that O = { E ∗ (0), L∗ (0), P ∗ (0), Nv∗ (0) } × T . Thus, from Hirsch et al. [43] (see Theorem 3.1), ( ) T = Sh∗ (0), Eh∗ (0), Ih∗ (0), Rh∗ (0), Ev∗ (0), Iv∗ (0) or T = (0, 0, 0, 0, 0, 0). Let us suppose that T = (0, 0, 0, 0, 0, 0). Then, it yields that u2 (t, φ2 ) → 0 as t → ∞. Hence, there exists a period ω such that for any ϵ∗ > 0 and t ≥ ω, (  )  E(t), L(t), P (t), Nv (t) − (E ∗ (t), L∗ (t), P ∗ (t), Nv∗ (t)) < ϵ∗ .

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

25

Hence, from system (22) we have ⎧ ∗ ˙ ⎪ ⎪ Eh (t) ≥ cvh β(t)(1 − ϵ )Iv (t) − (dh + α)Eh (t), ⎪ ⎪ ⎪ ⎪ I˙h (t) = αEh (t) − (κ + dh + rh )Ih (t), ⎪ ⎪ ⎨ ˙ Rh (t) = rh Ih (t) − (dh + γ)Rh (t), ( ) ⎪ Nv∗ (t) ⎪ ⎪ ˙ chv Ih (t) + c¯hv Rh (t) − (νv + γ4 (t))Ev (t), Ev (t) ≥ β(t) ∗ ⎪ ⎪ ⎪ Sh ⎪ ⎪ ⎩ ˙ Iv (t) = νv Ev (t) − γ4 (t)Iv (t). ( If R0 > 1, then, from Theorem 3.5, system (22) admits a globally attractive fixed point Sh∗ (0), Eh∗ (0), ) Ih∗ (0), Rh∗ (0), Ev∗ (0), Iv∗ (0) ≫ 0 for all φ2 in Y0 . So, system (22) is uniformly persistent, which leads to a ( ) contradiction. It then follows that T = E ∗ (0), L∗ (0), P ∗ (0), Nv∗ (0) . □ To finish this Section, we remark that when all coefficients are constant, system (5) reduces to an autonomous system. In this case, the vector reproduction ratio Rv and the basic reproduction ratio, R0 read as follows: bα1 α2 α3 , γ (γ + α1 )(γ2 + α2 )(γ3 + α3 ) √4 1 ( )( ) cvh ανv β 2 α3 KP 1 R0 = c¯hv rh + chv (dh + γ) 1 − . Q1 Q2 Rv

Rv =

(39) (40)

with Q1 = Sh∗ γ4 (νv + γ4 )(dh + γ)(dh + κ + rh )(dh + α), α1 KE + KL (α2 + γ2 ) Q2 = 1 + (α3 + γ3 ) . α1 α2 KL KE It should be noted from equality (40) that R0 exists if and only if Rv > 1. Moreover, it is obvious that the basic reproduction ratio R0 depends on the vector reproduction ratio Rv . Indeed, larger Rv is, larger R0 becomes. This result shows the importance to take into account the immature mosquitoes in the model because this parameter can be used as a control to eliminate the disease. 4. Numerical simulations In this section, we present a series of numerical simulations in order to support our theoretical results, to predict the trend of the disease, and to explore some control measures. We use the MATLAB technical computing software with the fourth-order Runge–Kutta method to perform our numerical results. 4.1. Estimation of the parameters We suppose that the life expectancy of humans for Burkina Faso in 2011 was 55 years. Thus, the human natural death rate dh can be calculated as follows: dh =

1 = 0.00152 per month. 12 × 55

In addition, in Burkina Faso, the total size of human population in 2011 was 16, 248, 558, then, the recruitment rate Λh is determined by Λh = dh × 16, 248, 558 = 24, 698 humans per month.

26

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

Table 3 Values for constant parameters for malaria model. Parameters

Values

Sources

Dimensions

Λh

24,698 0.00152 0.0028 3.04 0.0833 0.0146 2.523 0.022 0.48 0.048 8 000 000 6 000 000 5 500 000

See text See text [17] [17] [17] [17] [17] [9] [9] [9] [25] [25] [25]

/month /month /month /month /month /month /month – – – – – –

dh dp α rh γ νv cvh chv c ¯hv KE KL KP

Table 4 Monthly mean temperature of Burkina Faso (in

◦

C).

Month Temperature

January 25.8

February 28.77

March 31.11

April 30.75

May 29.54

June 27.50

Month Temperature

July 25.95

August 25.04

September 25.57

October 27.32

November 28.33

December 25.36

Further, values of some constant parameters for malaria transmission model (6) are listed in Table 3 with the references therein. Now, we estimate the periodic coefficients of the model. Thus, using the reported monthly mean temperature, for ASECNA in Burkina Faso, from January 2011 to December 2015, we can calculate the monthly mean temperature of Burkina Faso in Table 4. Moreover, assuming that the average number of mosquito bites depends on their gonotrophic cycle which is also a function of temperature θ, then, the temperature-dependent biting rate of mosquitoes reads as follows [18]: √ β(θ) = 0.00609T (θ − 11.7) 42.3 − θ per month. Similarly, the temperature-dependent death rate of adult mosquitoes and larvae reads respectively as follows [17]: ( ) 5−θ γ4 (θ) = 3 + 29.564 exp per month, 2.7035 30 γ2 (θ) = per month. −4.4 + 1.31θ − 0.03θ2 By considering that the egg-laying rate is proportional to the biting rate, eggs and pupae death rate are proportional to the death rate of larvae, then we have: b(T ) = 2.5β(θ), γ1 (θ) = 0.5γ2 (θ) and γ3 (θ) = 0.75γ2 (θ). Similarly, we have: α1 (θ) = 0.95b(θ), α2 (θ) = 0.5α1 (θ) and α3 (θ) = 0.8α2 (θ). Furthermore, assuming that temperature varies as a function of time, then we use the CFTOOL of MATLAB to get the best estimation of periodic functions β(t), γ4 (t), γ2 (t) of Burkina Faso as follows: ( ) ( )] 3 [ ∑ nπt nπt +a ˜n cos per month β(t) = a0 + an sin 6 6 n=1 with a0 = 10.16, a1 = 1.289,

a2 = −0.06678, a3 = 0.3301,

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

a ˜1 = 0.3614, a ˜2 = −1.097,

27

a ˜3 = −0.4723.

) ( )] ( 3 [ ∑ nπt nπt ˜ γ2 (t) = b0 + + bn cos per month bn sin 6 6 n=1 with b0 = 3.443, b1 = 0.3287, b2 = 0.03792, b3 = −0.004983, ˜b1 = 0.08803, ˜b2 = −0.2627, ˜b3 = −0.1211. γ4 (t) = c0 +

) ( )] ( 3 [ ∑ nπt nπt + c˜n cos per month cn sin 6 6 n=1

with c0 = 3.008, c1 = −0.004631, c2 = 0.0009815, c3 = −0.0026, c˜1 = −0.001408, c˜2 = 0.004167,

c˜3 = 0.001817.

4.2. Numerical results We use the values of Table 3 and the initial conditions: E(0) = 600, 000, L(0) = 300, 000, P (0) = 200, 000, Nv (0) = 24, 000, Sh (0) = 15, 000, 000, Eh (0) = 10, 000, Ih (0) = 50, 000, Rh (0) = 1, 188, 558, Ev (0) = 10, 000, Iv (0) = 5, 000 to get the following results: Figs. 4–7. Computing the two thresholds dynamics of the model, we get Rv = 4.8413 and R0 = 1.0364. Through Fig. 4, we remark that mosquitoes and malaria persist. In addition systems (10) and (22) admit respectively positive periodic solutions (E ∗ (t), L∗ (t), P ∗ (t), Nv∗ (t)) and (Sh∗ (t), Eh∗ (t), Ih∗ (t), Rh∗ (t), Ev∗ (t), Iv∗ (t)) which are globally asymptotically stable. That observation illustrates the results of our Theorems 3.3 and 3.5. Since the death rate induced by malaria in the country is very high, the government of Burkina Faso has decided to introduce some measures of control in order to slow down the disease evolution. These control measures concern three major points: vector control, prevention and treatment. Indeed, the vector control consists in reducing the breeding sites of eggs, larvae and pupae KE , KL , Kp , respectively in addition to the use of insecticides against adult mosquitoes. Concerning the prevention, it consists in avoiding the biting of mosquitoes β(t) by using bed nets. Besides, the health authorities have made available to the infectious populations some drugs to treat malaria. Therefore, we will study the effectiveness of prevention and vector control in the eradication of malaria. Firstly, we examine the impact of vector control on the transmission. Let us assume that 25% of eggs breeding sites KE , 34% of larvae breeding sites KL and 37% of pupae breeding sites KP have been destroyed by human populations through environmental sanitations. Then, we have the following results: Fig. 5. Computing the vector reproduction ratio and the basic reproduction ratio, we obtain Rv = 4.8413 and R0 = 0.8465. We remark that despite the persistence of mosquitoes population (see Fig. 5(a)), the disease will die out to both populations (see Fig. 5((c)–(f))). Furthermore, system (6) converges to the ( ) periodic disease-free equilibrium E ∗ (t), L∗ (t), P ∗ (t), Nv∗ (t), Sh∗ , 0, 0, 0, 0, 0 which is globally attractive. That illustrates the result of our Theorem 3.7. Now, let us replace the vector control by the prevention one. Thus, we assume that humans use bed nets to avoid the biting of mosquitoes. Hence, if the use of bed nets allows to reduce the biting rate to 25% then we obtain the following results: Fig. 6. Computing the two thresholds dynamics of the model, we get Rv = 3.0907 and R0 = 0.7612. We note that with a contact reduction rate of 25%, malaria will disappear from the country in the next thirty years (see Fig. 6((c)–(f))). Thus, increasing this control would be a very effective means for a short-term eradication of

28

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

Fig. 4. Global behavior of system (6) for Rv > 1 and R0 > 1.

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

29

Fig. 5. Global behavior of system (6) for KE = 6000 000, KL = 4000 000, KP = 3500 000, Rv > 1 and R0 < 1.

malaria in Burkina Faso. As a result, health authorities should be more involved in mass use of bed nets by populations. However, in practice, the financial cost of this control is very high for the country. Therefore, it would be more appropriate to combine the two types of control in order to obtain a more optimal result. Computing the vector reproduction ratio, we get Rv = 0.9551. Thus, when we reduce the contact rate by 61% and the larval sites by 50%, we notice that in less than 3 years the mosquitoes disappear from the country (see Fig. 7((a)–(c)-(d)). Since, Rv < 1, then system (6) converges globally to the disease-free equilibrium (0, 0, 0, 0, Sh∗ , 0, 0, 0, 0, 0). That illustrates our Theorem 3.6. 5. Conclusion An important issue in developing mathematical models is to identify which biological factors are necessary to include and which ones can be omitted. Thus, considering the impact of climate change on the mosquito

30

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

Fig. 6. Global behavior of system (6) for Rv > 1 and R0 < 1.

life cycle, we have formulated a deterministic mathematical model of malaria transmission by incorporating seasonality. The vector population has been divided into juvenile stage and adult stage. This structure enabled us to better understand the dynamics of mosquito population and malaria transmission dynamics. In addition, taking into account the juvenile stage offers many advantages in the control of the disease transmission. Two thresholds dynamics have been determined for the model: the vector reproduction ratio Rv and the basic ratio R0 and it emerges from our analysis that the global behavior of the model depends on those thresholds. Hence, we have shown that if: (i) Rv < 1, mosquitoes disappear. (ii) Rv > 1, and R0 < 1, mosquitoes persist but the disease dies out. (iii) Rv > 1, and R0 > 1 mosquitoes and disease persist with seasonal variation.

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

31

Fig. 7. Global behavior of system (6) for KE = 4000 000, KL = 3000 000, KP = 2750 000 and Rv = 0.9551.

Using the monthly mean temperature of Burkina Faso and the formula related to mosquito life cycle, we estimated the periodic coefficients of the model and then verified our theoretical results by numerical simulations. Our numerical results indicated that there exist two important parameters which permit to control malaria transmission in the country. The first control is the reduction of the periodic biting rate β(t) and the second control concerns the reduction of carrying capacities of immature mosquitoes KE , KL and KP . Indeed, if we reduce the biting rate to 25%, mosquitoes will persist but the disease will die out in the long run. That observation is the same if we reduce eggs, larvae and pupae breeding sites respectively to 25%, 34% and 37%. We remark that the reduction of carrying capacities has a great impact on the disease transmission. First, it reduces the transfer rates of immature stages and increases their death rates, because if there is no water, eggs cannot hatch. Furthermore, we showed that combining the two controls is more efficient and more optimal in the fight against malaria transmission. Following the previous remarks, the eradication of the disease can be achieved if the mosquito population is eradicated or if the basic reproduction ratio is less than one. However, even if there is no malaria epidemic, we need to be vigilant because climate change can lead to the re-emergence of malaria in some European and American countries. In our model, we neglected the effects of climate on the mosquitoes breeding sites. However, it must be noticed that these parameters depend on the climate profile. So, it would be more interesting to consider them as functions of time. In addition, due to the complexity of ours systems, we have not established the

32

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

global stability of the positive periodic solution (see hypothesis (H8)), nevertheless, numerical simulations indicate it clearly. As a future work we will investigate this open problem. Acknowledgments The authors are grateful to the editor and the anonymous referees for their insightful comments and suggestions which contributed to improve the original version of the manuscript. This work has been prepared with the support of the Department of Mathematics, Nazi BONI University. Disclosure The authors declare that there is no conflict of interest regarding the publication of this paper. References [1] J. Wang, Z. Teng, T. Zhang, Thresholds dynamics of a malaria transmission model in periodic environment, Commun. Nonlinear Sci. Numer. Simul. 18 (5) (2013) 1288–1303. [2] World Health Organisation (WHO), Global Malaria Programme, World Malaria Report, 2015. [3] Z. Ma, J. Li, Dynamical Modeling and Analysis of Epidemics, World Scientific Publishing Co. Pte. Ltd, Singapore. [4] S.E. Eikenberry, A.B. Gumel, Mathematical modeling of climate change and malaria transmission dynamics: a historical review, J. Math. Biol. 77 (4) (2018) 857–933. [5] J.C. Koella, On the use of mathematical models of malaria transmission, Acta Trop. 49 (1991) 1–25. [6] R. Ross, The Prevention of Malaria, John Murray, London, 1911. [7] G. Macdonald, The Epidemiology and Control of Malaria, Oxford University Press, London, 1957. [8] C. Chiyaka, J.M. Tchuenche, W. Garira, S. Dube, A mathematical analysis of the effects of control strategies on the transmission dynamics of malaria, Appl. Math. Comput. 195 (2) (2008) 641–662. [9] N. Chitnis, J.M. Hyman, J.M. Cushing, Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model, Bull. Math. Biol. 70 (5) (2008) 1272–1296. [10] G.A. Ngwa, W.S. Shu, A Mathematical model for endemic malaria with variable human and mosquito populations, Math. Comput. Modelling 32 (7–8) (2000) 747–763. [11] B. Traor´ e, B. Sangar´ e, S. Traor´ e, A mathematical model of malaria transmission in a periodic environment, J. Biol. Dyn. 12 (1) (2018) 400–432. [12] A. Ai, J. Lu, J. Li, Mosquito stage-structured malaria models and their global dynamics, SIAM J. Appl. Math. 72 (4) (2011) 1213–1237. [13] O. Koutou, B. Traor´ e, B. Sangar´ e, Mathematical modeling of malaria transmission global dynamics: taking into account the immature stages of the vectors, Adv. Difference Equ. 2018 (220) (2018) 1–34. [14] Y. Nakata, T. Kuniya, Global dynamics of a class of SEIRS epidemic models in a periodic environment, J. Math. Anal. Appl. 363 (1) (2010) 230–237. [15] X. Zhang, J. Jia, X. Song, Permanence and extinction for a nonautonomous malaria transmission model with distributed time delay, J. Appl. Math. 2014 (2014) 1–15, no. ID 139046. [16] K. Okuneye, A.B. Gumel, Analysis of a temperature-rainfall-dependent model for malaria transmission dynamics, Math. Biosci. 287 (2017) 72–92. [17] Y. Lou, X.-Q. Zhao, A climate-based malaria transmission model with structured vector population, SIAM J. Appl. Math. 70 (6) (2010) 2023–2044. [18] X.-Q. Wang, X. Zhao, A climate-based malaria model with the use of bed nets, J. Math. Biol. 77 (1) (2018) 1–25. [19] B. Traor´ e, B. Sangar´ e, S. Traor´ e, A mathematical model of malaria transmission withs structured vector population and seasonality, J. Appl. Math. 2017 (2017) 1–15, no. ID 6754097. [20] J. Tiana, J. Wang, Some results in Floquet theory, with application to periodic epidemic models, Appl. Anal. 94 (6) (2015) 1–25. [21] J.K. Hale, S.M. Verduyn Lunel, Introduction to Functional Differential Equations, Springer, New York. [22] X.-Q. Zhao, Dynamical systems in population biology, in: CMS Books in Mathematics/Ouvrages de Math´ ematiques de la SMC, Vol. 16, Springer-Verlag, New York, 2003. [23] L.M. Beck-Johnson, W.A Nelson, K.P. Paaijmans, A.F. Read, et al., The effect of temperature on Anopheles mosquito population dynamics and the potential for malaria transmission, PLoS One 8 (11) (2013) 1–12. [24] J. Lu, J. Li, Dynamics of stage-structured discrete mosquito population dynamics, J. Appl. Anal. Comput. 1 (1) (2011) 53–67. [25] D. Moulay, M.A. Aziz Alaoui, M. Cadivel, The chikungunya disease : Modeling, vector and transmission global dynamics, Math. Biosci. 229 (2011) 50–63. [26] A. Abdelrazec, A.B. Gumel, Mathematical assessment of the role of temperature and rainfall on mosquito population dynamics, J. Math. Biol. 74 (6) (2017) 1351–1395.

B. Traoré, O. Koutou and B. Sangaré / Nonlinear Analysis: Real World Applications 53 (2020) 103081

33

[27] A.M. Lutambi, M.A. Penny, T. Smith, N. Chitnis, Mathematical modelling of mosquito dispersal in a heterogeneous environment, Math. Biosci. 241 (2) (2013) 198–216. [28] J. Li, Simple stage-structured models for wild and transgenic mosquito populations, J. Difference Equ. Appl. 15 (4) (2009) 327–347. [29] F.P. Amerasinghe, T.S.B. Alagoda, Mosquito oviposition in bamboo traps, with special reference to Aedes albopictus, Aedes novalbopictus and Armigeres subalbatus, Int. J. Trop. Insect Sci. 5 (6) (1984) 493–500. [30] F. Pag` es, V. Corbel, C. Paupy, Aedes albopictus : Chronique d’un vecteur expansionniste, Med. Trop. 66 (3) (2006) 226–228. [31] J. Beehler, J. Millar, M. Mulla, Protein Hydrolysates and associated contamitants as ovipostion attractants for the mosquito, Med. Vet. Entomol. 8 (1994) 381–385. [32] J.A. Pickett, C.M. Woodcock, The role of mosquito olfaction in oviposition site location and in the avoidance of unsuitable hosts, in: Ciba Foundation Symposium 200-Olfaction in Mosquito-Host Interactions, 2007, pp. 109–123. [33] N. Chitnis, J.M. Cushing, J.M. Hyman, Bifurcation analysis of a mathematical model for malaria transmission, SIAM J. Appl. Math. 67 (1) (2006) 24–45. [34] P. Roop-O, W. Chinviriyasit, S. Chinviriyasit, The effect of incidence function in backward bifurcation for malaria model with temporaryimmunity, Math. Biosci. 265 (2015) 47–64. [35] H.L. Smith, P. Waltman, The Theory of the Chemostat, Cambridge University Press, 1995. [36] F. Zhang, X.-Q. Zhao, A periodic epidemic model in a patchy environment, J. Math. Anal. Appl. 325 (1) (2007) 496–516. [37] J.M. Bony, Principe du maximum,in´ egalit´ e de harnack et unicit´ e du probl` eme de Cauchy pour les op´ erateurs elliptiques d´ eg´ en´ er´ es, Ann. Inst. Fourier (Grenoble) 19 (1969) 277–304. [38] B. Traor´ e, O. Koutou, B. Sangar´ e, Global dynamics of a seasonal mathematical model of schistosomiasis transmission with general incidence function, J. Biol. Systems 27 (1) (2019) 19–49. [39] W. Wang, X.-Q. Zhao, Threshold dynamics for compartmental epidemic models in periodic environments, J. Dynam. Differential Equations 20 (3) (2008) 699–717. [40] P. Van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (2002) 29–48. [41] V. Lakshmikantha, S. Leela, A.A. Martynyuk, Stability Analysis of Nonlinear Systems, Marcel Dekker Inc., New York, Basel, 1989. [42] Z. Teng, Y. Liu, L. Zhang, Persistence and extinction of disease in non-autonomous SIRS epidemic models with disease-induced mortality, Nonlinear Anal. TMA 69 (8) (2008) 2599–2614. [43] M.W. Hirsch, H.L. Smith, X.-Q. Zhao, Chain transitivity, attractivity and strong repellors for semidynamical systems, J. Dynam. Differential Equations 13 (1) (2001) 107–131.