Mathematical model of malaria transmission dynamics with distributed delay and a wide class of nonlinear incidence rates

Abstract

Generally, the infection process of most vector-borne diseases involves a latent period in both human hosts and vectors. With regards to other publications, Tian and Song have recently proposed an SIR-SI model to analyze the effects of the incubation period on a vector-borne disease with nonlinear transmission rate. But they were silent on the fact that the partially immune individuals are slightly infective to mosquitoes. So, by considering that the partially immune individuals remain slightly infective to mosquitoes, a similar work has been done in this paper for malaria global transmission dynamics following an SIRS-SI pattern. The basic reproduction ratio has been calculated using the next-generation matrix method. Furthermore, using the characteristic equations and inequality analytical techniques, conditions are given under which the system exhibits threshold behavior as follows: when R₀ < 1, the disease-free equilibrium is globally asymptotically stable meaning that the disease will eventually die out; and the unique endemic equilibrium is globally asymptotically stable when R₀ > 1 meaning that the disease will persist. Finally, some numerical simulations have been performed to illustrate our theoretical results.

Keywords:

• malaria transmission dynamics
• delay differential equations
• basic reproduction number
• wide class of incidence rates
• numerical simulations
• limit cycles

JEL classification:

• 37N30
• 37N25
• 65L12
• 65U05

PUBLIC INTEREST STATEMENT

Malaria is one of the reemerging diseases which still remains a major public health issue in several low-income countries around the world. Although mathematical models are an abstract simplification of reality, they can be used to describe or predict the outcomes of epidemics or pandemics providing information that is crucial in informing public health policymakers. This allows them to optimize the use of their limited resources. Otherwise, the last few decades have witnessed an increase of importance of applying equations with distributed delay for understanding real-world phenomena. Therefore, by using delay differential equations, we provide in this paper a further understanding of the impact of incubation period and the role of partially immune individuals in the malaria transmission dynamics and its lasting negative effects. The paper covers a large information on mathematical models of malaria and provides interesting reading due to the original approach and the rich contents.

1. Introduction

Malaria is a vector-borne disease that has been a major public health issue. It is caused by different species of the Plasmodium parasite; the most common being Plasmodium falciparum. The disease is widespread in the tropical and subtropical regions that exist in a broad band around the equator. This includes much of Sub-Saharan Africa, Asia, and Latin America. In 2016, there were 216 million cases of malaria worldwide resulting in an estimated 731 thousand deaths. Approximately 90% of both cases and deaths occurred in Africa. Rates of the disease have decreased from 2000 to 2015 by 37%, but increased from 2014 during which there were 198 million cases. Malaria is commonly associated with poverty and has a major negative effect on economic development. In Africa, it is estimated to result in losses of $12$ billion dollars a year due to increased health-care costs, lost ability to work, and negative effects on tourism (Beretta & Kuang, 2002; Blyuss & Kyrychko, 2014; Traoré, Sangaré, & Traoré, 2018).

Mathematical models can project how infectious diseases progress to show the likely outcome of an epidemic and help inform public health interventions. Models use some basic assumptions and mathematics to find parameters for various infectious diseases and use those parameters to calculate the effects of possible interventions, like mass vaccination or drugs distribution programs. The standard technique for developing mathematical descriptions of infectious diseases is to model the system as a set of autonomous ordinary differential Equations (Birkhoff & Rota, 1982; Cheng, Wang, & Yang, 2012; Jordan & Smith, 2004; Ouedraogo, Ouedraogo, & Sangaré, 2018). This is an immensely powerful approach which has led to many insights into the factors that affect prevalence and control most of the emerging infectious diseases as dengue, malaria, chikungunya, yellow fever, zika virus, etc. (Diekmann, Getto, & Gyllenberg, 2008; Kermack & McKermack, 1927; Koutou, Traoré, & Sangaré, 2018; Olaniyi & Obabiyi, 2013).

The first mathematical model for understanding malaria transmission has been developed by Ross (Koutou et al., 2018; Moulay, Aziz-Alaoui, & Cadivel, 2011). The Ross’s model consists of two nonlinear differential equations in two state variables that correspond to the proportions of infected human beings and the infected mosquitoes. Macdonald (1957) improved Ross’s more famous differential equations model with some biological assumptions and entomological field data. Since then, multiple models have been done for malaria by several authors. For instance, Aron and May (Chitnis, Hyman, & Cushing, 2008; Khan, Khan, & Islam, 2015; Traoré et al., 2018) added various features of malaria to the model of Macdonald, such as an incubation period in the mosquito, super-infection and a period of immunity in human beings. Besides, the inclusion of acquired immunity proposed by (Chitnis, Cushing, & Hyman, 2006; Ruan, Xiao, & Beier, 2008) was a major point of malaria modeling. (Ngwa & Shu, 2000) proposed a model for malaria transmission in which human hosts follow an SEIRS-like pattern, and vector hosts follow the SEI pattern due to their short life cycle. Moreover, Chitnis et al.(2008) extended the model of Ngwa and Shu by assuming that although individuals in the recovered class are immune, in the sense that they do not suffer from serious illness and do not contract clinical malaria, they still have low levels of Plasmodium in their bloodstream and can infect the susceptible mosquitoes (Bai, 2015; Macdonald, 1957; Ngwa & Shu, 2000; Traoré, Sangaré, & Traoré, 2017).

However, many other aspects have aroused much attention and interest in recent decades such as the climate effects on the dynamics of the vector population and the biting rate from mosquitoes to humans (Khan et al., 2015; Zhang, Jia, & Song, 2014), the stage structure of the vector (Koutou et al., 2018; Moulay et al., 2011; Traoré et al., 2017), the incubation period in the hosts (Diekmann et al., 2008; Khan et al., 2015) and the form of incidence function. For instance, Traoré et al. (2018) and Koutou et al. (2018) have presented a nonautonomous model and an autonomous model for malaria transmission including the immature stages of the mosquitoes, respectively. Indeed, during the propagation of the epidemic, time delays exist because an individual may not be infectious until some time after becoming infected (Beretta & Kuang, 2002; Zhang et al., 2014). Some time is required before the infective organism develops in the vector to the level that is sufficient to pass the infection further (Khan et al., 2015; Posny & Wang, 2014; Van Den Driessche, 1996). When we consider the influence of delay on disease spreading, the models take a form of delay differential Equations (Diekmann & Gyllenberg, 2012; Kuang, 1993; Saker, 2010). There are many forms of delays among which we have the concentrated delay, the discrete delay, and the distributed delay. Generally, time delay is assumed to be a single-value. As it is well-known, a constant delay may be applied if the movement time is precisely known, which is not very realistic for a physical situation. In nature, it is obvious that the process depending on the past does not depend on the exact time $t-\tau$. But the dependence is somehow distributed. Besides, the length of the latent period varies from disease to disease. For example, according to Anderson and May (Khan et al., 2015), measles virus has an infectious period of $5-7$ days, and influenza has only $1-3$ days latent period, besides Acquired Immune Deficiency Syndrome (AIDS) has an incubation time ranging from a few months to years after the patient has been shown to get antibodies to the human immunodeficiency virus (HIV) (Omondi, Mbogo, & Luboobi, 2018; Wang, Cao, Alsaedi, & Ahmad, 2017). Accordingly, equations with distributed delay provide a more flexible and realistic description for real-world phenomena.

In the existing literatures, there exist several works about delay differential equations for infectious diseases especially malaria (Li, Shuai, & Wang, 2010; Nakata, Enatsu, & Muroya, 2011; Yang & Xiao, 2010), to name a few. In (Zhang et al., 2014), considering the intracellular delay between initial infection of a cell biting by an infected female anopheles mosquito and the release of new virus particles, Zhang et al. (2014) have proposed and analyzed a non-autonomous SI-SIRS model of malaria global transmission dynamics with distributed time delay. Otherwise, Tian and Song (2017) have more recently constructed and studied a model of a vector-borne disease with two delays and nonlinear transmission rate as follows:

