Electronic Journal of Differential Equations, Vol. 2022 (2022), No. 28, pp. 1–22. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu MATHEMATICAL MODELS FOR THE TRANSMISSION OF MALARIA WITH SEASONALITY AND IVERMECTIN ZHIHONG ZHAO, SHAOCHUN LI, YULAN LU Abstract. Ivermectin has shown good effects for malaria control in clinical trial stages because it can kill mosquitoes feeding on recently treated individu- als. In this article, we formulate and analyze a novel delay malaria transmission model taking into account seasonality and ivermectin. We show that the dy- namics of the model is totally determined by the basic reproduction ratio R0; that is, malaria will gradually die out if R0 < 1 and will persist if R0 > 1. Numerically, we verify the obtained theoretical results and evaluate the effect of ivermectin by related data of Kenya. We find that our simulation of the impact agrees with the prediction of the existing clinical trials in which it takes at least 25 years to eliminate malaria from Kenya with malaria control measures intact. 1. Introduction Malaria is an acute febrile disease caused by Plasmodium microorganisms spread to humans by infected adult female Anopheles mosquitoes. The disease accounts for 241 million infectious in 87 malaria endemic countries and 627,000 deaths world- wide, with about 80% of malaria deaths occurring in children under 5 years of age in 2020 [33]. Moreover, 95% and 96% of malaria cases and deaths occur in Africa. Mosquitoes as a major vector of malaria transmission have long been of interest to entomologists. Since the life cycle of mosquitoes is strongly correlated with season, then the trend toward malaria is most likely to follow the climate pattern [1, 10]. For instance, warmer temperature, which increases mosquito activity and lifespan, lead to mosquito bites more frequently. In addition, the temperature sensitivity of malaria parasites to mosquito hosts has long been established [15, 19]. After a long period of anti-malarial interventions, we found that Plasmodium falciparum infection prevalence in endemic Africa halved and the incidence of clin- ical disease fell by 40% between 2000 and 2015. Indoor residual spraying and insecticide-treated nets, the most widespread intervention, were by far the largest contributor (68% of cases averted) from 2000 to 2015 [4]. But still below target levels, there is an urgent need for additional tools to treat and control malaria, with the development of Plasmodium resistance to insecticides and drugs [20]. Ivermectin is a blocking drug that targets the vector itself, it works primarily by binding to glutamate-gated chlorine channels in nerve and muscle, leading to 2020 Mathematics Subject Classification. 92B05, 37N25. Key words and phrases. Malaria model; ivermectin; time delay; basic reproduction ratio; sensitivity analysis. ©2022. This work is licensed under a CC BY 4.0 license. Submitted November 16, 2021. Published April 7, 2022. 1 2 Z. ZHAO, S. LI, Y. LU EJDE-2022/28 hyperpolarization, paralysis and death of the invertebrate, including mosquitoes, and might inhibiting sporozoite development [6, 7]. It is the only avermectin class of endectocides that is available for human use from 1987 [6], and more than 416.8 million treatments have been distributed for mass drug administrations (MDA) to eliminate onchocerciasis and lymphatic filariasis in 2020 [16]. Higher doses is well tolerated in human beings up to 2000 µg/kg. It can be seen that ivermectin has an excellent safety profile [25]. Several experiments have shown that the mosquito- cidal effect of 150-200 µg/kg single doses of ivermectin are short-lived around 5-6 days [7]. Three doses of 300 µg/kg given over 3 days has a mosquitocidal effect in humans for 28 days against Anopheles gambiae sensu stricto [25]. It can be seen through those studies that ivermectin can reduce lifespan of mosquitoes, that means they are less likely to live long enough to complete sporogony and become in- fectious. Furthermore, unlike traditional vector control tools (e.g. indoor residual spraying and long-lasting insecticidal nets), ivermectin can reduce the likelihood of cross-resistance with existing insecticides [24, 25]. These results indicate that ivermectin has many attractive qualities as a novel malaria control tool, it targets mosquitoes regardless of feeding location or time. As mentioned above, we intend to formulate a mathematical model to investigate the impact of ivermectin on malaria transmission. The earliest malaria transmission model was the Ross-McDonald model [21]. On this basis, extensive research have been developed to study malaria transmission dynamics including different factors, such as seasonality, time- delay, impact of var- ious control strategies, spatial effects, stage structure of mosquitoes and humans and so on, see e.g. [2, 5, 9, 22, 30, 32] and references therein. Recently, [18] mod- elled the effect of ivermectin on malaria transmission control by ordinary differential equations and the results showed that ivermectin was significantly more effective in malaria control compared to the no-intervention state. In this paper, we formulate a novel delay malaria transmission model taking into account seasonality and iver- mectin, use the theoretical approach to analyze our model’s dynamical behavior and study the long-term effect of ivermectin on malaria transmission experiments in Kenya. We hope that our work will provide theoretical guidance for the control of malaria transmission using ivermectin in the future. The rest of this article is structured as follows. In Section 2, we derive a delay malaria transmission model with seasonality a and ivermectin, and present some properties, such as positivity and boundedness of solutions. The basic reproduction ratio R0 is discussed and use R0 to analyze the threshold dynamics of the model in Section 3. In Section 4, we demonstrate the validity of our theory by examining the long- term behavior of malaria transmission in Kenya and present an analysis of the effect of ivermectin for malaria control. Finally, in Section 5, a brief summary and discussion are given. 2. Model Formulation In this section, we propose a seasonal effect of delay malaria transmission model taking into account the treatment and ivermectin. First, we denote the total popu- lation size of humans and mosquitoes by Nh(t) and Nv(t), respectively. We classify the human population into five subclasses: susceptible Sh(t), exposed Eh(t), infec- tious Ih(t), treated Th(t) and recovered Rh(t) (those who recovered through treat- ment or recovered naturally, of which they were still slightly infectious). Mosquito EJDE-2022/28 MATHEMATICAL MODEL FOR THE TRANSMISSION OF MALARIA 3 populations are divided into two categories: susceptible Sv(t) and infectious Iv(t). Meanwhile, the time Plasmodium in completing its development in the mosquito and migrate to the salivary glands, known as the external incubation period (EIP). Let τ be the length of the EIP. Thus, Nh(t) = Sh(t) + Eh(t) + Ih(t) + Th(t) +Rh(t), Nv(t) = Sv(t) + Iv(t− τ). (2.1) If a susceptible human Sh is bitten by an infectious mosquito, then the human progresses through the exposed Eh, infectious Ih. The infected humans enter the recovery class Rh either through natural recovery or treatment. However, recov- ered humans becomes susceptible after losing immunity. Susceptible mosquitoes Sv become infected when they bite infectious, treated or recovered humans, and once infected they move into infectious class Iv. Then, the infection rates per susceptible human is cβ(t) Nv(t) Nh(t) Iv(t) Nv(t) = cβ(t) Iv(t) Nh(t) , where c is the probability of a mosquito infecting a human, β(t) represents the aver- age number of bites per mosquito at time t. By the same idea of model formulation as in [3], we see that the number of newly occurred infectious mosquitoes at time t is given by bβ(t− τ) Ih(t− τ) + σ1Th(t− τ) + σ2Rh(t− τ) Nh(t− τ) Sv(t− τ)e ∫ t t−τ µv(s)ds, where b is the probability of a human infecting a mosquito, σ1 is the ratio between the probability of transmission from a treated person to a susceptible mosquito and the probability of transmission from an infected person and σ2 is the ratio between the probability of transmission from a recovered person to a susceptible mosquito and the probability of transmission from an infected person. µv(t) represents nat- ural mortality rate of mosquitoes. Mosquitoes leave the total population through natural death and death caused by ivermectin. Mosquito mortality due to ivermectin is mainly determined by vaccination rate of ivermectin and the concentration of ivermectin [24]. Clinical trials investigating have shown that three doses of 300 µg/kg given over 3 days has a mosquitocidal effect in humans for 28 days against Anopheles gambiae sensu stricto [25]. Therefore, we suppose that the mosquito mortality due to ivermectin is a monthly periodic function, which is denoted by dv(t). Let κ be the vaccination rate of ivermectin. Then the mosquito mortality due to ivermectin can be described as dv(t)β(t) κNh(t) Nh(t) = dv(t)κβ(t). Following the above assumptions, we obtain the transmission diagram, see Figure 1: 4 Z. ZHAO, S. LI, Y. LU EJDE-2022/28 Figure 1. Transmission diagram of malaria among human and mosquito. Here M1 = cβ Iv(t) Nh(t) , M2 = bβ(t − τ) Ih(t−τ)+σ1Th(t−τ)+σ2Rh(t−τ) Nh(t−τ) Sv(t − τ)e− ∫ t t−τ µv(s)dz and M3 = dv(t)κβ(t) + µv(t). Accordingly, we obtain the following malaria transmission model with time delay dSh(t) dt = Λh + ρhRh(t)− dhSh(t)− cβ(t) Iv(t) Nh(t) Sh(t), dEh(t) dt = cβ(t) Iv(t) Nh(t) Sh(t)− dhEh(t)− νhEh(t), dIh(t) dt = νhEh(t)− dhIh(t)− δhIh(t)− γhIh(t)− αhIh(t), dTh(t) dt = αhIh(t)− dhTh(t)− ehTh(t), dRh(t) dt = γhIh(t) + ehTh(t)− dhRh(t)− ρhRh(t), dSv(t) dt = Λv(t)− bβ(t) H(t) Nh(t) Sv(t)− dv(t)κβ(t)Sv(t)− µv(t)Sv(t), dIv(t) dt = bβ(t− τ) H(t− τ) Nh(t− τ) Sv(t− τ)e− ∫ t t−τ µv(s)ds − dv(t)κβ(t)Iv(t) − µv(t)Iv(t), (2.2) where Λv(t), β(t), µ(t) are positive 12-month periodic continuous function. It is easy to see that e− ∫ t t−τ µv(s)ds is positive 12-month periodic function and H(t) = Ih(t)+σ1Th(t)+σ2Rh(t). Explicit description of the parameters of model (2.2) are given in Table 1. 2.1. Positivity and boundedness of solutions. Let C := C ( [−τ, 0],R7 ) , C+ = C ( [−τ, 0],R7 + ) . Define ‖φ‖ = ∑7 i=1 ‖φi‖∞, where ‖φi‖∞ = max−τ≤θ≤0 |φi(θ)| and φ = (φ1, φ2, φ3, φ4, φ5, φ6, φ7) ∈ C. Then, (C,C+) is an ordered Banach space , and C+ is an internally non-empty normal cone of C. For any given continuous function u : [−τ, σφ) → R7 with σφ > 0, we define ut ∈ C for t ≥ 0 by ut(θ) = u(t + θ) for all θ ∈ [−τ, 0]. EJDE-2022/28 MATHEMATICAL MODEL FOR THE TRANSMISSION OF MALARIA 5 Table 1. Biological description of model (2.2) parameters. Parameters Biological significance Λh Recruitment rate for human population ρh The probability of moving from recovered to susceptible dh Natural death rate for humans c Transmission probability of malaria from mosquitoes to susceptible humans β(t) Biting rate of mosquitoes νh The transmission rate of humans from the exposed state to the infectious state δv Malaria death rate for humans γh The probability of natural recovery of infectious humans rehabilitation αh The probability of infectious population receiving treatment eh Probability of recovery by receiving treatment Λv(t) Recruitment rate for mosquitoes population b Transmission probability of malaria from humans to susceptible mosquitoes σ1 Ratio between the probability of transmission from a treated person to a susceptible mosquito and the probability of transmission from an infected person σ2 Ratio between the probability of transmission from a recovered person to a susceptible mosquito and the probability of transmission from an infected person dv(t) Ivermectin-induced mortality of mosquitoes κ Ivermectin vaccination rate µv(t) Natural mortality rate of mosquitoes τ The time required for Plasmodium to develop and mature in mosquitoes Lemma 2.1. For any φ ∈ C+, model (2.2) has a unique non-negative solution u(t, φ) with u(0, φ) = u0 = φ such that ut(φ) ∈ C+ on t ∈ [0,∞), and all solutions are ultimately bounded. Proof. For any φ = (φ1, φ2, φ3, φ4, φ5, φ6, φ7) ∈ C+, system (2.2) can be written as u̇ = f(t, u), u(0, φ) = φ, (2.3) where the vector field f(t, u) is generated by the right side of system (2.2). Since f(t, u) is continuous in (t, u) ∈ R+×C+ and Lipschitz in u on each compact subset of R×C+, so it then follows from [11, Theorems 2.2.1 and 2.2.3] that model (2.2) has a unique solution u(t, φ) with u0 = φ on its maximum interval [0, σφ) of existence. For any φ ∈ C+ with φi(0) = 0, it is obvious that fi(t, φ) ≥ 0 for i=1,2,3,4,5,6,7. By [26, Theorem 5.2.1 and Remark 5.2.1], the unique solution u(t, φ) of model (2.2) with u0 = φ satisfies ut(φ) ∈ C+ for all t ∈ [0, σφ). 6 Z. ZHAO, S. LI, Y. LU EJDE-2022/28 From the first equation of (2.1) and model (2.2), we have dNh(t) dt = Λh − dhNh(t)− δhIh(t) ≤ Λh − dhNh(t), dSv(t) dt ≤ Λv(t)− dv(t)κβ(t)Sv(t)− µv(t)Sv(t), dIv(t) dt ≤ bβ(t− τ)Sv(t− τ)− dv(t)κβ(t)Iv(t)− µv(t)Iv(t), (2.4) when t ∈ [0, σφ). Thus, Sh(t), Eh(t), Ih(t), Th(t), Rh(t), Sv(t), Iv(t) are bounded on [0, σφ), which implies that σφ →∞ by [11, Theorem 2.3.1]. Note that the linear equation dN̂h(t) dt = Λh − dhN̂h(t) has a globally stable equi- librium point N̂∗h = Λh dh , it follows from the comparison principle that lim sup t→∞ (Sh(t) + Eh(t) + Ih(t) + Th(t) +Rh(t)) ≤ N̂∗h = Λh dh . This means that Sh(t), Eh(t), Ih(t), Th(t), Rh(t) are ultimately bounded. Similarly, we have Sv(t) is ultimately bounded by the second inequation of (2.4). Let gh = max t>τ g(t) and gl = min t>τ g(t), where g(t) is any bounded function on [τ,∞). From the seventh equation of model (2.2), we have dIv(t) dt ≤ bβhShv − ( dlvκβ l + µlv ) Iv(t), Then, Iv(t) is ultimately bounded for t > τ . This implies that all solutions of the model (2.2) are ultimately bounded. � 3. Threshold dynamics 3.1. Basic reproduction ratios. To deduce the basic reproduction ratio R0 for the model (2.2), we need to find the disease-free state of the model (2.2). Let Eh = Ih = Th = Rh = Iv = 0, we obtain dSh(t) dt = Λh − dhSh(t), dSv(t) dt = Λv(t)− dv(t)κβ(t)Sv(t)− µv(t)Sv(t). (3.1) Model (3.1) has a unique positive disease-free periodic solution E0 = (S∗h, 0, 0, 0, 0, 0, S ∗ v (t), 0), where S∗h = Λh dh , and S∗v (t) = [ ∫ t 0 Λ(r)e ∫ r 0 dv(s)κβ(s)+µv(s)dsdr + ∫ ω 0 Λ(r)e ∫ r 0 dv(s)κβ(s)+µv(s)dsdr e ∫ ω 0 dv(r)κβ(r)+µv(r)dr − 1 ] × e− ∫ t 0 dv(r)κβ(r)+µv(r)dr. EJDE-2022/28 MATHEMATICAL MODEL FOR THE TRANSMISSION OF MALARIA 7 Following the basic reproduction ratio theory developed by [34], we first linearize model (2.2) at the disease-free periodic solution E0, then we obtain dEh(t) dt = cβ(t)Iv(t)− ãEh(t), dIh(t) dt = νhEh(t)− b̃Ih(t), dTh(t) dt = αhIh(t)− c̃Th(t), dRh(t) dt = γhIh(t) + ehTh(t)− d̃Rh(t), dIv(t) dt = bβ(t− τ) H(t− τ) S∗h S∗v (t− τ)e− ∫ t t−τ µv(s)ds − g̃(t)Iv(t), (3.2) where ã = dh + νh, b̃ = dh + δh + γh + αh, c̃ = dh + eh, d̃ = dh + ρh, g̃(t) = dv(t)κβ(t) + µv(t). Let C̄ := C([−τ, 0],R5) and C̄+ := C([−τ, 0],R5 +), then (C̄, C̄+) is an ordered Banach space. Let F : R→ L(C,R5) be a map and V (t) be a 5×5 matrix function on R. For any t ∈ R and φ = (φ1, φ2, φ3, φ4, φ5)T ∈ C̄, one has F (t)φ =  cβ(t)φ5(0) 0 0 0 bβ(t− τ) Ih(−τ)+σ1Th(−τ)+σ2Rh(−τ) S∗ h S∗v (t− τ)e− ∫ t t−τ µv(s)ds  and V (t) =  ã 0 0 0 0 −νh b̃ 0 0 0 0 −αh c̃ 0 0 0 −γh −eh d̃ 0 0 0 0 0 g̃(t)  . Then the linear system (3.2) can be written as du(t) dt = (F (t)− V (t))u(t), t ≥ 0, where u(t) = (Eh(t), Ih(t), Th(t), Rh(t), Iv(t)) T . The internal evolution of individ- uals in the infectious compartments is governed by the linear ordinary differential system: du dt = −V (t)u(t). Let Φ(t, s), t ≥ s, be the evolution matrix of the above linear system. For each s ∈ R, Φ(t, s) satisfies dΦ(t, s) dt = −V (t)Φ(t, s), t ≥ s and Φ(s, s) = I, 8 Z. ZHAO, S. LI, Y. LU EJDE-2022/28 where I is the 5× 5 identity matrix. It then follows that Φ(t, s) =  e−ã(t−s) 0 0 0 0 a21e −b̃(t−s) e−b̃(t−s) 0 0 0 a31e −c̃(t−s) a32e −c̃(t−s) e−c̃(t−s) 0 0 a41e −d̃(t−s) a42e −d̃(t−s) a43e −d̃(t−s) e−d̃(t−s) 0 0 0 0 0 e− ∫ t s g̃(ξ)dξ  , where a21 = νh(t − s), a31 = 1 2αhνh(t− s)2, a32 = αh(t − s), a41 = 1 6αhνheh(t − s)3 + 1 2γhνh(t− s)2, a42 = 1 2αheh(t− s)2 + γh(t− s) and a43 = eh(t− s). Let Cω be the ordered Banach space of all ω-periodic continuous functions from R to R5, and equipped with the maximum norm ‖.‖∞ and the positive cone C+ ω := {υ ∈ Cω : υ(t) ≥ 0, for t ∈ R}. We suppose that υ(s) ∈ Cω is the initial distribution of infectious individuals in this periodic environment. Then for any given s ≥ 0, F (t − s)υt−s is the distribution of newly infected individuals at time t − s, which is produced by the infectious individuals who were introduced over the time interval [t − s − τ, t − s]. Then Φ(t, t − s)F (t − s)υt−s is the distribution of those infected individuals who were newly infected at time t−s and remain in the infected compartments at time t. It follows that∫ ∞ 0 Φ(t, t− s)F (t− s)υt−sds = ∫ ∞ 0 Φ(t, t− s)F (t− s)υ(t− s+ ·)ds is the distribution of accumulative new infections at time t produced by all those infected individuals υ(s) introduced at the previous time to t. We define a linear operator L : Cω → Cω as follows [Lυ](t) = ∫ +∞ 0 Φ(t, t− s)F (t− s)υ(t− s+ ·)ds, for any t ∈ R, υ ∈ Cω, where L is called the next infection operator. It then follows from [34] that the basic reproduction ratio R0 for system (3.2) is defined as the spectral radius of operator L, that is R0 = ρ(L). For a given t ≥ 0, let P (t) be the solution maps of system (3.2), that is P (t)φ = vt(φ) = v(t, φ), where v(t, φ) is the unique solution of system (3.2) with v0 = φ ∈ C. Then P := P (ω) is the Poincaré map associated with system (3.2). Let r(P ) be the spectral radius of P . By [34, Theorem 2.1], we have the following result. Lemma 3.1. R0 − 1 has the same sign as r(P )− 1. To study the relationship between the global properties of system (3.2) and R0, we need to prove the existence of the exponential positive solution of system (3.2). For this purpose, we first define a phase space: Y := R× [C((−τ, 0),R)]3 ×R, Y + := R× [C((−τ, 0),R+)]3 ×R. Similar to the proof of Lemma 2.1, we can obtain the following result. Lemma 3.2. For a φ ∈ Y +, system (3.2) has a unique non-negative solution v(t, φ) ∈ Y + with v0 = φ for all t ≥ 0. For a given t ≥ 0, let P̂ (t) be the solution maps of system (3.2) on Y . The fol- lowing lemma indicates that that the periodic semi-flow P̂ (t) is eventually strongly positive. EJDE-2022/28 MATHEMATICAL MODEL FOR THE TRANSMISSION OF MALARIA 9 Lemma 3.3. For φ ∈ Y + \ {0}, the solution v(t, φ) with v0 = φ of system (3.2) satisfies vi(t, φ) > 0 for all t > τ , i = 1, 2, 3, 4, 5 and hence, P̂ (t)φ � 0 for all t > τ . Proof. For a given φ = (φ1, φ2, φ3, φ4, φ5) ∈ Y + \ {0}, let v(t, φ) = (v1(t), v2(t), v3(t), v4(t), v5(t)). If v1(0) > 0, by solving system (3.2), for all t > 0, we obtain v1(t) = e−ãt ( v1(0) + ∫ t 0 cβ(ξ)v5(ξ)eãξdξ ) ≥ e−ãtv1(0) > 0, v2(t) = e−b̃t ( v2(0) + ∫ t 0 νhv1(ξ)eb̃ξdξ ) ≥ e−b̃t (∫ t 0 νhv1(ξ)eb̃ξdξ ) > 0. Similarly, we can obtain in turn v3(t) = e−c̃t ( v3(0) + ∫ t 0 αhv2(ξ)ec̃ξdξ ) > 0, v4(t) = e−d̃t ( v4(0) + ∫ t 0 (eh + γh)v3(ξ)ed̃ξdξ ) > 0, When t > τ , we find that v5(t) = e− ∫ t 0 g̃(ξ)dξv5(0) + e− ∫ t 0 g̃(ξ)dξ ∫ t 0 bβ(ξ − τ) Hh(ξ − τ) S∗h S∗v (ξ − τ)e− ∫ ξ ξ−τ µv(z)dze ∫ ξ 0 µv(z)dzdξ > 0. That implies (v1(t), v2(t), v3(t), v4(t), v5(t)) > (0, 0, 0, 0, 0) for all t > τ. (3.3) By similar discussion, we obtain (3.3) when vi(0) > 0 for i = 2, 3, 4, 5. Then we obtain P̂ (t) is strongly positive on Y + for t > τ . � Furthermore, P̂ := P̂ (ω) is the Poincaré map associated with system (3.2). Let r(P̂ ) be the spectral radius of P̂ . Choose an integer n0 > 0 such that n0ω > τ . By Lemma 3.3, we see that P̂n0 = P̂ (n0ω) is strongly positive. Following [11, Theorem 3.6.1], the linear operator P̂n0 is compact on Y +. Since r(P̂n0) = r(P̂ )n0 and according to the Krein-Rutman theorem, we obtain that r(P̂ ) is a simple eigenvalue and having a strongly positive eigenvector. By [13, Lemma 3.8], we have r(P ) = r(P̂ ). From [35], we have the following result. Lemma 3.4. Let µ = ln r(P ) ω . Then there exists a positive ω-periodic function v∗(t) such that u∗(t) = eµtv∗(t) is a positive solution of linear system (3.2). 3.2. Equilibrium stability. In this part, we establish a threshold-type result on the global dynamics of model (2.2) in terms of R0. Let X := C([−τ, 0],R7 +), X0 := {φ ∈ X : φi(0) > 0 for any i ∈ 2, 3, 4, 5, 7}, ∂X0 := X\X0 = {φ ∈ X : φ2(0)φ3(0)φ4(0)φ5(0)φ7(0) = 0}, 10 Z. ZHAO, S. LI, Y. LU EJDE-2022/28 where φ = (φ1, φ2, φ3, φ4, φ5, φ6, φ7). Theorem 3.5. If R0 > 1, then model (2.2) exists a positive ω -periodic solution, and there exist a positive real number η > 0 such that the solution (Sh(t), Eh(t), Ih(t), Th(t), Rh(t), Sv(t), Iv(t)) ∈ X with Sh(0) > 0, Eh(0) > 0, Ih(0) > 0, Th(0) > 0, Rh(0) > 0, Sv(0) > 0, Iv(0) > 0 satisfies lim inf t→∞ (Eh(t), Ih(t), Th(t), Rh(t), Iv(t)) ≥ (η, η, η, η, η). Proof. Let Q(t) : X → X be the solution maps of model (2.2), that is Q(t)φ = u(t, φ) = ut(φ), t ≥ 0, where u(t, φ) is the unique solution of the model (2.2) with u0 = φ ∈ X. Then Q := Q(ω) is the Poincaré map associated with model (2.2), and Qn = Q(nω). It is easy to find that Q(t) (X0) ⊂ X0 for all t ≥ 0. Now, we prove that Q(t) is uniformly persistent with respect to (X0, ∂X0). Let M1 = (S∗h, 0, 0, 0, 0, S ∗ v0, 0), where S∗v0(θ) = S∗v (θ) for θ ∈ [−τ, 0]. Then Q(t)M1 = M1 for all t ≥ 0. Since limφ→M1 ‖Q(t)φ−Q(t)M1‖ = 0 uniformly for t ∈ [0, ω], for any given ε > 0, there exists a positive real number δ = δ(ε) such that for any φ satisfying ‖φ−M1‖ < δ, we have ‖Q(t)φ−Q(t)M1‖ ≤ ε, for t ∈ [0, ω]. We proceed with the following two claims. Claim 1. lim supn→∞ ‖Qn(φ)−M1‖ ≥ δ for all φ ∈ X0. Suppose by contradiction that lim supn→∞ ‖Qn(ψ)−M1‖ < δ for some ψ ∈ X0. Then, there exists an integer N ≥ 1 such that‖Qn(ψ) − M1‖ < δ for all n ≥ N . For any t ≥ Nω, letting t = nω + t′ with n ≥ N and t′ ∈ [0, ω), we have ‖Q(t)ψ −Q(t)M1‖ = ‖Q(t′)(Qn(ψ))−Q(t′)M1‖ < ε. Then for all t ≥ Nω, S∗h − ε < Sh(t, ψ) < S∗h + ε, S∗v (t)− ε < Sv(t, ψ) < S∗v (t) + ε, 0 < Eh(t, ψ), Ih(t, ψ), Th(t, ψ), Rh(t, ψ), Iv(t, ψ) < ε. We also have the inequalities Sh(t, ψ) Nh(t, ψ) ≥ S∗h − ε S∗h + 5ε = 1− 6ε S∗h + 5ε and Sv(t, ψ) Nh(t, ψ) ≥ S∗v (t)− ε S∗h + 5ε . Then from model (2.2), for t ≥ Nω + τ , we have dEh(t) dt ≥ cβ(t) ( 1− 6ε S∗h + 5ε ) Iv(t)− (dh + νh)Eh(t), dIh(t) dt = νhEh(t)− (dh + δh + γh + αh) Ih(t), dTh(t) dt = αhIh(t)− (dh + eh)Th(t), dRh(t) dt = γhIh(t) + ehTh(t)− (dh −m)Rh(t), dIv(t) dt ≥ bβ(t− τ) H(t− τ) S∗h + 5ε (S∗v (t− τ)− ε)e− ∫ t t−τ µv(s)dz − (dv(t)κβ(t) + µv(t))Iv(t). (3.4) EJDE-2022/28 MATHEMATICAL MODEL FOR THE TRANSMISSION OF MALARIA 11 Let Pε be the solution map of the perturbed linear periodic system dĒh(t) dt = cβ(t) ( 1− 6ε S∗h + 5ε ) Īv(t)− (dh + νh)Ēh(t), dĪh(t) dt = νhĒh(t)− (dh + δh + γh + αh)Īh(t), dT̄h(t) dt = αhĪh(t)− (dh + eh)T̄h(t), dR̄h(t) dt = γhĪh(t) + ehT̄h(t)− (dh +m)R̄h(t), dĪv(t) dt = bβ(t− τ) H̄(t− τ) S∗h + 5ε (S∗v (t− τ)− ε)e− ∫ t t−τ µv(s)dz − (dv(t)κβ(t) + µv(t))Īv(t), with Pε := Pε(ω). Since R0 > 1, limε→0+ r(Pε) = r(P ) > 1, there exists a suf- ficiently small ε > 0 such that ε < min{mint∈[0,ω] S ∗ v (t), S∗h} and r(Pε) > 1. By Lemma 3.4, there is a positive ω-periodic function v∗ε (t) such that v̄ε(t) = eλtv∗ε (t) is a solution of above perturbed linear periodic system, where λ = ln r(Pε) ω > 0. Then we have v̄ε(t) → +∞ as t → +∞. And because of the system (3.4), the comparison principle implies that lim t→∞ (Eh(t, ψ), Ih(t, ψ), Th(t, ψ), Rh(t, ψ), Ih(t, ψ)) = (∞,∞,∞,∞,∞,∞), which leads to a contradiction. By claim 1 above, we can see that M1 is an isolated invariant set for Q in X, and W s(M1) ∩X0 = ∅, where W s(M1) is the stable set of M1 for Q. Claim 2. M∂ = {φ ∈ ∂X0 : φi(0) = 0, i = 2, 3, 4, 5, 7}, where M∂ := {φ ∈ ∂X0 : Qn(φ) ∈ ∂X0, for n ≥ 0}. Clearly, it suffices to prove that for any ϕ ∈M∂ , (Eh(t, ϕ), Ih(t, ϕ), Th(t, ϕ), Rh(t, ϕ), Ih(t, ϕ)) = (0, 0, 0, 0, 0) holds for all t ≥ 0. Suppose not, then there exists some t0 ≥ 0 such that Eh(t0, ϕ) > 0 or Ih(t0, ϕ) > 0 or Th(t0, ϕ) > 0 or Rh(t0, ϕ) > 0 or Iv(t0, ϕ) > 0. Assuming that Eh(t0, ϕ) > 0, then by the second equation of the model (2.2) we know that dEh dt ≥ −(dh + νh)Eh(t) and Eh(t, ϕ) ≥ Eh(t0, ϕ)e(dh+νh)(t0−t) > 0, for t > t0. Following the third equation of the model (2.2) we know that Ih(t, ϕ) > ∫ t t0 νhEh(s, ϕ)e(s−t)(dh+δh+γh+αh)ds > 0, for all t > t0. By the same argument, Th(t, ϕ) > 0 and Rh(t, ϕ) > 0 for any t > t0. Similar, it is easy to see that Iv(t, ϕ) > ∫ t t0 bβ(s− τ) H(s− τ, ϕ) Nh(s− τ) Sv(s− τ, ϕ)e ∫ s s−τ µ(z)dze ∫ s t (dvκβ(z)+µ(z))dzds > 0, when t > t0 + τ , where Nh(s − τ) = Sh(s − τ, ϕ) + Eh(s − τ, ϕ) + Ih(s − τ, ϕ) + Th(s − τ, ϕ) + Rh(s − τ, ϕ). Thus, (Eh(t, ϕ), Ih(t, ϕ), Th(t, ϕ), Rh(t, ϕ), Ih(t, ϕ)) > (0, 0, 0, 0, 0) for all t > t0 + τ . Then we can find some n > 0 with nω > t0 + τ such that (Eh(nω, ϕ), Ih(nω, ϕ), Th(nω, ϕ), Rh(nω, ϕ), Ih(nω, ϕ)) /∈M∂ 12 Z. ZHAO, S. LI, Y. LU EJDE-2022/28 for all t > t0 + τ , which is a contradiction. Similarly, when Ih(0) > 0 or Th(0) > 0 or Rh(0) > 0 or Iv(0) > 0, we can obtain the same contradiction. Therefore, Claim 2 is proved. Moreover, from system (3.1), we have Sh(t, φ)→ S∗h, Sv(t, φ)→ S∗v (t) as t→∞, i.e., Qn(φ)→M1. Thus, cupψ∈M∂ ω(ψ) = {M1}. That means M1 cannot form a cycle in ∂X0. By Claims 1 and 2 and the acyclicity theorem on uniform persistence for maps [35, Theorem 1.3.1 and Remark 1.3.1], it follows that Q : X → X is uniformly persistent with respect to (X0, ∂X0). We can prove the practical uniform persistence, that is, there exists an η > 0 such that lim inf t→∞ min(Eh(t, φ), Ih(t, φ), Th(t, φ), Rh(t, φ), Iv(t, φ)) = lim inf t→∞ p(Q(t)φ) ≥ η, for all φ ∈ X0. The proof is quite standard, a more detailed display of a similar reasoning can be found in [12]. � Theorem 3.6. If R0 < 1 and δh = 0, the disease-free periodic solution E0 is globally attractive for model (2.2) in X. Proof. In view of model (2.2) and (3.1), there exists a sufficiently large integer n > 0 with nω ≥ τ and a such that Nh(t) ≥ N∗h − ε = Λh dh − ε, Sv(t) ≤ S̄∗v (t) + ε, for all t ≥ nω − τ. Therefore, for all t ≥ nω, we have dEh(t) dt ≤ cβ(t)Iv(t)− (dh + νh)Eh(t), dIh(t) dt = νhEh(t)− (dh + γh + αh)Ih(t), dTh(t) dt = αhIh(t)− (dh + eh)Th(t), dRh(t) dt = γhIh(t) + ehTh(t)− (dh +m)Rh(t), dIv(t) dt ≤ bβ(t− τ) H(t− τ) N∗h − ε (S∗v (t− τ) + ε)e− ∫ t t−τ µv(s)ds − (dv(t)κβ(t) + µv(t))Iv(t). EJDE-2022/28 MATHEMATICAL MODEL FOR THE TRANSMISSION OF MALARIA 13 Let Pε be the Poincaré map of the auxiliary system dẼh(t) dt = cβ(t)Ĩv(t)− (dh + νh)Ẽh(t), dĨh(t) dt = νhẼh(t)− (dh + γh + αh)Ĩh(t), dT̃h(t) dt = αhĨh(t)− (dh + eh)T̃h(t), dR̃h(t) dt = γhĨh(t) + ehT̃h(t)− (dh +m)R̃h(t), dĨv(t) dt = bβ(t− τ) H̃(t− τ) N∗h − ε (S∗v (t− τ) + ε)e− ∫ t t−τ µv(s)ds − (dv(t)κβ(t) + µv(t))Ĩv(t). (3.5) Since limε→0 r(Pε) = r(P ) < 1, we fix a sufficiently small ε ∈ (0, N∗h) and r(Pε) < 1. In a similar manner, there is a positive ω-periodic function v∗ε (t) such that v̄ε(t) = eλtv∗ε (t) is a positive solution of (3.5), where λ = ln r(Pε) ω < 0. Assuming that E∗0 = (Eh(t, ψ), Ih(t, ψ), Th(t, ψ), Rh(t, ψ), Ih(t, ψ)) is a positive ω-periodic solution of system (3.2). we can choose a sufficiently large constant K > 0 such that E∗0 ≤ Kv̄ε(t) for all t ∈ [nω, nω+ τ ]. It follows from [26, Theorem 5.1.1], we can get E∗0 ≤ Kv̄ε(t), for all t ≥ nω + τ . Thus, we have lim t→∞ E∗ = (Eh(t), Ih(t), Th(t), Rh(t), Iv(t)) = (0, 0, 0, 0, 0) when R0 < 1. Meanwhile, we can obtain the following limit system for model (2.2): dSh(t) dt = Λh − dhSh(t), dSv(t) dt = Λv(t)− dv(t)κβ(t)Sv(t)− µv(t)Sv(t). Because S∗h is globally asymptotically stable and S∗v (t) is globally attractive, we have lim t→∞ ( Sh(t), Eh(t), Ih(t), Th(t), Rh(t), Sv(t), Iv(t) ) = ( S∗h, 0, 0, 0, 0, 0, S ∗ v (t), 0 ) by the chain transitive sets arguments [35, Theorem 1.2.1]. � 4. Study case In this section, based on malaria data and clinical trial data of ivermectin in malaria transmission control from Kenya, we estimate the values of parameters related to malaria and invermectin-induced mortality of mosquitoes dv(t), Then, the numerical fitted curve of malaria transmission cases is shown in Figure 3 and 4. Sensitivity analysis of the ivermectin vaccination rate κ and mosquito bite rate β(t), and analysis of the effect of ivermectin vaccination rate, number of and timing between vaccine intervention rounds on malaria transmission control are given in Figure 7 and Figure 11, respectively. 4.1. Parameter estimation. According to the information provided by World Health Organization [28], the total population of Kenya in 2013 is Nh = 45519986 14 Z. ZHAO, S. LI, Y. LU EJDE-2022/28 and the average life expectancy of humans is 63.419 years. Thus, the human natural death rate dh can be calculated as follows : dh = 1 63.419× 12 ≈ 0.0013 Month−1. The recruitment rate Λh is: Λh = dh × 45519986 ≈ 59813 Humans×Month−1. The values of constant parameters for system (3.2) that do not heavily depend on temperature are listed in Table 2 from [8, 14, 23, 27]. Next, we evaluate the periodic parameters µ(t) and β(t) in system (3.2) by using the monthly mean temperature data from 2002 to 2020 from [29], which is shown in Table 3. Table 2. Parameter values. parameters Values Dimension Sources Λh 59175 Humans See text ρh 0.01672 Month−1 [8, 14] dh 0.0013 Month−1 See text c 0.01 Dimensionless [8] β(t) To be evaluated Month−1 See text νh 3.344 Month−1 [8] δh 0.002736 Month−1 [8] αh 0.53 Dimensionless [23] eh 0.85 Dimensionless [27] γh 0.13984 Month−1 [8] Λv(t) To be evaluated Month−1 See text b 0.2 Dimensionless [8, 14] σ1 0.5 Dimensionless Assumed σ2 0.1 Dimensionless Assumed µv(t) To be evaluated Month−1 See text dv(t) To be evaluated Dimensionless Assumed κ [0,1] Dimensionless Assumed τ 9/30.4 Dimensionless [8] Table 3. Monthly mean temperature of Kenya (in ◦C). Month Jan. Feb. Mar. Apr. May Jun. Temperature 25.50 26.28 26.74 26.11 25.05 24.06 Month Jul. Aug. Sep. Oct. Nov. Dec. Temperature 23.44 23.83 24.59 25.37 25.08 25.03 It follows from [14] that the temperature dependent mosquito biting rate can be expressed as: β(C) = 30.4 107.204− 13.3523C + 0.677509C2 − 0.0159732C3 + 0.000144876C4 Month−1 (4.1) EJDE-2022/28 MATHEMATICAL MODEL FOR THE TRANSMISSION OF MALARIA 15 where C represents temperature in ◦C. Substituting the temperatures in Table 3 into (4.1), the biting rate of mosquitoes can be approximated by β(t) = 7.989 + 0.4487 cos(πt/6) + 0.7684 sin(πt/6)− 0.4603 cos(2πt/6) − 0.2658 sin(2πt/6)− 0.09456 cos(3πt/6)− 0.03793 sin(3πt/6) + 0.01472 cos(4πt/6) + 0.02339 sin(4πt/6) + 0.004283 cos(5πt/6) + 0.05418 sin(5πt/6) Month−1. Similarly, the temperature-dependent death rate of adult mosquitoes is given by [14], µv(C) = 30.4 + 29.564e(−C−278◦K 2.7035 ) Month−1. (4.2) Then the death rate of adult mosquitoes can be fitted as µv(t) = 3.058− 0.004607 cos(πt/6)− 0.005952 sin(πt/6) + 0.003278 cos(2πt/6) + 0.003186 sin(2πt/6) + 0.0008082 cos(3πt/6)− 0.0008577 sin(3πt/6) − 0.0001089 cos(4πt/6)− 000007509 sin(4πt/6) + 0.0001428 cos(5πt/6) − 0.0003749 sin(5πt/6) Month−1. To estimate the maturation function of the mosquito, we suppose that the egg deposition rate is a linear function of the biting rate [14], Λv(t) = 5× β(t)×Nh Mosquitoes×Month−1. In 2018, [25] studied safety and mosquitocidal efficacy of high-dose ivermectin in Kenya adults. They found ivermectin 300 µg/kg per day for 3 days provided a good balance between efficacy and tolerability, and reduced mosquito survival for at least 28 days after treatment. It can be seen that this drug shows promise as a potential new tool for malaria elimination. By fitting the data of mosquito survival post treatment in [25], we obtain the ivermectin-induced mortality of mosquitoes on post-treatment days 0, 2 + 4h, 7, 10, 14, 21 and 28, see Table 4. Table 4. Ivermectin-induced mortality of mosquitoes. Ivermectin test times (Day) 2+4h 7 10 14 21 28 Drug-induced mortality 0.7817 0.6782 0.4364 0.3156 0.1479 0.1105 And we can get the following fitted function for the mortality rate of mosquito (see Figure 2) dv(t) = 0.9668e−0.07611t. 0 5 10 15 20 25 30 Time (Day) 0 0.2 0.4 0.6 0.8 1 T he m or ta lit y ra te o f m os qu ito Figure 2. Ivermectin-induced mortality of mosquitoes. 16 Z. ZHAO, S. LI, Y. LU EJDE-2022/28 4.2. Model validation. From [17], We can get data on the monthly malaria cases generated in Kenya from January 2013 to December 2020. From Figure 3, the monthly reported numbers of malaria has pronounced seasonality in Kenya. Based on the estimated parameter values above and the initial values: Sh(0) = 34850000, Eh(0) = 13960, Ih(0) = 359300, Th(0) = 517200, Rh(0) = 12889000, Sv(0) = 247100000, Iv(0) = 1236000, we fit the Kenyan malaria cases by model (2.2). The reported data, the third-order Fourier fitted function for these data and the simu- lation result in Kenya from January 2013 to December 2015 are shown in Figure 3. Jan Apr Jul Oct Jan Apr Jul Oct Jan Apr Jul Oct Month 1 2 3 4 5 6 7 8 M al ar ia C as es ×105 Monthly Reported Cases Monthly Reported Case Fit Curves Model prediction curve Figure 3. Comparison between the reported malaria cases from 2013 to 2015 and the simulation cases form model (2.2). The numbers of malaria in the future several years in Kenya is shown in Figure 4 with no further effective control measure is taken and the initial values: Sh(0) = 37930000, Eh(0) = 53770, Ih(0) = 415200, Th(0) = 817200, Rh(0) = 8820000, Sv(0) = 246300000, Iv(0) = 1597000. (4.3) 2017 2018 2019 2020 2021 2022 2023 Time (Year) 1 2 3 4 5 6 7 8 M al ar ia C as es ×105 Monthly Reported Cases Monthly Reported Case Fit Curves Model prediction curve Figure 4. Malaria development trend by forecasting model (2.2). 4.3. Long-term behavior. In this subsection, we verify of the theoretical results by computing the basic reproduction R0 and simulating the long-term behavior of the model (2.2) under the same set of parameter values as Figure 3 and initial values (4.3) Assuming that vaccination rate κ = 0.2, we can obtain R0 = 1.0691 > 1, and the images of Ih and Iv are shown in Figure 5. In this case, the disease persists and eventually shows stable periodical fluctuations. If vaccination rate is increased to κ = 0.7, then R0 = 0.7487 < 1. In this case, the long-term behavior of infectious EJDE-2022/28 MATHEMATICAL MODEL FOR THE TRANSMISSION OF MALARIA 17 mosquitoes and humans are shown in Figure 6, which implies that malaria will eventually die out. These simulations are consistent with the results of Theorem 3.5 and Theorem 3.6. 0 100 200 300 400 500 Time (Month) 0 1 2 3 4 In fe ct io us h um an s ×105 0 100 200 300 400 500 Time (Month) 0 0.5 1 1.5 2 2.5 3 3.5 In fe ct io us m os qu ito es ×106 (a) (b) Figure 5. Long-term behaviors of the infectious compartments in model (2.2) when R0 = 1.0691. 0 100 200 300 400 500 Time (Month) 0 1 2 3 4 5 In fe ct io us h um an s ×105 0 100 200 300 400 500 Time (Month) 0 0.5 1 1.5 2 In fe ct io us m os qu ito es ×106 (a) (b) Figure 6. Long-term behaviors of the infectious compartments in model (2.2) when R0 = 0.7487. 4.4. Sensitivity analysis of R0. To explore the effectiveness of ivermectin vaccine in the control of malaria, it is important to analyze the relationship between some parameters of model (2.2) and R0. We mainly consider the effect of ivermectin vaccination rate κ and mosquito bite rate β(t) on R0. We use the parameter values in Table 2 and the initial values (4.3). Firstly, we discuss the effect of ivermectin vaccination rate κ on R0. By keeping the other parameter values the same as those in Table 2, we observe that R0 is a decreasing function of κ and R0 < 1 when κ > 0.2837 (see Figure 7a). Therefore, increasing the ivermectin vaccination rate can be an effective way to control malaria transmission. Then, to simulate the effect of the bite rate β(t) on R0, we make κ = 0 and replace β(t) with β̂(t) = (1− q)β(t), where q can be considered as the efficiency of people’s reduced mosquito bites. Under other parameter values the same as those in Table 2, we observe that R0 is a decreasing function of q and R0 < 1 when q > 0.1570 (see Figure 7b). That means malaria transmission gradually decreases 18 Z. ZHAO, S. LI, Y. LU EJDE-2022/28 and will eventually disappear completely by controlling the rate of mosquito bites. Figure 8 gives the graph of R0 as a function of the parameters κ and q. It can be seen that increasing the values of κ and q is helpful for reducing R0. 0 0.2 0.4 0.6 0.8 1 κ 0.6 0.8 1 1.2 1.4 R 0 0 0.2 0.4 0.6 0.8 1 q 0 0.2 0.4 0.6 0.8 1 1.2 1.4 R 0 (a) (b) Figure 7. Sensitivity analysis of R0. a Relationship between R0 and κ. b Relationship between R0 and q. 0 0.2 0 0 0.4 0.6 0.1 0.1 0.8 R 0 1 0.2 0.2 1.2 1.4 0.3 0.30.4 0.4 q κ 0.5 0.50.6 0.60.7 0.70.8 0.80.9 0.91 1 Figure 8. Influence of κ and q on R0. 4.5. Analysis of the effect of ivermectin. In this section, we explore the effect of ivermectin on malaria transmission in various regimens of use by model (2.2). We mainly discuss the effect of ivermectin vaccination rate, number of and timing between vaccine intervention rounds on malaria transmission control. Figure 9 shows the effect of different ivermectin vaccination rate κ on Ih and Iv by 1 month apart of vaccine intervention. We find that the higher value of κ is, the lower level Ih and Iv can be reduced to, and when κ = 1, malaria disappears after only about 300 months. Since malaria transmission is seasonal in some areas, adjusting the frequency of administration in line with the malaria transmission season could have an effect on malaria transmission [24]. In addition, seasonal vaccine intervention is easier to administer and also effect in reducing malaria outbreaks. We consider four ivermectin regimens: 3 rounds given 1 month apart, 3 rounds given 2 month apart, 3 rounds given 3 month apart and 4 rounds given 1 month apart (three consecutive daily doses of 300 µg/kg per day) and ivermectin vaccination rate κ = 0.7. As shown in Figure 10, ivermectin has a significant effect on malaria transmission season, the effect of continuous vaccines interventions is greater than that of intermittent EJDE-2022/28 MATHEMATICAL MODEL FOR THE TRANSMISSION OF MALARIA 19 0 100 200 300 400 500 Time (Month) 0 1 2 3 4 5 6 7 In fe ct io us h um an s ×105 κ = 0 κ = 0.2 κ = 0.6 κ = 1 0 100 200 300 400 500 Time (Month) 0 2 4 6 8 10 In fe ct io us m os qu ito es ×106 κ = 0 κ = 0.2 κ = 0.6 κ = 1 Figure 9. Different κ correspond to different infection humans. vaccines interventions and increases in the vaccines interventions frequency could reduce malaria transmission season. As Kenya is a seasonal malaria transmission country, Figure 10 also shows that the interventions will still have an impact in year 2 and the impact trend is similar to Figure 3.8 in the Appendix of [24] (with different doses of the drug). 0 5 10 15 20 Time (Month) 0 2 4 6 8 10 In fe ct io us h um an s ×105 No Ivermectin Three rounds of ivermectin 1 month apart Three rounds of ivermectin 2 month apart Three rounds of ivermectin 3 month apart Four rounds of ivermectin 1 month apart Continuous use of ivermectin 1 month apart Figure 10. Impact of different ivermectin regimens on malaria control. Using the four ivermectin regimens above, we simulate the long-term behaviors of the infectious humans in Figure 11. We observe that the higher frequency of vaccines interventions is, the lower level Ih(t) can be reduced to. 0 50 100 150 200 250 300 350 400 450 500 Time (Month) 0 2 4 6 8 10 12 In fe ct io us h um an s ×105 No Ivermectin Three rounds of ivermectin 1 month apart Three rounds of ivermectin 2 month apart Three rounds of ivermectin 3 month apart Four rounds of ivermectin 1 month apart Continuous use of ivermectin 1 month apart Figure 11. Long-term impact of different ivermectin regimens on malaria control. 20 Z. ZHAO, S. LI, Y. LU EJDE-2022/28 5. Discussion Ivermectin is the only drug in the avermectin class of endectocides that is avail- able for human use. It can reduce lifespan of mosquitoes, and also causes sec- ondary behavioural and reproductive disturbances that could affect mosquito sur- vival. Based on the experiment results obtained [6, 7, 25], we developed a delay malaria transmission model incorporating seasonality and ivermectin. The basic reproduction ratio R0 is derived by the theory developed in [34]. By appealing to the theory of persistence of dynamical systems and the theory of chain transitive sets, we obtained R0 is the threshold parameter for the extinction and persistence of malaria. That is, if R0 < 1, then infective compartments approach zero even- tually; if R0 > 1, then malaria will persist. Numerically, we have estimated all constant and periodic parameters from some published data and studied malaria transmission in Kenya. We verified our theoretical results by simulating the long- term behavior of the solution. By showing a graph of how R0 varies with ivermectin vaccination rate and bite rate, we found that it is possible to eliminate malaria from Kenya when we combined ivermectin with tools to control bite rates. Furthermore, it takes at least 25 years to eliminate malaria from Kenya with malaria control measures intact (see Figure 9). We also simulated four ivermectin regimens and found that the higher the dosing frequency is, the lower level Ih(t) can be decreased to. Furthermore, our simulation of the effect in year 2 is similar to [24]. Thus, the model is more realistic for the control of malaria transmission. From the above analysis, we found that ivermectin may be effective in a malaria transmission seasonal. At present, there are many ivermectin-related clinical trials being conducted, for example: one in Guinea-Bissau, one in Thailand and a mul- tisite study in Mozambique and Tanzania [24]. We expected to study the effect of ivermectin by constructing a reasonable mathematical model with these relevant clinical and entomological data, and combining with different transmission settings, different malaria vectors, and different control measures in the future. Acknowledgments. This research was supported by grants FRF-TP-20-068A1 from the Fundamental Research Funds for the Central Universities. We are grateful to the anonymous referees for their valuable comments. References [1] G. J. Abiodun, et al.; Modelling the influence of temperature and rainfall on the population dynamics of Anopheles arabiensis, Malar J., 15 (2016), 364–379. [2] F. B. Agusto, N. Marcus, K. O. Okosun; Application of optimal control to the epidemiology of malaria, Electron. J. Differ. Eq., 81 (2012), 1–22. [3] I. Al-Darabsah, Y. Yuan; A periodic disease transmission model with asymptomatic carriage and latency periods, J. Math. Biol., 77 (2018), 343–376. [4] S. Bhatt, D. Weiss, E. Cameron, et al.; The effect of malaria control on Plasmodium falci- parum in Africa between 2000 and 2015, Nature, 526 (2015), 207–211. [5] C. Boëte, F. B. Agusto, R. G. Reeves; Impact of mating behaviour on the success of malaria control through a single inundative release of transgenic mosquitoes, J. Theoret. Biol., 347 (2014), 33–43. [6] C. J. Chaccour, et al.; Ivermectin to reduce malaria transmission: a research agenda for a promising new tool for elimination, Malar J., 12 (2013), 153–161. [7] C. J. Chaccour, et al.; Establishment of the Ivermectin Research for Malaria Elimination Network: updating the research agenda, Malar J., 14 (2015), 243–251. EJDE-2022/28 MATHEMATICAL MODEL FOR THE TRANSMISSION OF MALARIA 21 [8] 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 (2008), 1272–1296. [9] B. Dembele, A. -A. Yakubu; Optimal treated mosquito bed nets and insecticides for eradi- cation of malaria in Missira, Discrete Contin. Dyn. Syst. Ser. B, 17 (2017), 1831–1840. [10] D. A. Ewing, C. A. Cobbold, B. V. Purse, M. A. Nunn, S. M. White; Modelling the effect of temperature on the seasonal population dynamics of temperate mosquitoes, J. Theoret. Biol., 400 (2016), 65–79. [11] J. K. Hale, S. M. Verduyn Lunel; Introduction to functional-differential equations, Appl. Math. Sci., 99, 1993. [12] F. Li, J. Liu, X. Q. Zhao; A West Nile virus model with vertical transmission and periodic time delays, J. Nonlinear Sci., 30 (2020), 449–486. [13] Y. Lou, X. -Q. Zhao; A theoretical approach to understanding population dynamics with seasonal devel opmental durations, J. Nonlinear Sci., 27 (2017), 573–603. [14] Y. Lou, X. -Q. Zhao; A climate-based malaria transmission model with structured vector population, SIAM J. Appl. Math., 6 (2010), 2023–2044. [15] G. Macdonald; The epidemiology and control of malaria, Oxford University Press, London, UK, 1957. [16] Mectizan Donation Program. Annual highlights 2020. [17] Ministry Of Health Data Analytics Platform (MOH DAP). [18] G.A. Ngwa, et al.; On a three-stage structured model for the dynamics of malaria transmis- sion with human treatment, adult vector demographics and one aquatic stage, J. Theoret. Biol., 481 (2019), 202–222. [19] K. Okuneye, A. B. Gumel; Analysis of a temperature- and rainfall-dependent model for malaria transmission dynamics, Math. Biosci., 287 (2017), 72–92. [20] R. N. Price, L. Von Seidlein, N. Valecha, F. Nosten, J. K. Baird, N. J. White; Global extent of chloroquine-resistant Plasmodium vivax: a systematic review and meta-analysis, Lancet Infect. Dis., 14 (2014), 982–991. [21] R. Ross; The Prevention of Malaria, John Murray, London, 1911. [22] S. G. Ruan, D. M. Xiao, J. C. Beier; On the delayed Ross-Macdonald model for malaria transmission. Bull. Math. Biol., 70 (2008), 1098–1114. [23] I. Saran, E. M. Maffioli, D. Menya, et al.; Household beliefs about malaria testing and treatment in Western Kenya: the role of health worker adherence to malaria test results, Malar J., 16 (2017), 349. [24] H. C. Slater, B. D. Foy, K. Kobylinski, et al.; Ivermectin as a novel complementary malaria control tool to reduce incidence and prevalence: a modelling study, Lancet Infect. Dis., 20 (2020), 498–508. [25] M. R. Smit, E. O. Ochomo, G. Aljayyoussi, et al.; Safety and mosquitocidal efficacy of high- dose ivermectin when co-administered with dihydroartemisinin-piperaquine in Kenyan adults with uncomplicated malaria (IVERMAL): a randomised, double-blind, placebo-controlled trial, Lancet Infect. Dis., 18 (2018), 615–639. [26] H. L. Smith; Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems, American Mathematical Society, Providence, RI, 1995. [27] P. O. Sumba, S. L. Wong, H. K. Kanzaria, et al.; Malaria treatment-seeking behaviour and recovery from malaria in a highland area of Kenya, Malar. J., 7 (2008), 245. [28] The World Bank, Kenya, [29] The World Bank, Climate Change Knowledge Portal, https://climateknowledgeportal.worldbank.org/country/Kenya/climate-data-historical [30] B. Tsanou, J. C. Kamgang, J. M. S. Lubuma, et al.; Modeling pyrethroids repellency and its role on the bifurcation for a bed net malaria model, Chaos Soliton Fractals, 136 (2020), 109809. [31] H. Wang, G. Tian; Propagating interface in reaction-diffusion equations with distributed delay, Electron. J. Differential Equations, 2021 (2021), no. 54, 1-22. [32] W. A. Woldegerima, R. Ouifki, J. Banasiak; Mathematical analysis of the impact of transmission-blocking drugs on the population dynamics of malaria, Appl. Math. Comput., 400 (2021), 126005. [33] World Health Organisation (WHO). https://www.who.int/health-topics/malaria#tab=tab 1 22 Z. ZHAO, S. LI, Y. LU EJDE-2022/28 [34] X. Q. Zhao; Basic reproduction ratios for periodic compartmental models with time delay. J. Dynam. Differential Equations , 29 (2015), 67–82. [35] X. Q. Zhao; Dynamical Systems in Population Biology, Second edition, CMS Books in Math- ematics, Springer, Cham, 2017. Zhihong Zhao School of Mathematics and Physics, University of Science and Technology Beijing, Beijing 100083, China Email address: zzh@ustb.edu.cn Shaochun Li School of Mathematics and Physics, University of Science and Technology Beijing, Beijing 100083, China Email address: leeschun@foxmail.com Yulan Lu School of Mathematics and Physics, University of Science and Technology Beijing, Beijing 100083, PR China Email address: luyulan@ustb.edu.cn 1. Introduction 2. Model Formulation 2.1. Positivity and boundedness of solutions 3. Threshold dynamics 3.1. Basic reproduction ratios 3.2. Equilibrium stability 4. Study case 4.1. Parameter estimation 4.2. Model validation 4.3. Long-term behavior 4.4. Sensitivity analysis of R0 4.5. Analysis of the effect of ivermectin 5. Discussion Acknowledgments References