(1) $\left\{\begin{array}{c}\begin{array}{c}\frac{S_h(t)}{dt}=\mathcal{A}-f_1(S_h(t),I_v(t)-{\mu }_h{s}_h(t)\\ \frac{I_h(t)}{dt}={\int }_0^{\tau h}gh(\tau )f_1(s_h(t-\tau ),I_v(t-\tau ))d\tau -({\mu }_h+{\alpha }_h)I_h(t)\end{array}\\ \frac{R_h(t)}{dt}={\alpha }_h{I}_h(t)-{\mu }_v{R}_v(t)\\ \frac{S_v(t)}{dt}=\mathcal{B}-f_1(s_v(t),I_h(t))-{\mu }_v{S}_v(t)\\ \frac{I_v(t)}{dt}={\int }_0^{\tau v}gv(\tau )f_2(S_v(t-\tau ),I_h(t-\tau ))d\tau -{\mu }_v{I}_v(t)\end{array}\right.$(1)

where $\tau \ge 0$, $max\left\{{\tau }_h,{\tau }_v\right\}>0$, $t\ge 0$ and $S_v$, $I_v$, $S_h$, $I_h$, $R_h$ denote, respectively, the susceptible mosquitoes, the infectious mosquitoes, the susceptible humans, the infectious humans, and the partially immune individuals.

Indeed, malaria differs from many other vector-borne diseases because for it, the partially immune individuals remain slightly infectious to the susceptible mosquitoes (Chiyaka, Tchuenche, Garira, & Dube, 2008). Even though they do not show symptoms of malaria, they are still reservoirs of parasites. However, due to their size in an endemic area, their role in the disease transmission cannot be neglected. But both of the previous studies firstly proposed by Zhang et al. (2014) and recently by Tian and Song (2017), respectively, failed to mention this aspect. Motivated by all these works and more precisely by system (1), we have proposed an autonomous SIRS-SI epidemic model to describe malaria transmission dynamics with a wide class of incidence rates and distributed delay. We have then extended the study proposed in (Tian & Song, 2017) to malaria including a specific biological assumption that partially immune individuals are infective to susceptible mosquitoes.

The remaining part of the paper is organized as follows. Section 2 concerns the model formulation. In this Section, we have also presented some useful preliminaries and assumptions for the mathematical analysis of the model. In Section 3, we have derived some results about the basic reproduction number, the disease-free equilibrium, and the endemic equilibrium. The global stability of the disease-free equilibrium and the endemic equilibrium were, respectively, studied in Section 4 and 5. Section 6 is devoted to the numerical simulations. Section 7 concludes the paper with some remarks.

2. Model formulation and preliminaries

2.1. Model formulation

In this part, we formulate a model of malaria spread dynamics relying on some biological aspects. The first compartmental model for malaria transmission has been proposed by Ross. Since many other types of models have been constructed by several authors. Among these types of models, we have the SIRS-SI one that is used in this study. In that model, the human hosts are supposed to follow SIRS and the vector hosts follow SI. The flow between the different compartments is shown by the chart of Figure 1. For humans, the three compartments represent the total population of humans at a given time $t\ge 0$, who are susceptible $S_H$, infectious $I_H$ and partially immune $R_H$. At the infectious stage, a host may die from the disease or may recover with acquired partial immunity. Partially immune hosts still have protective antibodies and other immune effectors at low levels (Chiyaka et al., 2008; Dietz, Molineaux, & Thomas, 1974; Macdonald, 1957). For mosquitoes, the two compartments represent susceptible mosquitoes $S_m$ and infectious mosquitoes $I_m$ at a given time $t\ge 0$. Infectious mosquitoes are those that have sporozoites in their salivary glands and can infect the human hosts. Due to the mosquito’s short life cycle, it cannot be recovered from the infection.

Figure 1. The dashed arrows indicate the direction of the infection and the solid arrows represent the transition from one class to another.

From the above descriptions, the total sizes of human population $N_H(t)$ and vector population $N_m(t)$ are given by

$N_H(t)=S_H(t)+I_H(t)+R_H(t)$

and

$N_m(t)=S_m(t)+I_m(t)$

respectively.

Functions $f$ and $g$ denote the incidence functions and the constant values $b$, $\alpha$, $d$, ${\mathrm{\Lambda }}_H$, ${\mathrm{\Lambda }}_m$ $\beta$, $\delta$, $\gamma$, $\mu$ and $\sigma$ are the model parameters and can be described as follows:

• $b$ is the natural death rate of the mosquito population,

• $\alpha$ is the rate of contact from humans to mosquitoes when susceptible mosquitoes get into contact with the infective humans and the number of new infections is given as follows:

(2) $\bar{f}=\alpha {\int }_0^h\eta (\tau )f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))d\tau$(2)

• $d$ is the disease-induced death rate of the mosquito population,

• ${\mathrm{\Lambda }}_H$ is considered as the immigration rate of human population and new individuals which are born into the susceptible state with a constant birth rate,

• ${\mathrm{\Lambda }}_m$ is the constant recruitment rate into the susceptible class. It includes the rate of the immature mosquitoes that become adults,

• $\beta$ is the contact rate from mosquito to human when susceptible humans make contact with infected mosquitoes and we write the number of new infections as follows:

(3) $\bar{g}=\beta {\int }_0^h\eta (\tau )g(S_H(t),I_m(t-\tau ))d\tau$(3)

• $\delta$ is the constant natural death rate of human hosts,

• $\gamma$ is the waning immunity rate of human population. In fact, it is assumed that a recovered immune individual can lose their immunity, so that a member of the state $R_H$ returns to the susceptible state with the constant rate $\gamma \ge 0$,

• $\mu$ is the transfer rate from infectious to immune class. It is assumed that the immune individuals are slightly infectious to susceptible mosquitoes,

• $\sigma$ is the disease-induced death rate of the human population.

The disease transmission dynamic is described through the flowchart in Figure 1. The mathematical model is the following system of delay differential equations

(4) $\left\{\begin{array}{c}\frac{d{S}_m(t)}{dt}={\mathrm{\Lambda }}_m-\alpha {\int }_0^h\eta (\tau )f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))d\tau -b{S}_m(t)\\ \hspace{0.17em}\\ \frac{d{I}_m(t)}{dt}=\alpha {\int }_0^h\eta (\tau )f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))d\tau -(d+b)I_m(t)\\ \hspace{0.17em}\\ \frac{d{S}_H(t)}{dt}={\mathrm{\Lambda }}_H-\beta {\int }_0^h\eta (\tau )g(S_H(t),I_m(t-\tau ))d\tau -\delta S_H(t)+\gamma R_H(t)\\ \hspace{0.17em}\\ \frac{d{I}_H(t)}{dt}=\beta {\int }_0^h\eta (\tau )g(S_H(t),I_m(t-\tau ))d\tau -(\delta +\mu +\sigma )I_H(t)\\ \hspace{0.17em}\\ \frac{d{R}_H(t)}{dt}=\mu I_H(t)-(\delta +\gamma )R_H(t).\end{array}\right.$(4)

The constant $h$ represents the superior limit of incubation times. The incubation period distribution is $\eta (\tau )$, and denotes the fraction of the hosts in which the time taken to become infectious is $\tau$. It is assumed to be continuous on $[0,h]$ satisfying ${\int }_0^h\eta (\tau )d\tau =1$ and $\eta (\tau )\ge 0$ for $0\le \tau \le h$.

Remark 2.1. The two incidence functions $f$ and $g$ explain the contact between two species: humans and mosquitoes (Diekmann & Gyllenberg, 2012).

2.2. Preliminaries

In this section, we introduce some initial conditions for system (4) as follows:

For any $\theta \in [-h,0]$,

(5) $S_m(\theta )={\phi }_1(\theta ),I_m(\theta )={\phi }_2(\theta ),S_H(\theta )={\phi }_3(\theta ),I_H(\theta )={\phi }_4(\theta ),R_H(\theta )={\phi }_5(\theta )$(5)

where $\phi =({\phi }_1,{\phi }_2,{\phi }_3,{\phi }_4,{\phi }_5{)}^T\in \mathcal{C}$ such as ${\phi }_i(\theta )\ge 0,$ $(i=1,...,5)$ for all $\theta \in [-h,0]$.

$\mathcal{C}$ denotes the Banach space $C([-h,0],{\mathbb{R}}_{+})$ of continuous functions mapping the interval $[-h,0]$ into ${\mathbb{R}}_{+}$ is equipped with the sup norm: $\left|\left|\phi \right|\right|=\underset{-h\le \theta \le 0}{sup}\left|\phi (\theta )\right|$.

For a biological meaning Zhang et al. (2014) we suppose that, ${\phi }_i(\theta )>0$ for all $i=1,...,5$, and we consider system (4) in the phase space ${\mathcal{C}}^5$. It is clear that all solutions of (4) in ${\mathcal{C}}^5$ with initial conditions given by (5), remain non-negative.

In the special case, when only discrete experimental data of the past have influence on the present rate of change of the state, the corresponding mathematical equations can be either retarded or neutral differential-difference equations (RDDE or NDDE, respectively). However, the implementation of the complex-step method and the MatLab solver DDE23 are appropriate to compute sensitivities for delay differential equations that have discontinuities in their history functions. This provides a very good approximation of sensitivities as long as discontinuities in the initial data do not cause loss of smoothness in the solution (Bellman & Cooke, 1963; Tian & Song, 2017).

The incidence rates $f$ and $g$ are non-negative differentiable functions on the non-negative quadrant. Furthermore, $f$ and $g$ are positive if and only if both arguments are positive. Note also that when there is no one infected in the human population and the mosquito population, then the incidence functions are equal to zero.

To make biological sense, we assume that the incidence functions $f$ and $g$ satisfy the following conditions for all $t\ge 0$ and $\tau \in [0,h]$:

• (H1) $f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))>0$ and $g(S_H(t),I_m(t-\tau ))>0$

• (H2) $f(S_m(t),0,0)=0$ and $g(S_H(t),0)=0$

• (H3) $\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H(t)}(S_m(t),I_H(t-\tau ),R_H(t-\tau ))>0$, $\frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H(t)}(S_m(t),I_H(t-\tau ),R_H(t-\tau ))>0$ and $\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m(t)}(S_H(t),$ $I_m(t-\tau ))>0$.

The conditions in assumption (H1) state that the infection occurs when there is contact between infected mosquitoes and susceptible humans and also between susceptible mosquitoes and infected or partially immune humans. The assumption (H2) ensures that the model has a unique disease-free equilibrium; and (H3) means that the rate of new infection increases both with not only the infected population size, but also with partially immune individuals (Posny & Wang, 2014).

Theorem 2.1. (Bony, 1969) Let $L:{\mathbb{R}}^n\to \mathbb{R}$ be a differentiable function and let $X(x)$ be a vector field. Let $a\in \mathbb{R}$ and let $G(x)=\left\{x\in {\mathbb{R}}^n:L(x)\le a\right\}$ be such that $\mathrm{\nabla }L(x)\ne 0$ for all $x\in L^{-1}(a)=\left\{x\in {\mathbb{R}}^n:L(x)=a\right\}$. If $⟨X(x)|\mathrm{\nabla }L(x)⟩\le 0$ for all $x\in L^{-1}(a)$ then, the set $G$ is positively invariant.

Lemma 2.1. For any initial conditions given in (5), system (4) admits a unique positive solution $(S_m(t),I_m(t),S_H(t),I_H(t),R_H(t))$ for $t\ge 0$. Furthermore,

$\underset{t\to +\mathrm{\infty }}{lim}(S_m(t)+I_m(t))\le \frac{{\mathrm{\Lambda }}_m}{b}and\underset{t\to +\mathrm{\infty }}{lim}(S_H(t)+I_H(t)+R_H(t))\le \frac{{\mathrm{\Lambda }}_H}{\delta }$

Proof. Note that the right hand side of (4) is completely continuous and locally Lipschitzian on $\mathcal{C}$. So, from (Hale, 1977; Kuang, 1993) we deduce that the solution $(S_m(t),I_m(t),S_H(t),I_H(t),R_H(t))$ of (4) exists and is unique on $[0,t_{max})$, where $t_{max}$ represents the maximal time.

Let $x=(S_m,I_m,S_H,I_H,R_H)$. We will prove that the set $\left\{S_m\ge 0\right\}$ is positively invariant. Then, let $L(x)=-S_m$. $L$ is differentiable and $\mathrm{\nabla }L(x)=(-1,0,0,0,0)\ne 0_{{\mathbb{R}}^5}$ for all $x\in L^{-1}(x)=$ $\left\{x\in {\mathbb{R}}^5:L(x)=0\right\}$.

The vector field on the set $\left\{S_m=0\right\}$ is given by

$X(x)=\left(\begin{array}{cccc}{\mathrm{\Lambda }}_m& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ -(d+b)I_m& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ {\mathrm{\Lambda }}_H-\beta {\int }_0^h\eta (\tau )g(S_H(t),I_m(t-\tau ))d\tau -\delta S_H(t)+\gamma R_H(t)& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \beta {\int }_0^h\eta (\tau )g(S_H(t),I_m(t-\tau ))d\tau -(\delta +\mu +\sigma )I_H(t)& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \mu I_H(t)-(\delta +\gamma )R_H(t)& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\end{array}\right)$

Then, $⟨X(x)|\mathrm{\nabla }L(x)⟩=-{\mathrm{\Lambda }}_m<0$; this proves that the set $\left\{S_m\ge 0\right\}$ is positively invariant. Similarly, we prove that the sets $\left\{I_m\ge 0\right\}$, $\left\{S_H\ge 0\right\}$, $\left\{I_H\ge 0\right\}$, $\left\{R_H\ge 0\right\}$ are positively invariant. So, $\left\{(S_H,I_H,R_H,S_m,I_m)\in {\mathbb{R}}^5:S_H\ge 0,I_H\ge 0,R_H\ge 0,S_m\ge 0,I_m\ge 0\right\}$ is positively invariant.

However, adding the first two equations and adding the last three equations of system (4), we, respectively, obtain

(6) $\frac{d(S_m(t)+I_m(t))}{dt}={\mathrm{\Lambda }}_m-b(S_m(t)+I_m(t))-d{I}_m(t)$(6)

(7) $\frac{d(S_H(t)+I_H(t)+R_H(t))}{dt}={\mathrm{\Lambda }}_H-\delta (S_H(t)+I_H(t)+R_H(t))-\sigma I_H(t).$(7)

Since $I_m(t),I_H(t)\ge 0$, it follows that

(8) $\frac{d(S_m(t)+I_m(t))}{dt}\le {\mathrm{\Lambda }}_m-b(S_m(t)+I_m(t))$(8)

(9) $\frac{d(S_H(t)+I_H(t)+R_H(t))}{dt}\le {\mathrm{\Lambda }}_H-\delta (S_H(t)+I_H(t)+R_H(t))$(9)

Hence, using the differential inequality theorem given in (Birkhoff & Rota, 1982), it follows that

(10) $S_m(t)+I_m(t)\le \frac{{\mathrm{\Lambda }}_m}{b}+\left(S_m(0)+I_m(0)-\frac{{\mathrm{\Lambda }}_m}{b}\right)exp(-bt)$(10)

(11) $S_H(t)+I_H(t)+R_H(t)\le \frac{{\mathrm{\Lambda }}_H}{\delta }+\left(S_H(0)+I_H(0)+R_H(0)-\frac{{\mathrm{\Lambda }}_H}{\delta }\right)exp(-\delta t)$(11)

Thus,

$\underset{t\to +\mathrm{\infty }}{lim}(S_m(t)+I_m(t))\le \frac{{\mathrm{\Lambda }}_m}{b}$

and

$\underset{t\to +\mathrm{\infty }}{lim}(S_H(t)+I_H(t)+R_H(t))\le \frac{{\mathrm{\Lambda }}_H}{\delta }.$

□

Furthermore, model (4) is mathematically and biologically well-defined in the following set

$\mathrm{\Delta }=\left\{\left(S_m(t),I_m(t),S_H(t),I_H(t),R_H(t)\right)\in {\mathcal{C}}^5|\begin{array}{c}0\le S_m(t)+I_m(t)\le \frac{{\mathrm{\Lambda }}_m}{b}\\ 0\le S_H(t)+I_H(t)+R_H(t)\le \frac{{\mathrm{\Lambda }}_H}{\delta }\end{array}\right\}.$

Thanks to Lemma 2.1, system (4) has the following limit system for any $\tau \ge 0$, $h>0$, $t\ge 0$, (Nakata et al., 2011)

(12) $\left\{\begin{array}{c}\frac{d{S}_m(t)}{dt}={\mathrm{\Lambda }}_m-\alpha {\int }_0^h\eta (\tau )f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))d\tau -b{S}_m(t)\\ \hspace{0.17em}\\ \frac{d{I}_m(t)}{dt}=\alpha {\int }_0^h\eta (\tau )f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))d\tau -(d+b)I_m(t)\\ \hspace{0.17em}\\ \frac{d{S}_H(t)}{dt}={\mathrm{\Lambda }}_H-\beta {\int }_0^h\eta (\tau )g(S_H(t),I_m(t-\tau ))d\tau -\delta S_H(t)+\gamma R_H(t)\\ \hspace{0.17em}\\ \frac{d{I}_H(t)}{dt}=\beta {\int }_0^h\eta (\tau )g(S_H(t),I_m(t-\tau ))d\tau -(\delta +\mu +\sigma )I_H(t)\\ \hspace{0.17em}\\ \frac{d{R}_H(t)}{dt}=\mu \left(\frac{{\mathrm{\Lambda }}_H}{\delta }-S_H(t)-R_H(t)\right)-(\delta +\gamma )R_H(t)\end{array}\right.$(12)

The limit system (12) (see Theorem 3.1 in (Beretta & Takeuchi, 1995)) is derived by using the relations

$I_H(t)=N_H(t)-S_H(t)-R_H(t)$

and

$\underset{t\to +\mathrm{\infty }}{lim}{N}_H(t)=\frac{{\mathrm{\Lambda }}_H}{\delta }.$

3. The model equilibria and the basic reproduction number

3.1. The disease-free equilibrium

Theorem 3.1. Model (4) admits two disease-free equilibria $e_0$ and $E_0$ respectively given by

$e_0=(0,0,0,\frac{{\mathrm{\Lambda }}_H}{\delta },0)$

and

(13) $E_0=\left(\frac{{\mathrm{\Lambda }}_m}{b},0,0,\frac{{\mathrm{\Lambda }}_H}{\delta },0\right).$(13)

Proof. On the one hand when we set all the infected classes to zero, the following equations yield $S_m^0=\frac{{\mathrm{\Lambda }}_m}{b}$ and $S_H^0=\frac{{\mathrm{\Lambda }}_H}{\delta }$. We then obtain the second infection-free equilibrium point $E_0$.

□

On the other hand, setting the susceptible mosquitoes class to zero, we deduce the first infection-free equilibrium point $e_0$.

Remark 3.1. The disease-free equilibrium $e_0$ represents the case where there is no disease and the area is completely devoid of mosquitoes.

The disease-free equilibrium $E_0$ represents the case where the mosquito population exists but there is no malaria.

The first case is difficult to obtain in most areas where malaria is intensive because of the difficulties and the environmental effects due to the complete eradication of the mosquito population. Hence, in this paper, we focus our study on the last one $E_0$ because it is more biologically realistic (Traoré et al., 2018).

3.2. The basic reproduction number

The basic reproduction number is now unanimously recognized as a key concept in mathematical epidemiology. It is commonly denoted by ${\mathcal{R}}_0$ and represents the average number of new infections produced directly and indirectly by a single infective individual, when that individual is introduced into a completely susceptible population during the infective period. For about 20 years, ${\mathcal{R}}_0$ has been part of the majority of research articles using mathematical modeling. For classical epidemic models, it is common that the basic reproduction number is a threshold in the sense that, when it is less than unity, on average each infected individual infects fewer than one individual, and the disease dies out. If ${\mathcal{R}}_0>1$, on average, each infected individual infects more than one other individual, so we would expect the disease to spread (Van den Driessche, 1996; Van den Driessche & Watmough, 2002).

In fact, new infections occur during contacts between susceptible individuals and infectious individuals. Here, only the following three classes $I_m$, $I_H$ and $R_H$ are considered to be disease states. Therefore, using the technique described in (Van den Driessche & Watmough, 2002), we compute the basic reproduction number as follows.

Linearizing the infected equations of system (4) at disease-free equilibrium $E_0$, we obtain the following system:

(14) $\left[\begin{array}{c}{\dot{I}}_m\\ {\dot{I}}_H\\ {\dot{R}}_H\end{array}\right]=\left(F-V\right)\left[\begin{array}{c}{I}_m\\ I_H\\ R_H\end{array}\right]$(14)

where the matrices $F$ and $V$ are given by

(15) $F=\left(\begin{array}{ccc}0& \alpha \frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A)& \alpha \frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}(A)\\ \beta \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B)& 0& 0\\ 0& 0& 0\end{array}\right)$(15)

and

(16) $V=\left(\begin{array}{ccc}(d+b)& 0& 0\\ 0& (\delta +\mu +\sigma )& 0\\ 0& -\mu & (\delta +\gamma )\end{array}\right).$(16)

The matrix $F{V}^{-1}$ is called the next-generation operator, and it is known that the basic reproduction is ${\mathcal{R}}_0=\rho (F{V}^{-1})$. For system (4), it is given by

(17) $\begin{array}{rl}& {\mathcal{R}}_0^2=\frac{\alpha \beta }{(\delta +\mu +\sigma )(d+b)}\times \frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}\left(\frac{{\mathrm{\Lambda }}_m}{b},0,0\right)\times \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}\left(\frac{{\mathrm{\Lambda }}_H}{\delta },0\right)\\ & +\frac{\alpha \beta \mu }{(\delta +\mu +\sigma )(d+b)(\delta +\gamma )}\times \frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}\left(\frac{{\mathrm{\Lambda }}_m}{b},0,0\right)\times \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}\left(\frac{{\mathrm{\Lambda }}_H}{\delta },0\right).\end{array}$(17)

From above, it is obvious that the basic reproduction number varies directly to the number of partially immune individuals who remain infective to mosquito population.

3.3. The endemic equilibrium

Note that the delay systems and the ordinary differential equation systems share the same equilibria (Beretta & Kuang, 2002; Yang & Xiao, 2010).

Lemma 3.1. If ${\mathcal{R}}_0>1$, there exists a unique positive solution $(S_m^{\ast },I_H^{\ast },R_H^{\ast },S_H^{\ast },I_m^{\ast })$ of model (4) with initial conditions (5) satisfying to the following system

(18) $\left\{\begin{array}{c}{\mathrm{\Lambda }}_m-\alpha f(S_m^{\ast },I_H^{\ast },R_H^{\ast })-b{S}_m^{\ast }=0\\ \hspace{0.17em}\\ \alpha f(S_m^{\ast },I_H^{\ast },R_H^{\ast })-(d+b)I_m^{\ast }=0\\ \hspace{0.17em}\\ {\mathrm{\Lambda }}_H-\beta g(S_H^{\ast },I_m^{\ast })-\delta S_H^{\ast }+\gamma R_H^{\ast }=0\\ \hspace{0.17em}\\ \beta g(S_H^{\ast },I_m^{\ast })-(\delta +\mu +\sigma )I_H^{\ast }=0\\ \hspace{0.17em}\\ \mu I_H^{\ast }-(\delta +\gamma )R_H^{\ast }=0.\end{array}\right.$(18)

Proof. At a fixed point $(S_m,I_m,S_H,I_H,R_H)$ e the incidence functions and the constant valu

(19) $\left\{\begin{array}{c}{\mathrm{\Lambda }}_m-\alpha f(S_m,I_H,R_H)-b{S}_m=0\\ \hspace{0.17em}\\ \alpha f(S_m,I_H,R_H)-(d+b)I_m=0\\ \hspace{0.17em}\\ {\mathrm{\Lambda }}_H-\beta g(S_H,I_m)-\delta S_H+\gamma R_H=0\\ \hspace{0.17em}\\ \beta g(S_H,I_m)-(\delta +\mu +\sigma )I_H=0\\ \hspace{0.17em}\\ \mu I_H-(\delta +\gamma )R_H=0.\end{array}\right.$(19)

From the last equation of the system we obtain

$R_H^{\ast }=\frac{\mu }{\delta +\gamma }{I}_H^{\ast }.$

Adding the first two equations of (19), we obtain that

$S_m^{\ast }=\frac{{\mathrm{\Lambda }}_m-(d+b)I_m^{\ast }}{b}.$

Moreover, adding the third and the fourth equations of (19), it then follows that

$S_H^{\ast }=\frac{{\mathrm{\Lambda }}_H-(\delta +\mu +\sigma -\frac{\mu \gamma }{\delta +\gamma })I_H^{\ast }}{\delta }.$

Let us consider the following function

(20) $H(I_m^{\ast },I_H^{\ast })=\left(H_1(I_m^{\ast },I_H^{\ast }),H_2(I_m^{\ast },I_H^{\ast })\right)$(20)

with

(21) $\begin{array}{rl}& H_1(I_m^{\ast },I_H^{\ast })=\frac{\alpha \beta (\delta +\gamma )}{\delta +\gamma +\mu \gamma }\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}\left(\frac{{\mathrm{\Lambda }}_H}{\delta },0\right)f\left(\frac{{\mathrm{\Lambda }}_m-(d+b)I_m^{\ast }}{b},I_H^{\ast },\frac{\mu \gamma }{\delta +\gamma }{I}_H^{\ast }\right)\\ & -\frac{1}{2}(\delta +\mu +\sigma )(d+b)I_H^{\ast }\end{array}$(21)

and

(22) $\begin{array}{rl}& H_2(I_m^{\ast },I_H^{\ast })=\alpha \beta \mu \frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}\left(\frac{{\mathrm{\Lambda }}_m}{b},0,0\right)g\left(\frac{{\mathrm{\Lambda }}_H-(\delta +\mu +\sigma -\frac{\mu \gamma }{\delta +\gamma })I_H^{\ast }}{\delta },I_m^{\ast }\right)\\ & -\frac{1}{2}(\delta +\mu +\sigma )(d+b)(\delta +\gamma )I_m^{\ast }.\end{array}$(22)

We can see that the function $H$ is continuous and $H(0,0)=0$.

Any zero $(I_m^{\ast },I_H^{\ast })$ of $H$ is in the set

$\left[0,\frac{{\mathrm{\Lambda }}_m}{d+b}\right]\times \left[0,\frac{{\mathrm{\Lambda }}_H}{\delta +\mu +\sigma -\frac{\mu \gamma }{\delta +\gamma }}\right]$

where $I_m^{\ast }$ and $I_H^{\ast }$ correspond to endemic equilibrium components.

Moreover,

$H\left(\frac{{\mathrm{\Lambda }}_m}{d+b},\frac{{\mathrm{\Lambda }}_H}{\delta +\mu +\sigma -\frac{\mu \gamma }{\delta +\gamma }}\right)\le 0$

and then, in order to have a solution to the equation

$H(I_m^{\ast },I_H^{\ast })=0$

on the set

$\left[0,\frac{{\mathrm{\Lambda }}_m}{d+b}\right]\times \left[0,\frac{{\mathrm{\Lambda }}_H}{\delta +\mu +\sigma -\frac{\mu \gamma }{\delta +\gamma }}\right],$

the function $H$ must be increasing at $(0,0)$.

Positing $Z=(I_m^{\ast },I_H^{\ast })$ leads to the following condition

$\frac{dH}{dZ}(0,0)>0.$

It follows that

(23) $\alpha \beta \frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}\left(\frac{{\mathrm{\Lambda }}_m}{b},0,0\right)\times \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}\left(\frac{{\mathrm{\Lambda }}_H}{\delta },0\right)-\frac{1}{2}(\delta +\mu +\sigma )(d+b)>0,$(23)

and

(24) $\alpha \beta \mu \frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}\left(\frac{{\mathrm{\Lambda }}_m}{b},0,0\right)\times \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}\left(\frac{{\mathrm{\Lambda }}_H}{\delta },0\right)-\frac{1}{2}(\delta +\mu +\sigma )(d+b)(\delta +\gamma )>0.$ (24)

Hence,

(25) $\frac{\alpha \beta }{(\delta +\mu +\sigma )(d+b)}\times \frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}\left(\frac{{\mathrm{\Lambda }}_m}{b},0,0\right)\times \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}\left(\frac{{\mathrm{\Lambda }}_H}{\delta },0\right)>\frac{1}{2},$(25)

and

(26) $\frac{\alpha \beta \mu }{(\delta +\mu +\sigma )(d+b)(\delta +\gamma )}\times \frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}\left(\frac{{\mathrm{\Lambda }}_m}{b},0,0\right)\times \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}\left(\frac{{\mathrm{\Lambda }}_H}{\delta },0\right)>\frac{1}{2}.$(26)

Adding the inequalities (25) and (26), we observe that the condition $\frac{dH}{dZ}(0,0)>0$ is equivalent to ${\mathcal{R}}_0>1$. □

4. The local stability of equilibria

In this section, we study the local stability of our model equilibria.

Theorem 4.1. If ${\mathcal{R}}_0<1$, and the assumption (H2) holds, then the disease-free equilibrium $E_0$ is locally asymptotically stable.

Proof. We begin by linearizing Equation (4) at $E_0=(A,B)$, with $A=\left(\frac{{\mathrm{\Lambda }}_m}{b},0,0\right)$ and $B=\left(\frac{{\mathrm{\Lambda }}_H}{\delta },0\right)$.

From hypothesis (H2), we have for all positive $S_m(t)$,

$f(S_m(t),0,0)=0$

and then

$\frac{\mathrm{\partial }f}{\mathrm{\partial }{S}_m(t)}\left(\frac{{\mathrm{\Lambda }}_m}{b},0,0\right)=0.$

Moreover, for all positive $S_H(t)$, the following condition yields

$g(S_H(t),0)=0,$

and consequently,

$\frac{\mathrm{\partial }g}{\mathrm{\partial }{S}_H(t)}\left(\frac{{\mathrm{\Lambda }}_H}{\delta },0\right)=0.$

The linearization of system (4) at the disease-free equilibrium is given by

(27) $\left\{\begin{array}{c}\frac{d{s}_m(t)}{dt}=-b{s}_m(t)-\alpha {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H(t)}(A)i_H(t-\tau )d\tau -\alpha {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H(t)}(A)r_H(t-\tau )d\tau \\ \hspace{0.17em}\\ \frac{d{i}_m(t)}{dt}=\alpha {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H(t)}(A)i_H(t-\tau )d\tau +\alpha {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H(t)}(A)r_H(t-\tau )d\tau -(d+b)i_m(t)\\ \hspace{0.17em}\\ \frac{d{s}_H(t)}{dt}=-\delta s_H(t)+\gamma r_H(t)-\beta {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m(t)}(B)i_m(t-\tau )d\tau \\ \hspace{0.17em}\\ \frac{d{i}_H(t)}{dt}=\beta {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m(t)}(B)i_m(t-\tau )d\tau -(\delta +\mu +\sigma )i_H(t)\\ \hspace{0.17em}\\ \frac{d{r}_H(t)}{dt}=\mu i_H(t)-(\delta +\gamma )r_H(t)\end{array}\right.$ (27)

Here, we apply the principle of linearized stability method established in (Diekmann et al., 2008; Diekmann & Gyllenberg, 2012).

Looking as in (Blyuss & Kyrychko, 2014; McCluskey, 2010), the solution of the following form,

(28) $(s_m(t),i_m(t),s_H(t),i_H(t),r_H(t))=e^{\lambda t}(s_m(0),i_m(0),s_H(0),i_H(0),r_H(0))$(28)

where $s_m(0),i_m(0),s_H(0),i_H(0),r_H(0)\in {\mathbb{R}}_{+}^{\ast }$ and $\lambda \in \mathbb{C}$.

Therefore, plugging the relation (28) into (27) and canceling $e^{\lambda t}$ from each term, yields the following matrix equation:

(29) $0=A_0\left[\begin{array}{c}{i}_m(0)\\ i_H(0)\\ r_H(0)\end{array}\right]$(29)

with the next-generation matrix $A_0$ given by

(30) $A_0=\left(\begin{array}{cccc}-(\lambda +d+b)& \alpha {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A)e^{-\lambda \tau }d\tau & \alpha {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}(A)e^{-\lambda \tau }d\tau & \hspace{0.17em}\\ \beta {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B)e^{-\lambda \tau }d\tau & -(\lambda +\delta +\mu +\sigma )& 0& \hspace{0.17em}\\ 0& \mu & -(\lambda +\delta +\gamma )& \hspace{0.17em}\end{array}\right).$(30)

There exist non-zero solutions if and only if $det(A_0)=0$. Thus, we obtain the characteristic equation as follows:

$-(\lambda +d+b)(\lambda +\delta +\mu +\sigma )(\lambda +\delta +\gamma )$

(31) $+\alpha \beta (\lambda +\delta +\gamma ){\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A)e^{-\lambda \tau }d\tau \times {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B)e^{-\lambda \tau }d\tau$(31)

$+\alpha \beta \mu {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}(A)e^{-\lambda \tau }d\tau \times {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B)e^{-\lambda \tau }d\tau =0.$

Consider the quantities ${\mathcal{Q}}_1$ and ${\mathcal{Q}}_2$ defined by

(32) ${\mathcal{Q}}_1=\alpha \beta (\lambda +\delta +\gamma ){\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A)e^{-\lambda \tau }d\tau \times {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B)e^{-\lambda \tau }d\tau$(32)

(33) ${\mathcal{Q}}_2=\alpha \beta \mu {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}(A)e^{-\lambda \tau }d\tau \times {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B)e^{-\lambda \tau }d\tau .$(33)

Then,

$\left|{\mathcal{Q}}_1+{\mathcal{Q}}_2\right|\le \left|{\mathcal{Q}}_1\right|+\left|{\mathcal{Q}}_2\right|$

$\le \alpha \beta \left|\lambda +\delta +\gamma \right|\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A)\times \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B)+\alpha \beta \mu \frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}(A)\times \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B)$

$={\mathcal{R}}_0^2\left|\lambda +d+b\right|\left|\lambda +\delta +\mu +\sigma \right|\left|\lambda +\delta +\gamma \right|$

$<\left|\lambda +d+b\right||\lambda +\delta +\mu +\sigma |\left|\lambda +\delta +\gamma \right|,$

and so, $\lambda$ cannot be the solution of (31). Hence, all characteristic roots have a negative real part, and therefore $E_0$ is locally asymptotically stable. □

Suppose that the incidence functions $f$ and $g$ satisfy the following conditions at endemic equilibrium

• (H4) $\alpha \beta \frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A^{\ast })\times \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B^{\ast })<\frac{1}{2}(d+b)(\delta +\mu +\sigma )$,

• (H5) $\alpha \beta \frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A^{\ast })\times \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B^{\ast })<\frac{1}{2}(d+b)(\delta +\mu +\sigma )$.

Theorem 4.2. If ${\mathcal{R}}_0>1$, (H4) and (H5) hold, then the endemic equilibrium is locally asymptotically stable.

Proof. The condition ${\mathcal{R}}_0>1$ guaranteeing the existence of an endemic equilibrium

$E^{\ast }=(A^{\ast },B^{\ast })$ with $A^{\ast }=(S_m^{\ast },I_H^{\ast },R_H^{\ast })$ and $B^{\ast }=(S_H^{\ast },I_m^{\ast })$.

The linearization of system (4) about the endemic equilibrium $E^{\ast }=(S_m^{\ast },I_H^{\ast },R_H^{\ast },S_H^{\ast },I_m^{\ast })$ is given by

(34) $\left\{\begin{array}{c}\frac{d{s}_m(t)}{dt}=-b{s}_m(t)-\alpha {\int }_0^h\eta (\tau )\left(\frac{\mathrm{\partial }f}{\mathrm{\partial }{S}_m}(A^{\ast })s_m(t)+\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A^{\ast })i_H(t-\tau )+\frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}(A^{\ast })r_H(t-\tau )\right)d\tau \\ \hspace{0.17em}\\ \frac{d{i}_m(t)}{dt}=\alpha {\int }_0^h\eta (\tau )\left(\frac{\mathrm{\partial }f}{\mathrm{\partial }{S}_m}(A^{\ast })s_m(t)+\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A^{\ast })i_H(t-\tau )+\frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}(A^{\ast })r_H(t-\tau )\right)d\tau -(d+b)i_m(t)\\ \hspace{0.17em}\\ \frac{d{s}_H(t)}{dt}=-\delta s_H(t)-\beta {\int }_0^h\eta (\tau )\left(\frac{\mathrm{\partial }g}{\mathrm{\partial }{S}_H}(B^{\ast })s_H(t)+\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B^{\ast })i_m(t-\tau )\right)d\tau +\gamma r_H(t)\\ \hspace{0.17em}\\ \frac{d{i}_H(t)}{dt}=\beta {\int }_0^h\eta (\tau )\left(\frac{\mathrm{\partial }g}{\mathrm{\partial }{S}_H}(B^{\ast })s_H(t)+\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B^{\ast })i_m(t-\tau )\right)d\tau -(\delta +\mu +\sigma )i_H(t)\\ \frac{d{r}_H(t)}{dt}=\mu i_H(t)-(\delta +\gamma )r_H(t).\end{array}\right.$(34)

We demonstrate that all the roots of the characteristic equation have a negative real part. So, as in the previous proof, using Equation (28) into (34) and canceling $e^{\lambda t}$ from each term and rearranging the following matrix equation yields

(35) $0=A_1\left[\begin{array}{c}{i}_m(0)\\ i_H(0)\\ r_H(0)\end{array}\right]$(35)

with

(36) $A_1=\left(\begin{array}{cccc}-(\lambda +d+b)& \alpha {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A^{\ast })e^{-\lambda \tau }d\tau & \alpha {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}(A^{\ast })e^{-\lambda \tau }d\tau & \hspace{0.17em}\\ \beta {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B^{\ast })e^{-\lambda \tau }d\tau & -(\lambda +\delta +\mu +\sigma )& 0& \hspace{0.17em}\\ 0& \mu & -(\lambda +\delta +\gamma )& \hspace{0.17em}\end{array}\right).$(36)

There exist non-zero solutions if and only if $det(A_1)=0$. Hence, we obtain the characteristic equation as follows:

$-(\lambda +d+b)(\lambda +\delta +\mu +\sigma )(\lambda +\delta +\gamma )$

(37) $+\alpha \beta (\lambda +\delta +\gamma ){\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A^{\ast })e^{-\lambda \tau }d\tau \times {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B^{\ast })e^{-\lambda \tau }d\tau$(37)

$+\alpha \beta \mu {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}(A^{\ast })e^{-\lambda \tau }d\tau \times {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B^{\ast })e^{-\lambda \tau }d\tau =0.$

Consider the following quantities

(38) ${\mathcal{Q}}_3=\alpha \beta (\lambda +\delta +\gamma ){\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A^{\ast })e^{-\lambda \tau }d\tau \times {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B^{\ast })e^{-\lambda \tau }d\tau$(38)

(39) ${\mathcal{Q}}_4=\alpha \beta \mu {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}(A^{\ast })e^{-\lambda \tau }d\tau \times {\int }_0^h\eta (\tau )\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B^{\ast })e^{-\lambda \tau }d\tau .$(39)

Hence, thanks to the hypothesis (H4) and (H5),

$\left|{\mathcal{Q}}_3+{\mathcal{Q}}_4\right|\le \left|{\mathcal{Q}}_3\right|+\left|{\mathcal{Q}}_4\right|$

$\le \alpha \beta \left|\lambda +\delta +\gamma \right|\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A^{\ast })\times \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B^{\ast })+\alpha \beta \mu \frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}(A^{\ast })\times \frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B^{\ast })$

$<\left|\lambda +d+b\right||\lambda +\delta +\mu +\sigma |\left|\lambda +\delta +\gamma \right|,$

and so the characteristic Equation (37) has solutions $\lambda$ with negative real part. In fact, the characteristic Equation (37) has solutions with non-negative real part if and only if all of the above inequalities are equations. Thus, the endemic equilibrium $E^{\ast }$ is locally asymptotically stable.

Now, let us analyze the global behavior of the equilibria.

5. The global stability of the equilibria

In this section, we study the global behavior of the endemic equilibrium. First, let us make the following useful assumptions:

(H6) For all $(S_m,I_H,R_H,S_H,I_m)\in {\mathbb{R}}_{+}^5$,

$f(S_m,I_H,R_H)\le I_H\frac{\mathrm{\partial }f}{\mathrm{\partial }{I}_H}(A)+R_H\frac{\mathrm{\partial }f}{\mathrm{\partial }{R}_H}(A),$

and

$g(S_H,I_m)\le I_m\frac{\mathrm{\partial }g}{\mathrm{\partial }{I}_m}(B).$

Theorem 5.1. The disease-free equilibrium $E_0=\left(\frac{{\mathrm{\Lambda }}_m}{b},0,0,\frac{{\mathrm{\Lambda }}_H}{\delta },0\right)$ is globally asymptotically stable in $\mathrm{\Delta }$ whenever ${\mathcal{R}}_0<1$.

Proof. The proof is based on using a comparison theorem (Lakshmikantham, Leela, & Martynyuk, 1989). Note that the equations of the infected components in system (4) can be expressed as in (Posny & Wang, 2014) like

(40) $\left[\begin{array}{c}{\dot{I}}_m\\ {\dot{I}}_H\\ {\dot{R}}_H\end{array}\right]\le \left(F-V\right)\left[\begin{array}{c}{I}_m\\ I_H\\ R_H\end{array}\right],$(40)

where the matrixes $F$ and $V$ are respectively given by (15) and (16).

Using the fact that all the eigenvalues of the matrix $F-V$ have negative real parts, it follows that the linearized differential inequality system (40) is stable whenever ${\mathcal{R}}_0<1$. Indeed, we can see that all eigenvalues of the matrix $F-V$ satisfy Equation (31). By the same reasoning as the proof of the Theorem 4.1, we deduce that the eigenvalues of the matrix $F-V$ have negative real parts since ${\mathcal{R}}_0<1$, (Van den Driessche & Watmough, 2002; Yang & Xiao, 2010).

Thus, $\left(I_m(t),I_H(t),R_H(t)\right)\to (0,0,0)$ as $t\to \mathrm{\infty }$ for system (40). Consequently, by the standard comparison theorem (Lakshmikantham et al., 1989), $\left(I_m(t),I_H(t),R_H(t)\right)\to (0,0,0)$ as $t\to \mathrm{\infty }$ and substituting $I_m(t)=I_H(t)=R_H(t)=0$ into system (4) we obtain $S_m(t)\to \frac{{\mathrm{\Lambda }}_m}{b}$ and $S_H(t)\to \frac{{\mathrm{\Lambda }}_H}{\delta }$ as $t\to \mathrm{\infty }$.

So, $\left(S_m(t),I_H(t),R_H(t),S_H(t),I_m(t)\right)\to \left(\frac{{\mathrm{\Lambda }}_m}{b},0,0,\frac{{\mathrm{\Lambda }}_H}{\delta },0\right)$ as $t\to \mathrm{\infty }$ for ${\mathcal{R}}_0<1$. Therefore, the disease-free equilibrium $E_0$ is globally asymptotically stable. ◴

Let us consider as in (Nakata et al., 2011), the positive function, $G(x)=x-1-lnx$ for all $x>0$ and the following region

${\mathrm{\Delta }}_1=\left\{(S_H,I_H,R_H,S_m,I_m)\in {\mathbb{R}}_{+}^5|\begin{array}{c}0<S_m(t)\le S_m^{\ast },0<I_m\le I_m^{\ast }\\ 0<S_H(t)\le S_H^{\ast },0<I_H(t)\le I_H^{\ast },0<R_H(t)\le R_H^{\ast }\end{array}\right\}.$

Theorem 5.2. If ${\mathcal{R}}_0>1$, then the endemic equilibrium $E^{\ast }$ is globally asymptotically stable in ${\mathrm{\Delta }}_1$.

Proof. Let us consider the following Lyapunov function

(41) $\mathcal{U}(t)={\mathcal{U}}_H(t)+{\mathcal{U}}_m(t)$(41)

with

(42) ${\mathcal{U}}_H(t)=\frac{S_H^{\ast }}{\beta g(S_H^{\ast },I_m^{\ast })}{U}_{S_H}(t)+\frac{I_H^{\ast }}{\beta g(S_H^{\ast },I_m^{\ast })}{U}_{I_H}(t)+U_{R_H}(t)$(42)

(43) ${\mathcal{U}}_m(t)=\frac{S_m^{\ast }}{\alpha f(S_m^{\ast },I_H^{\ast },R_H^{\ast })}{U}_{S_m}(t)+\frac{I_H^{\ast }}{\alpha f(S_m^{\ast },I_H^{\ast },R_H^{\ast })}{U}_{I_m}(t)$(43)

where $U_{S_H}(t),U_{I_H}(t)$, $U_{R_H}(t)$, $U_{S_m}(t)$, and $U_{I_m}(t)$ are respectively defined as follows

$U_{S_H}(t)=G\left(\frac{S_H(t)}{S_H^{\ast }}\right),U_{I_H}(t)=G\left(\frac{I_H(t)}{I_H^{\ast }}\right),U_{R_H}(t)=\frac{1}{2}{\left(R_H(t)-R_H^{\ast }\right)}^2,$

and

$U_{S_m}(t)=G\left(\frac{S_m(t)}{S_m^{\ast }}\right),U_{I_m}(t)=G\left(\frac{I_m(t)}{I_m^{\ast }}\right)$

Now, let us show that $\frac{d\mathcal{U}(t)}{dt}\le 0$.

First, we calculate $\frac{d{U}_{S_H}(t)}{dt}$ as:

$\frac{d{U}_{S_H}(t)}{dt}=\frac{S_{H^{\prime }}(t)}{S_H^{\ast }}-\frac{S_{H^{\prime }}(t)}{S_H(t)}$

Substituting $S_{H^{\prime }}(t)$ by its expression from the third equation of system (4) and also

${\mathrm{\Lambda }}_H=\delta S_H^{\ast }+\beta g(S_H^{\ast },I_m^{\ast })-\gamma R_H^{\ast },$

it then follows that

$\frac{d{U}_{S_H}(t)}{dt}=-\delta \left(1-\frac{S_H^{\ast }}{S_H(t)}\right)\left(\frac{S_H(t)}{S_H^{\ast }}-1\right)+\frac{\gamma R_H^{\ast }}{S_H^{\ast }}\left(1-\frac{S_H^{\ast }}{S_H(t)}\right)\left(\frac{R_H(t)}{R_H^{\ast }}-1\right)$

$+\beta \frac{g(S_H^{\ast },I_m^{\ast })}{S_H^{\ast }}{\int }_0^h\eta (\tau )\left(1-\frac{S_H^{\ast }}{S_H(t)}\right)\left(1-\frac{g(S_H(t),I_m(t-\tau )}{g(S_H^{\ast },I_m^{\ast })}\right)d\tau .$

We will use the following notations if need be

(44) $x_t=\frac{S_H(t)}{S_H^{\ast }},y_t=\frac{I_H(t)}{I_H^{\ast }},z_t=\frac{R_H(t)}{R_H^{\ast }},w_{t,\tau }=\frac{g(S_H(t),I_m(t-\tau ))}{g(S_H^{\ast },I_m^{\ast })}$(44)

and so $\frac{d{U}_{S_H}(t)}{dt}$ can be rewritten as

$\begin{array}{rl}& \frac{d{U}_{S_H}(t)}{dt}=-\delta \left(1-\frac{1}{x_t}\right)\left(x_t-1\right)+\frac{\gamma R_H^{\ast }}{S_H^{\ast }}\left(1-\frac{1}{x_t}\right)\left(z_t-1\right)\\ & +\beta \frac{g(S_H^{\ast },I_m^{\ast })}{S_H^{\ast }}{\int }_0^h\eta (\tau )\left(1-\frac{1}{x_t}\right)\left(1-w_{t,\tau }\right)d\tau \end{array}$

Secondly, we calculate $\frac{d{U}_{I_H}(t)}{dt}$. Indeed,

$\frac{d{U}_{I_H}(t)}{dt}=\frac{I_{H^{\prime }}(t)}{I_H^{\ast }}-\frac{I_{H^{\prime }}(t)}{I_H(t)}=\left(\frac{I_H(t)-I_H^{\ast }}{I_H^{\ast }{I}_H(t)}\right)I_{H^{\prime }}(t)$

At endemic equilibrium, thanks to the fourth equation of the system (4), we obtain that

$\delta +\mu +\sigma =\beta \frac{g(S_H^{\ast },I_m^{\ast })}{I_H^{\ast }}.$

Thus,

$\frac{d{U}_{I_H}(t)}{dt}=\frac{I_H(t)-I_H^{\ast }}{I_H^{\ast }{I}_H(t)}\left(\beta {\int }_0^h\eta (\tau )g(S_H(t),I_m(t-\tau ))d\tau -\beta \frac{g(S_H^{\ast },I_m^{\ast })}{I_H^{\ast }}{I}_H(t)\right)$

$=\beta \frac{g(S_H^{\ast },I_m^{\ast })}{I_H^{\ast }}{\int }_0^h\eta (\tau )\left(1-\frac{I_H^{\ast }}{I_H(t)}\right)\left(\frac{g(S_H(t),I_m(t-\tau ))}{g(S_H^{\ast },I_m^{\ast })}-\frac{I_H(t)}{I_H^{\ast }}\right)d\tau$

Using the above notations given by (44), it follows that

$\frac{d{U}_{I_H}(t)}{dt}=\beta \frac{g(S_H^{\ast },I_m^{\ast })}{I_H^{\ast }}{\int }_0^h\eta (\tau )\left(1-\frac{1}{y_t}\right)\left(w_{t,\tau }-y_t\right)d\tau .$

The above notations given by (44) lead to the following equations

$\begin{array}{rl}& \left(1-\frac{1}{x_t}\right)\left(1-w_{t,\tau }\right)+\left(1-\frac{1}{y_t}\right)\left(w_{t,\tau }-y_t\right)=2-\frac{1}{x_t}+\frac{w_{t,\tau }}{x_t}-y_t-\frac{w_{t,\tau }}{y_t}\\ & =-G\left(\frac{1}{x_t}\right)+\left(-ln\left(\frac{1}{x_t}\right)+1+\frac{w_{t,\tau }}{x_t}-y_t-\frac{w_{t,\tau }}{y_t}\right)\end{array}$

In addition,

$\begin{array}{rl}& -ln\left(\frac{1}{x_t}\right)+1+\frac{w_{t,\tau }}{x_t}-y_t-\frac{w_{t,\tau }}{y_t}+G\left(\frac{w_{t,\tau }}{x_t}\right)+G\left(\frac{w_{t,\tau }}{y_t}\right)=-ln\left(\frac{1}{x_t}\right)+1+2\frac{w_{t,\tau }}{x_t}-y_t-2\\ & -ln\frac{w_{t,\tau }}{x_t}-ln\frac{w_{t,\tau }}{y_t}\\ & =2\left(\frac{w_{t,\tau }}{x_t}-1-ln\frac{w_{t,\tau }}{x_t}\right)-\left(y_t-1-ln{y}_t\right)\\ & =2G\left(\frac{w_{t,\tau }}{x_t}\right)-G\left(y_t\right)\end{array}$

Consequently,

$-ln\left(\frac{1}{x_t}\right)+1+\frac{w_{t,\tau }}{x_t}-y_t-\frac{w_{t,\tau }}{y_t}=-G\left(\frac{w_{t,\tau }}{y_t}\right)-G\left(y_t\right)+G\left(\frac{w_{t,\tau }}{x_t}\right).$

And then,

$\left(1-\frac{1}{x_t}\right)\left(1-w_{t,\tau }\right)+\left(1-\frac{1}{y_t}\right)\left(w_{t,\tau }-y_t\right)=-G\left(\frac{1}{x_t}\right)-G\left(\frac{w_{t,\tau }}{y_t}\right)-G\left(y_t\right)+G\left(\frac{w_{t,\tau }}{x_t}\right)\le 0.$

Besides, we have

$-\delta \left(1-\frac{1}{x_t}\right)\left(x_t-1\right)+\frac{\gamma R_H^{\ast }}{S_H^{\ast }}\left(1-\frac{1}{x_t}\right)\left(z_t-1\right)\le 0$

Let us calculate now $\frac{d{U}_{R_H}(t)}{dt}$.

At endemic equilibrium point, we have

$\frac{{\mathrm{\Lambda }}_H}{\delta }=S_H^{\ast }-\frac{\gamma }{\delta }{R}_H^{\ast }+\frac{\beta }{\delta }g\left(S_H^{\ast },I_m^{\ast }\right)$

and

$g\left(S_H^{\ast },I_m^{\ast }\right)=\frac{(\delta +\mu +\sigma )(\delta +\gamma )}{\beta \mu }{R}_H^{\ast }.$

Hence,

$\frac{d{U}_{R_H}(t)}{dt}=(R_H(t)-R_H^{\ast })\frac{d{R}_H(t)}{dt}$

$=(R_H(t)-R_H^{\ast })\left(\mu \left(\frac{{\mathrm{\Lambda }}_H}{\delta }-S_H-R_H\right)-(\delta +\gamma )R_H\right)$

$=(R_H(t)-R_H^{\ast })\left(\mu \left(S_H^{\ast }-S_H(t)\right)+(\delta +\mu +\gamma )\left(R_H^{\ast }-R_H(t)\right)+\sigma \left(1+\frac{\gamma }{\delta }\right)R_H^{\ast }\right)$

$=-\mu \left(S_H^{\ast }{R}_H^{\ast }\right)\left(x_t-1\right)\left(z_t-1\right)-(\delta +\mu +\gamma ){\left(R_H^{\ast }\right)}^2{\left(z_t-1\right)}^2+\sigma \left(1+\frac{\gamma }{\delta }\right){\left(R_H^{\ast }\right)}^2\left(z_t-1\right)$

$\le 0.$

In summary, all these calculations lead to the following result,

(45) $\frac{d{\mathcal{U}}_H(t)}{dt}=\frac{d{U}_{S_H}(t)}{dt}+\frac{d{U}_{I_H}(t)}{dt}+\frac{d{U}_{R_H}(t)}{dt}\le 0$(45)

Secondly, we calculate $\frac{d{\mathcal{U}}_m(t)}{dt}$. As previously, let us begin by calculating $\frac{d{U}_{S_m}(t)}{dt}$.

$\frac{d{U}_{S_m}(t)}{dt}=\frac{S_{m^{\prime }}(t)}{S_m^{\ast }}-\frac{S_{m^{\prime }}(t)}{S_m(t)}$

$=\frac{S_m(t)-S_m^{\ast }}{S_m^{\ast }{S}_m(t)}\left({\mathrm{\Lambda }}_m-\alpha {\int }_0^h\eta (\tau )f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))d\tau -b{S}_m\right)$

For the sake of simplicity, we adopt the following notations

(46) $u_t=\frac{S_m(t)}{S_m^{\ast }};v_t=\frac{I_m(t)}{I_m^{\ast }},s_{t,\tau }=\frac{f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))}{f(S_m^{\ast },I_H^{\ast },R_H^{\ast })}$(46)

At the endemic state, from the first equation of system (4), we have

${\mathrm{\Lambda }}_m=\alpha f(S_m^{\ast },I_H^{\ast },R_H^{\ast })+b{S}_m^{\ast }.$

Hence

$\frac{d{U}_{S_m}(t)}{dt}=\frac{1}{S_m^{\ast }}\left(1-\frac{1}{u_t}\right)\left(\alpha f(S_m^{\ast },I_H^{\ast },R_H^{\ast }){\int }_0^h\eta (\tau )\left(1-s_{t,\tau }\right)d\tau +\frac{b}{S_m^{\ast }}\left(1-u_t\right)\right)$

$=\frac{\alpha f(S_m^{\ast },I_H^{\ast },R_H^{\ast })}{S_m^{\ast }}{\int }_0^h\eta (\tau )\left(1-\frac{1}{u_t}\right)\left(1-s_{t,\tau }\right)d\tau +\frac{b}{{(S_m^{\ast })}^2}\left(1-\frac{1}{u_t}\right)\left(1-u_t\right).$

Now, let us calculate $\frac{d{U}_{I_m}(t)}{dt}$

$\frac{d{U}_{I_m(t)}}{dt}=\frac{I_m^{\prime }(t)}{I_m^{\ast }}-\frac{I_m^{\prime }(t)}{I_m(t)}$

$=\frac{I_m(t)-I_m^{\ast }}{I_m^{\ast }{I}_m(t)}\left(\alpha {\int }_0^h\eta (\tau )f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))d\tau -(d+b)I_m(t)\right)$

At endemic equilibrium, thanks to the second equation of system (4), the following equality yields

$d+b=\frac{\alpha f(S_m^{\ast },I_H^{\ast },R_H^{\ast })}{I_m^{\ast }}.$

Consequently,

$\frac{d{U}_{I_m}(t)}{dt}=\frac{1}{I_m^{\ast }}\left(1-\frac{I_m^{\ast }}{I_m(t)}\right)\left(\alpha {\int }_0^h\eta (\tau )f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))d\tau -\frac{\alpha f(S_m^{\ast },I_H^{\ast },R_H^{\ast })}{I_m^{\ast }}{I}_m(t)\right)$

$=\frac{\alpha f(S_m^{\ast },I_H^{\ast },R_H^{\ast })}{I_m^{\ast }}{\int }_0^h\eta (\tau )\left(1-\frac{I_m^{\ast }}{I_m(t)}\right)\left(\frac{f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))}{f(S_m^{\ast },I_H^{\ast },R_H^{\ast })}-\frac{I_m(t)}{I_m^{\ast }}\right)d\tau$

$=\frac{\alpha f(S_m^{\ast },I_H^{\ast },R_H^{\ast })}{I_m^{\ast }}{\int }_0^h\eta (\tau )\left(1-\frac{1}{v_t}\right)\left(s_{t,\tau }-v_t\right)d\tau .$

Moreover, as from above, we show that

$\left(1-\frac{1}{u_t}\right)\left(1-s_{t,\tau }\right)+\left(1-\frac{1}{v_t}\right)\left(s_{t,\tau }-v_t\right)=-G\left(\frac{1}{u_t}\right)-G\left(\frac{s_{t,\tau }}{v_t}\right)-G\left(v_t\right)+G\left(\frac{s_{t,\tau }}{u_t}\right)\le 0$

And so,

(47) $\frac{d{\mathcal{U}}_m(t)}{dt}=\frac{d{U}_{S_m}(t)}{dt}+\frac{d{U}_{I_m}(t)}{dt}\le 0$(47)

In conclusion, thanks to the inequalities (45) and (47), we then have $\frac{d\mathcal{U}(t)}{dt}\le 0$. Thus, thanks to the LaSalle’s invariance principle (Koutou et al., 2018; Moulay et al., 2011), the endemic equilibrium is globally asymptotically stable. □

6. Numerical simulations

In this section, we use numerical tools to illustrate the theoretical results obtained in earlier sections for our model. For that, we consider as in (Nakata et al., 2011), the following corresponding discrete delay model of (4)

(48) $\left\{\begin{array}{c}\frac{d{S}_m(t)}{dt}={\mathrm{\Lambda }}_m-\alpha f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))-b{S}_m(t)\\ \hspace{0.17em}\\ \frac{d{I}_m(t)}{dt}=\alpha f(S_m(t),I_H(t-\tau ),R_H(t-\tau ))d\tau -(d+b)I_m(t)\\ \hspace{0.17em}\\ \frac{d{S}_H(t)}{dt}={\mathrm{\Lambda }}_H-\beta g(S_H(t),I_m(t-\tau ))-\delta S_H(t)+\gamma R_H(t)\\ \hspace{0.17em}\\ \frac{d{I}_H(t)}{dt}=\beta g(S_H(t),I_m(t-\tau ))-(\delta +\mu +\sigma )I_H(t)\\ \hspace{0.17em}\\ \frac{d{R}_H(t)}{dt}=\mu I_H(t)-(\delta +\gamma )R_H(t)\end{array}\right.$(48)

It is easy to see that (48) also has the same endemic equilibrium of (4) if and only if the condition ${\mathcal{R}}_0>1$ holds (Nakata et al., 2011). Moreover, in this computation, the functions $f$ and $g$ are chosen as follows $f(S_m,I_H,R_H)=S_m{I}_H+S_m{R}_H$ and $g(S_H,I_m)=S_H{I}_m$.

Figures 2–8 show the numerical simulations of model (4), with different values of the model parameters given in the corresponding captions, with the following initial conditions $S_m(0)=2$, $I_m(0)=0.12$, $S_H(0)=3$, $I_H(0)=0.15$, $R_H(0)=0.1$.

Figure 2. Density of infected mosquitoes and infected humans when, $\gamma =0.4$, $\delta =0.0187$, $\alpha =0.3$, $d=0.1$, $\mu =0.018$, $b=0.11$, ${\mathrm{\Lambda }}_H=0.29$, ${\mathrm{\Lambda }}_m=0.2$, $\beta =0.021$, $\sigma =0.011$ and $\tau =0.34$ and $\tau =1.34$ respectively for the sub-figures (a) and (b).

Figure 3. Periodic solutions with delay $\tau =0.1$ and $\gamma =0.01$, $\delta =0.011$, $\alpha =0.19$, $d=0.22$, $\mu =0.022$, $b=0.11$, ${\mathrm{\Lambda }}_H=0.4$, ${\mathrm{\Lambda }}_m=0.29$, $\beta =0.29$, $\sigma =0.011$.

Figure 4. Periodic solutions with a length incubation period $\tau =1.6$ and $\gamma =0.01$, $\delta =0.011$, $\alpha =0.19$, $d=0.22$, $\mu =0.022$, $b=0.11$, ${\mathrm{\Lambda }}_H=0.4$, ${\mathrm{\Lambda }}_m=0.29$, $\beta =0.29$, $\sigma =0.011$.

Figure 5. Density of infected mosquitoes versus susceptible mosquitoes; and susceptible mosquitoes versus infected mosquitoes when $\tau =0.1$ and $\gamma =0.01$, $\delta =0.011$, $\alpha =0.19$, $d=0.22$, $\mu =0.022$, $b=0.11$, ${\mathrm{\Lambda }}_H=0.4$, ${\mathrm{\Lambda }}_m=0.29$, $\beta =0.29$, $\sigma =0.011$.

Figure 6. Density of infected humans versus susceptible humans; and susceptible humans versus partially immune individuals when $\tau =0.1$ and $\gamma =0.01$, $\delta =0.011$, $\alpha =0.19$, $d=0.22$, $\mu =0.022$, $b=0.11$, ${\mathrm{\Lambda }}_H=0.4$, ${\mathrm{\Lambda }}_m=0.29$, $\beta =0.29$, $\sigma =0.011$.

Figure 7. Density of partially immune individuals versus susceptible; and infected humans versus partially immune individuals when $\tau =0.34$ and $\gamma =0.01$, $\delta =0.011$, $\alpha =0.19$, $d=0.22$, $\mu =0.022$, $b=0.11$, ${\mathrm{\Lambda }}_H=0.4$, ${\mathrm{\Lambda }}_m=0.29$, $\beta =0.29$, $\sigma =0.011$.

Figure 8. Limit cycles appearance when $\tau =1.34$ and $\gamma =0.01$, $\delta =0.011$, $\alpha =0.19$, $d=0.22$, $\mu =0.022$, $b=0.11$, ${\mathrm{\Lambda }}_H=0.4$, ${\mathrm{\Lambda }}_m=0.29$, $\beta =0.29$, $\sigma =0.011$.

Figure 2 shows that the disease-free equilibrium is global asymptotically stable whenever the threshold parameter ${\mathcal{R}}_0$ is less than unity and this, for any value of time-lag. For Figure 2, the parameter values which are chosen lead to ${\mathcal{R}}_0=0.6236$ which is less than unity. This result illustrates our Theorem 4.1. While we observe that the length of the incubation period decreases the number of infectious hosts. For instance, with $\tau =0.34$ there are about $3500$ infected mosquitoes and $1000$ infected humans (see sub-figure (a)). When we increase the incubation period to $\tau =1.34$, the number of infected hosts decreases and we obtain about $2500$ for mosquitoes and $900$ for humans (see sub-figure (b)).

Figures 3 – 8 illustrate the existence and the global stability of the unique endemic equilibrium of the model. The parameter values which are chosen lead to ${\mathcal{R}}_0=27.6087$ which is greater than unity. Theorem 4.2 then holds.

However, Figures 3 and 4 show that for some particular values of time-lag $\tau$, oscillations occur which can destabilize the system, and periodic solutions can arise due to Hopf bifurcation (see Figures 3 and 4). Moreover, one biological significance of oscillations is that each peak in the parasite marks a clinical episode of the disease for the host whatever it is (Ncube, 2013; Rihan & Anwar, 2012). This type of information would be essential for disease control and management strategies.

While Figures 5 – 8 show that a little incubation period such as $\tau =0.1$, the endemic steady of the model is exactly reached. Increasing this value can destabilize the system and give rise to a stable limit cycle. Indeed, the limit cycle is a trajectory for which the dynamics of the system would be constant over a cycle. This means that the epidemic gets a limit progression (Rasmussen, Wake, & Donaldson, 2003).

7. Discussions and conclusion

In this paper, we have studied in detail the persistence, the extinction, and the global asymptotic stability of an autonomous malaria transmission model with distributed time delay. The effect of the time-lag parameter $\tau$ on the behavior of the infection has been investigated. We determined the basic reproduction number ${\mathcal{R}}_0$ of the model and proved that it is not complete without considering the partially immune individual’s effects on the vectors and the malaria infection rate. And, we have shown that it is the threshold dynamics between the persistence and the extinction of the disease. More precisely, we have proven that if ${\mathcal{R}}_0<1$, the disease-free equilibrium is globally asymptotically stable. In this case, the epidemic will be extinct gradually. Similarly, if ${\mathcal{R}}_0>1$, the endemic equilibrium point is globally asymptotically stable. Hence, the disease will persist. Our theoretical studies deliver some new insights for understanding the propagation of the diseases and the information on biological networks, especially malaria.

Finally, a series of numerical simulations have been presented to show that the model agrees with the theoretical results. Furthermore, that numerical simulations show the impact of the delay on the model’s global behavior as follows: a great value of delay slows down the spreading of the disease, while a small value of time delay accelerates it. For particular values of time-lag $\tau$, oscillations occur which can destabilize the system, and periodic solutions can arise due to Hopf bifurcation. Although the underlying DDE model is simple, it displays very rich dynamics, such as quasi-periodic and chaotic patterns, and is suitable for large population densities (Sekiguchi & Ishiwata, 2010; Wang, Zhang, & Liu, 2018; Zhang & Teng, 2009).

Moreover, it emerges from our study that the incubation period is an important factor in the dynamics of malaria transmission. So, it can be used as a good parameter in fighting against malaria propagation.

Cover Image

Source: Author.

Additional information

Funding

The authors received no direct funding for this research.

Notes on contributors

Ousmane Koutou

Ousmane Koutou is a vigorous researcher PhD student at University Nazi BONI. His research interests lie in the field of Mathematical Biology and the present paper is a portion of his thesis under the supervision of Dr. Boureima SANGARÉ. He has published a number of research papers in reputed international journals, including, Advance in Difference Equation (Springer) and others journals.

Bakary Traoré

Bakary Traoré is a PhD student at University Nazi BONI under the supervision of Dr. Boureima SANGARÉ. A vigorous researcher in the area of Mathematical modelling, he has presented many papers in conferences and has published a number of research papers in reputed international journals, including Journal of Biological Dynamics (Taylor’s and Francis), Advance in Difference Equation (Springer) and others journals.

Boureima Sangaré

Dr Boureima Sangaré research interests lie in the field of Mathematical Modelling and four PhD students are currently working with him. He has published his research contributions in some internationally renowned journals whose publishers are Springer, Taylor & Francis and other journals.