# Analysis, state estimation and control of a malaria transmission model with semi-immune compartment for humans

Abderrahman Iggidr, Ali Moussaoui, Mohamed Oumoun, Souad Yacheur

► To cite this version:

Abderrahman Iggidr, Ali Moussaoui, Mohamed Oumoun, Souad Yacheur. Analysis, state estimation and control of a malaria transmission model with semi-immune compartment for humans. Journal of Mathematical Analysis and Applications, 2024, 535 (2), pp.128162. 10.1016/j.jmaa.2024.128162. hal-04038287

HAL Id: hal-04038287

https://inria.hal.science/hal-04038287

Submitted on 20 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analysis, state estimation and control of a malaria transmission model with semi-immune compartment for humans

Abderrahman Iggidr¹*, Ali Moussaoui², Mohamed Oumoun³, Souad Yacheur¹,²

¹ Université de Lorraine, CNRS, Inria, IECL, F-57000 Metz, France

² Department of Mathematics, University of Tlemcen, Algeria

³ Université Cadi Ayyad, École Nationale des Sciences Appliquées, Marrakech, Maroc

## Abstract

The aim of this paper is to study a mathematical model describing the dynamics of malaria transmission. The model tracks that humans gain a temporary immunity but they can still transmit the disease to susceptible mosquitoes with a probability less than that of infected humans. Furthermore, we take in consideration both return from the semi-immune class as a probability of loss of immunity and the return from the infected directly to the susceptible class. We investigate the existence of both trivial and non trivial equilibria and their stability. We show how to estimate the size of different human populations using the host incidence which is, in general, the only available information. To this end, we construct some observers or "software sensors" that allow to dynamically estimate the non-measured state variables. We also provide some strategies to control the disease using Lyapunov control functions.

**Mathematics Subject Classification:** 92D25, 92D30, 93B52, 93B53, 34H15.

**Keywords:** Epidemic models; Transmission dynamics; Stability analysis; Observer; State estimation; Control strategies.

## 1 Introduction:

Malaria is a vector-borne disease transmitted to humans primarily through the bite of infected female Anopheles. It is considered one of the most dangerous diseases in the world. The World Health Organization (WHO) estimates that there were about 219 million cases and 435000 related deaths in 2017 [27]. Malaria is one of the biggest risks that travelers can face while abroad, they could be at risk of this infection in 87 countries around the world, mainly in Africa, Asia and the

<footnote>*Corresponding author (email: abderrahman.iggidr@inria.fr).</footnote>
Americas [28]. Unfortunately, people at the highest risk of malaria infection include much more of infants, children under 5 years of age, pregnant women and people living with HIV/AIDS, as we know, this categories of people has low immunity. Many researches were done in the case of vector borne disease and specially Malaria models, let us recall some previously known results related to the present study. Since the first model proposed by Ross in 1911 [24], several models were elaborated by including different factors. MacDonald in 1957 [20] was the first one that introduced the latent compartment for the mosquitoes population, where Anderson and May in 1991 [1] had introduced this phase for the human population.

Environment, migration, immunity function, socio-economic factors and many other factors where introduced to better represent the reality of the malaria transmission process [24, 4, 29, 20].

Experimental evidence show that (60–90%) of humans in endemic area are asymptomatic carriers of the parasites (this is well explained in Ducrot et al. 2009 [10]; see also Chitnis et al. 2006 [7]; Chiyaka et al. 2007 [8]).

The model we consider in this paper is originated from the work [3], where the authors, consider the class of semi-immune as asymptomatic humans carriers who are less infectious to mosquitoes than symptomatic carriers. In the considered model, the human population is divided into three classes: susceptible individuals $S_h$ that represent people free of malaria capable of being infected, infectious $I_h$ bringing the parasite in the form *gametocyte* and semi-immune individuals $R_h$, which also carries the parasite in the form *gametocyte* but with one difference they have no symptoms and are less infectious to mosquitoes than symptomatic carriers $I_h$. The mosquito population is divided into two classes: susceptible $S_v$ and infected $I_v$ that carry the parasite in its form *sporozoite*. We denote by $H(t) = S_h(t) + I_h(t) + R_h(t)$ the total size of the human population. The total mosquito population will be denoted $V(t) = S_v(t) + I_v(t)$.

<table><caption>Table 1: The model parameters.</caption><thead><tr><th>Parameters</th><th>Description</th></tr></thead><tbody><tr><td>Λ<sub>h</sub></td><td>Recruitment into the human susceptible class</td></tr><tr><td>β<sub>vh</sub></td><td>Transmission probability from infectious mosquitoes to susceptible human.</td></tr><tr><td>β<sub>hv</sub></td><td>Transmission probability from infectious human to susceptible mosquitoes.</td></tr><tr><td>β̂<sub>hv</sub></td><td>Transmission probability from semi-immune human to susceptible mosquitoes.</td></tr><tr><td>γ<sub>h</sub></td><td>Rate of progression from the infectious to the semi-immune class.</td></tr><tr><td>α<sub>h</sub></td><td>Rate of recovery from being infectious.</td></tr><tr><td>δ<sub>h</sub></td><td>Disease induced death rate.</td></tr><tr><td>μ<sub>h</sub></td><td>Natural human death rate.</td></tr><tr><td>Λ<sub>v</sub></td><td>Recruitment into the vector susceptible class.</td></tr><tr><td>μ<sub>v</sub></td><td>Natural mosquitoes death rate.</td></tr><tr><td>V</td><td>The total mosquitoes population.</td></tr><tr><td>H</td><td>The total human population.</td></tr></tbody></table>

We take into account the return to the susceptible class $S_h$ from $I_h$ class as well as from the $R_h$ class. The model is then given as follows (using standard notations):
<figure><img src="image_2.png" /><figcaption>Figure 1: Compartmental diagram of the malaria disease model.</figcaption></figure>

$$
\left\{
\begin{aligned}
\frac{dS_h}{dt} &= \Lambda_h + \rho_h R_h + \alpha_h I_h - \beta_{vh} \frac{S_h}{H} I_v - \mu_h S_h, \\
\frac{dI_h}{dt} &= \beta_{vh} \frac{S_h}{H} I_v - \alpha_h I_h - \gamma_h I_h - \mu_h I_h - \delta_h I_h = \beta_{vh} \frac{S_h}{H} I_v - \epsilon_h I_h, \\
\frac{dR_h}{dt} &= \gamma_h I_h - \rho_h R_h - \mu_h R_h = \gamma_h I_h - \theta_h R_h, \\
\frac{dS_v}{dt} &= \Lambda_v - \left( \beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H} \right) S_v - \mu_v S_v, \\
\frac{dI_v}{dt} &= \left( \beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H} \right) S_v - \mu_v I_v.
\end{aligned}
\right.
\qquad (I)
$$

With initial conditions $S_h(0), S_v(0) > 0; I_h(0), I_v(0), R_h(0) \ge 0$.

To simplify the writing we denote by: $\epsilon_h = \delta_h + \alpha_h + \mu_h + \gamma_h$ and $\theta_h = \rho_h + \mu_h$. All the parameters of the model are positive except $\delta_h \ge 0$.

The mathematical analysis of the model proposed in [3] has addressed the local stability of the disease free equilibrium but not the global behavior of the system because of the potential occurrence of a backward bifurcation. In their analysis, they did not consider the behavior of the system when the basic reproduction number is greater than one.

Recently many results concerning the behavior of models similar to (I) were published, see for instance [22, 5, 16, 23, 2, 11]. In most of these references, it has been pointed that a backward bifurcation may occur when there is a disease-induced death rate (i.e., $\delta_h > 0$). In [2], a sufficient condition for the global asymptotic stability of the disease free equilibrium (DFE) has been given. The authors of [11] gave sufficient conditions for the global stability of the DFE as well as for the endemic equilibrium for a similar model but with $\hat{\beta}_{hv} = 0$. In [22] the authors proposed a new malaria transmission model with vector-bias effect is developed, they divided the human population into four classes SEIR and the mosquitoes population into three classes SEI. The crucial difference with Model (I) is that they didn't consider the possibility of passing from the infections class to the susceptible class nor the possibility of transmission of the disease from the "semi-recovered" $R_h$ to the susceptible mosquitoes $S_v$.

In this paper, besides the behavior of the model, we are interested in how to estimate the state variables and how to control the evolution of the disease in order to achieve some goals. A study of the estimation of the size of the populations $S_h$, $I_h$ and $R_h$ is made in order to be able to measure the complete state at any given time, using the only measure available in reality which is the number of newly infected humans per day, our goal is to take advantage of this information to give dynamic estimates of $S_h(t)$, $I_h(t)$ and $R_h(t)$ by employing an observer or a state estimator developed in the theory of automatic control. We also propose two control strategies, the first one consists in the use of a rate of over-mortality due to pesticides, which leads to a disappearance of the disease. The second one by applying a treatment to the infected individuals, knowing that
the treatment rate will be considered as a control u, we will show how to compute the treatment rate u as a function of the state so as to make the malaria-free equilibrium globally asymptotically stable, we provide a formula for the stabilizing feedback using a control Lyapunov function that we compute explicitly.

The rest of this article is organized as follows. Section 2 is dedicated to the study of the dynamical properties of Model (I): an analytic formula is derived for the basic reproduction number $R_0$, the disease free equilibrium (DFE) is shown to be globally asymptotically stable when $R_0$ is smaller than some constant which is smaller than one when there is a disease induced mortality and which is equal to one when there is no disease induced mortality, it is also shown that the system is uniformly persistent when $R_0 > 1$ and that the model has a unique endemic equilibrium whose local and global stability are studied. Section 3 is devoted to the state estimation problem. Some control strategies to fight against the disease are given in Section 4. A brief discussion is given in Section 5.

## 2 The dynamical model behavior

The evolution of the total human population is governed by $\frac{dH}{dt} = \Lambda_h - \mu_h H - \delta_h I_h$. When $H = \frac{\Lambda_h}{\mu_h + \delta_h}$, we have $\dot{H} = \Lambda_h - \mu_h \frac{\Lambda_h}{\mu_h + \delta_h} - \delta_h I_h \ge \Lambda_h - \mu_h \frac{\Lambda_h}{\mu_h + \delta_h} - \delta_h \frac{\Lambda_h}{\mu_h + \delta_h} \ge 0$, and if $H = \frac{\Lambda_h}{\mu_h}$ then $\dot{H} \le 0$. On the other hand the total vector population V satisfies $\dot{V} = \Lambda_v - \mu_v V$. Thus we have the following.

**Lemma 1.** The compact set $\Omega = \{ \frac{\Lambda_h}{\mu_h + \delta_h} \le H \le \frac{\Lambda_h}{\mu_h}, 0 \le S_v + I_v \le \frac{\Lambda_v}{\mu_v} = V^* \}$ is a positively invariant and attractive compact set for system (I).

### 2.1 The Disease-free equilibrium (DFE) and the basic reproduction number ($R_0$)

System (I) has a unique-disease free equilibrium (DFE) that will be denoted $E_0 = (\frac{\Lambda_h}{\mu_h}, 0, 0, \frac{\Lambda_v}{\mu_v}, 0)$.

Using the method described in [9] and denoting $x = (I_h, R_h, I_v)$ and $y = (S_h, S_v)$, we can rewrite our main system in the following form

$$ \begin{cases} \dot{x}_i = \mathcal{F}_i(x, y) - \mathcal{V}_i(x, y) & 1 \le i \le 3, \\ \dot{y}_j = g_j(x, y) & 1 \le j \le 2. \end{cases} \quad (1) $$

$\mathcal{F}(x, y)$ contains all the new infection terms:

$$ \mathcal{F} = \left( \beta_{vh} \frac{S_h}{H} I_v, 0, (\beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H}) S_v \right)^T, $$

and $\mathcal{V}$ contains all other transition terms:

$$ \mathcal{V} = \begin{pmatrix} -\epsilon_h I_h \\ \gamma_h I_h - \theta_h R_h \\ -\mu_v \end{pmatrix}. $$
Let $F = DF|_{E_0}$ and $V = DV|_{E_0}$ the Jacobian matrices of the maps $F$ and $V$, evaluated at the DFE,

$$F = \begin{pmatrix} 0 & 0 & \beta_{vh} \\ 0 & 0 & 0 \\ \beta_{hv} \frac{\Lambda_v}{\mu_v} \frac{\mu_h}{\Lambda_h} & \hat{\beta}_{hv} \frac{\Lambda_v}{\mu_v} \frac{\mu_h}{\Lambda_h} & 0 \end{pmatrix}, \quad V = \begin{pmatrix} -\epsilon_h & 0 & 0 \\ \gamma_h & -\theta_h & 0 \\ 0 & 0 & -\mu_v \end{pmatrix}. \qquad (2)$$

The basic reproduction number $R_0$ is the spectral radius of the next generation matrix

$$-FV^{-1} = \begin{pmatrix} 0 & 0 & \beta_{vh} \frac{1}{\mu_v} \\ 0 & 0 & 0 \\ \frac{\Lambda_v}{\mu_v} \frac{\mu_h}{\Lambda_h} \frac{1}{\epsilon_h} \left(\beta_{hv} + \frac{\gamma_h}{\theta_h} \hat{\beta}_{hv}\right) & \frac{\hat{\beta}_{hv}}{\theta_h} \frac{\Lambda_v}{\mu_v} \frac{\mu_h}{\Lambda_h} & 0 \end{pmatrix}.$$

Hence, we obtain

$$R_0 = \sqrt{\frac{\beta_{vh}}{\mu_v \epsilon_h} \frac{\Lambda_v}{\mu_v} \frac{\mu_h}{\Lambda_h} \left( \beta_{hv} + \frac{\gamma_h}{\theta_h} \hat{\beta}_{hv} \right)}, = \sqrt{\frac{\beta_{vh}}{\mu_v(\rho_h + \alpha_h + \mu_h + \delta_h)} \frac{\Lambda_v}{\mu_v} \frac{\mu_h}{\Lambda_h} \left( \beta_{hv} + \frac{\gamma_h}{\rho_h + \mu_h} \hat{\beta}_{hv} \right)}. \quad (3)$$

### 2.1.1 Global stability of the DFE

The total vector population $V$ satisfies $\dot{V} = \Lambda_v - \mu_v V$ that admits a unique equilibrium point $V^* = \frac{\Lambda_v}{\mu_v}$ which is globally asymptotically stable (GAS). Therefore it is possible to use [26, Theorem 3.1] and so the stability properties of System (I) are the same as those of the following system:

$$\left\{
\begin{aligned}
\frac{dH}{dt} &= \Lambda_h - \mu_h H - \delta_h I_h, \\
\frac{dI_h}{dt} &= \beta_{vh} \frac{I_v}{H} (H - I_h - R_h) - \epsilon_h I_h, \\
\frac{dR_h}{dt} &= \gamma_h I_h - \theta_h R_h, \\
\frac{dI_v}{dt} &= \left( \beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H} \right) (V^* - I_v) - \mu_v I_v, \text{ with } V^* = \frac{\Lambda_v}{\mu_v}.
\end{aligned}
\right.
\qquad (4)$$

Using the coordinates $(H, I_h, R_h, I_v)$, the DFE is given by $E_0 = (\frac{\Lambda_h}{\mu_h}, 0, 0, 0)$. The Jacobian at the DFE is given by:

$$J = \begin{pmatrix} -\mu_h & -\delta_h & 0 & 0 \\ 0 & -\epsilon_h & 0 & \beta_{vh} \\ 0 & \gamma_h & -\theta_h & 0 \\ 0 & \frac{\beta_{hv}V^*}{H^*} & \frac{\hat{\beta}_{hv}V^*}{H^*} & -\mu_v \end{pmatrix}, \text{ with } H^* = \frac{\Lambda_h}{\mu_h}, V^* = \frac{\Lambda_v}{\mu_v}.$$
The eigenvalues of $J$ are $-\mu_h$ and the eigenvalues of $J_1$ which is a Metzler matrix so it is stable if and only if it is invertible and $J_1^{-1} < 0$. We have

$$ \det(J_1) = \frac{-\epsilon_h \theta_h \mu_v \beta_{hv} H^* + V^* \beta_{vh} (\hat{\beta}_{hv} \gamma_h + \beta_{hv} \theta_h)}{H^*} \\ = (\mu_h + \alpha_h + \gamma_h + \delta_h) (\mu_h + \rho_h) \mu_v \beta_{hv} (\mathcal{R}_0^2 - 1) $$

$$ J_1^{-1} = \frac{1}{\det(J_1)} \begin{pmatrix} \mu_v \beta_{hv} (\mu_h + \rho_h) & \frac{\beta_{vh} \hat{\beta}_{hv} V^*}{H^*} & (\mu_h + \rho_h) \beta_{vh} \\ \gamma_h \mu_v \beta_{hv} & \frac{(\mu_h + \alpha_h + \gamma_h + \delta_h) \mu_v H^* - \beta_{hv} \beta_{vh} V^*}{H^*} & \beta_{vh} \gamma_h \\ \frac{V^* (\hat{\beta}_{hv} \gamma_h + (\mu_h + \rho_h) \beta_{hv})}{H^*} & \frac{(\mu_h + \alpha_h + \gamma_h + \delta_h) \hat{\beta}_{hv} V^*}{H^*} & (\mu_h + \alpha_h + \gamma_h + \delta_h) (\mu_h + \rho_h) \end{pmatrix} $$

If $\mathcal{R}_0^2 < 1$ then $J_1$ is invertible and $J_1^{-1} < 0$. If $\mathcal{R}_0^2 > 1$ then $J_1^{-1} > 0$ and so $J_1$ is not Hurwitz. Hence we have proved the following:

**Proposition 1.** The disease free equilibrium DFE $E_0$ is locally asymptotically stable if $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$.

To explore the global stability properties of the DFE, thanks to Proposition 1 it is sufficient to study its global attraction. Considering the infected individuals $I(t) = (I_h(t), R_h(t), I_v(t))^T$, we look for a condition under which

$$ I(t) \underset{t \to \infty}{\longrightarrow} 0 \quad \forall (S_h(0), I_h(0), R_h(0), S_v(0), I_v(0)) \in \Omega. $$

**Proposition 2.** If $\mathcal{R}_0 \le \sqrt{\frac{\mu_h}{\mu_h + \delta_h}}$ then the DFE is globally asymptotically stable.

*Proof.* Rewriting the equations of $I_H$, $R_H$, $I_V$ (with $\epsilon_h = \mu_h + \alpha_h + \gamma_h + \delta_h$ and $\theta_h = \mu_h + \rho_h$), and using the fact that $S_h \le H$, $S_v \le V^*$ and $\frac{\Lambda_h}{\mu_h + \delta_h} \le H$, we have

$$ \left\{ \begin{aligned} \frac{dI_h}{dt} &= \beta_{vh} \frac{I_v}{H} S_h - \epsilon_h I_h & \le \beta_{vh} I_v - \epsilon_h I_h \\ \frac{dR_h}{dt} &= \gamma_h I_h - \theta_h R_h \\ \frac{dI_v}{dt} &= \left( \beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H} \right) S_v - \mu_v I_v & \le \left( \beta_{hv} \frac{I_h}{m} + \hat{\beta}_{hv} \frac{R_h}{m} \right) V^* - \mu_v I_v \end{aligned} \right. $$

with $m = \frac{\Lambda_h}{\mu_h + \delta_h}$. Let $A$ be the following matrix.

$$ A = \begin{pmatrix} -\epsilon_h & 0 & \beta_{vh} \\ \gamma_h & -\theta_h & 0 \\ \beta_{hv} \frac{V^*}{m} & \hat{\beta}_{hv} \frac{V^*}{m} & -\mu_v \end{pmatrix}, $$

then we have the following relation: $\frac{dI(t)}{dt} \le A I(t)$. The matrix $A$ is a Metzler matrix. To show that $A$ is Hurwitz it is necessary and sufficient to show that $A$ is invertible and that $-A^{-1} > 0$.
We have $det(A) = \frac{V^*\gamma\hat{\beta}_{hv}\beta_{vh} + V^*\theta_h\beta_{hv}\beta_{vh} - \epsilon_h m\theta_h\mu_v}{m}$, and

$$A^{-1} = \frac{m}{\beta_{vh}V^*(\gamma\hat{\beta}_{hv} + \theta_h\beta_{hv}) - \epsilon_h m\theta_h\mu_v} \begin{pmatrix} \theta_h\mu_v m & \beta_{vh}\hat{\beta}_{hv}V^* & \theta_h\beta_{vh}m \\ \gamma\mu_v m & \epsilon_h m\mu_v - V^*\beta_{hv}\beta_{vh} & \gamma\beta_{vh}m \\ V^*(\gamma\hat{\beta}_{hv} + \theta_h\beta_{hv}) & \epsilon_h\hat{\beta}_{hv}V^* & \epsilon_h m\theta_h \end{pmatrix}$$

Therefore, $A^{-1}$ exists and $-A^{-1} > 0$ if and only if

$$\beta_{vh}V^*(\gamma\hat{\beta}_{hv} + \theta_h\beta_{hv}) - \epsilon_h m\theta_h\mu_v < 0 \text{ and } \epsilon_h m\mu_v - V^*\beta_{hv}\beta_{vh} \geq 0.$$

Using the expressions $\mathcal{R}_0^2 = \frac{\beta_{vh}}{\mu_v\epsilon_h}\frac{\Lambda_v}{\mu_v}\frac{\mu_h}{\Lambda_h}(\beta_{hv} + \frac{\gamma_h}{\theta_h}\hat{\beta}_{hv})$, $V^* = \frac{\Lambda_v}{\mu_v}$, and $m = \frac{\Lambda_h}{\mu_h + \delta_h}$, we obtain

$\beta_{vh}V^*(\gamma\hat{\beta}_{hv} + \theta_h\beta_{hv}) - \epsilon_h m\theta_h\mu_v = \Lambda_h \left(\frac{\mathcal{R}_0^2}{\mu_h} - \frac{1}{\mu_h + \delta_h}\right) < 0$ if $\mathcal{R}_0^2 < \frac{\mu_h}{\mu_h + \delta_h}$. On the other hand,

$$\epsilon_h m\mu_v - V^*\beta_{hv}\beta_{vh} = \epsilon_h \mu_v \frac{\Lambda_h}{\mu_h + \delta_h} - \beta_{hv}\beta_{vh}\frac{\Lambda_v}{\mu_v} = \frac{\epsilon_h \mu_v \Lambda_h}{\mu_h + \delta_h} \left(1 - \frac{\beta_{hv}\beta_{vh}\Lambda_v \mu_h}{\epsilon_h \mu_v^2 \Lambda_h} \frac{\mu_h + \delta_h}{\mu_h}\right).$$

We have $\frac{\beta_{hv}\beta_{vh}\Lambda_v\mu_h}{\epsilon_h \mu_v^2\Lambda_h} = \mathcal{R}_0^2 - \frac{\beta_{vh}}{\mu_v\epsilon_h}\frac{\Lambda_v}{\mu_v}\frac{\mu_h}{\Lambda_h}\frac{\gamma_h}{\theta_h}\hat{\beta}_{hv}$, so we can write

$$1 - \frac{\beta_{hv}\beta_{vh}\Lambda_v\mu_h}{\epsilon_h \mu_v^2 \Lambda_h} \frac{\mu_h + \delta_h}{\mu_h} = 1 - \mathcal{R}_0^2 \frac{\mu_h + \delta_h}{\mu_h} + \frac{\beta_{vh}}{\mu_v\epsilon_h}\frac{\Lambda_v}{\mu_v}\frac{\mu_h}{\Lambda_h}\frac{\gamma_h}{\theta_h}\hat{\beta}_{hv}\frac{\mu_h + \delta_h}{\mu_h} \geq 0 \text{ if } \mathcal{R}_0^2 \leq \frac{\mu_h}{\mu_h + \delta_h}.$$

Hence, $A^{-1}$ exists and $-A^{-1} > 0$ if and only if $\mathcal{R}_0^2 < \frac{\mu_h}{\mu_h + \delta_h}$. $\square$

When $\mathcal{R}_0^2 = \frac{\mu_h}{\mu_h + \delta_h}$, the above method does not allow to conclude since in this case the matrix A is not invertible. So we propose an alternative proof using a Lyapunov function.

**A Lyapunov proof.** We consider the candidate Lyapunov function for system (4) on the compact set $\Omega$:

$$W = \frac{(\rho + \mu_h) \mu_v^2 \Lambda_h}{\mu_h \Lambda_v \beta_{vh} \hat{\beta}_{hv}} \mathcal{R}_0^2 I_h + R_h + \frac{(\rho + \mu_h) \mu_v \Lambda_h}{\hat{\beta}_{hv} \Lambda_v (\mu_h + \delta)} I_v$$

The derivative of W along the solutions of system (4) is

$$\dot{W} = \underbrace{\left( - \frac{(\rho + \mu_h) \mu_v^2 \Lambda_h \mathcal{R}_0^2 (\alpha + \gamma + \mu_h + \delta)}{\mu_h \Lambda_v \beta_{vh} \hat{\beta}_{hv}} + \gamma + \frac{(\rho + \mu_h) \Lambda_h \beta_{hv}}{\hat{\beta}_{hv} (\mu_h + \delta) H} \right)}_{B} I_h - \underbrace{\left( \rho + \mu_h + \frac{(\rho + \mu_h) \Lambda_h}{(\mu_h + \delta) H} \right)}_{C} R_h + \underbrace{\left( \frac{(\rho + \mu_h) \mu_v^2 \Lambda_h \mathcal{R}_0^2}{\mu_h \Lambda_v \hat{\beta}_{hv}} - \frac{(\rho + \mu_h) \mu_v^2 \Lambda_h}{\hat{\beta}_{hv} \Lambda_v (\mu_h + \delta)} \right)}_{D} I_v - \underbrace{\left( \frac{(\rho + \mu_h) \mu_v^2 \Lambda_h \mathcal{R}_0^2}{\mu_h \Lambda_v \hat{\beta}_{hv} H} + \frac{(\rho + \mu_h) \mu_v \Lambda_h \beta_{hv}}{\hat{\beta}_{hv} \Lambda_v (\mu_h + \delta) H} \right)}_{E} I_h + \underbrace{\left( \frac{(\rho + \mu_h) \mu_v^2 \Lambda_h \mathcal{R}_0^2}{\mu_h \Lambda_v \hat{\beta}_{hv} H} + \frac{(\rho + \mu_h) \mu_v \Lambda_h}{\Lambda_v (\mu_h + \delta) H} \right)}_{F} R_h \quad I_v$$

We have

$$A = - \frac{(\rho + \mu_h) \mu_v^2 \Lambda_h \mathcal{R}_0^2 (\alpha + \gamma + \mu_h + \delta)}{\mu_h \Lambda_v \beta_{vh} \hat{\beta}_{hv}} + \gamma + \frac{(\rho + \mu_h) \Lambda_h \beta_{hv}}{\hat{\beta}_{hv} (\mu_h + \delta) H} \\ = - \frac{(\rho + \mu_h) \beta_{hv} ((\mu_h + \delta) H - \Lambda_h)}{\hat{\beta}_{hv} (\mu_h + \delta) H} \leq 0 \text{ since } \frac{\Lambda_h}{\mu_h + \delta_h} \leq H.$$
$$B = - \left( \rho + \mu_h + \frac{(\rho + \mu_h) \Lambda_h}{(\mu_h + \delta) H} \right) = -(\rho + \mu_h) \left( 1 - \frac{\Lambda_h}{(\mu_h + \delta)H} \right) \le 0 \text{ since } \frac{\Lambda_h}{\mu_h + \delta_h} \le H.$$

$$C = \frac{(\rho + \mu_h)\mu_v^2 \Lambda_h \mathcal{R}_0^2}{\mu_h \Lambda_v \hat{\beta}_{hv}} - \frac{(\rho + \mu_h)\mu_v^2 \Lambda_h}{\hat{\beta}_{hv} \Lambda_v (\mu_h + \delta)} = \frac{(\rho + \mu_h)\mu_v^2 \Lambda_h}{\mu_h \Lambda_v \hat{\beta}_{hv}} \left( \mathcal{R}_0^2 - \frac{\mu_h}{\mu_h + \delta} \right) \le 0 \text{ if } \mathcal{R}_0 \le \sqrt{\frac{\mu_h}{\mu_h + \delta}}.$$

$D > 0$ and $E > 0$.

Hence, $\dot{W} \le 0$ if $\mathcal{R}_0 \le \sqrt{\frac{\mu_h}{\mu_h + \delta}}$, and $\dot{W} = 0$ if and only if $(I_h = R = I_v = 0)$ or $(\mathcal{R}_0 = \sqrt{\frac{\mu_h}{\mu_h + \delta}}$ and $I_h = R = 0)$. In the latter case, it is easy to show that the largest invariant set contained in $\{x \in \Omega : \dot{W}(x) = 0\}$ is reduced to the DFE. Thus, by LaSalle Invariance Principle, we conclude that the DFE is GAS. $\square$

**Remark 1.** Proposition 2 could be proved using Theorem 4.3 in [15].

Indeed, with the notations of [15], we have

$$\dot{x}_2 = \begin{pmatrix} \dot{I}_h \\ \dot{R}_h \\ \dot{I}_v \end{pmatrix} = \underbrace{\begin{pmatrix} -\epsilon_h & 0 & \beta_{vh} \frac{S_h}{H} \\ \gamma_h & -\theta_h & 0 \\ \beta_{hv} \frac{S_v}{H} & \hat{\beta}_{hv} \frac{S_v}{H} & -\mu_v \end{pmatrix}}_{A_2(x)} \begin{pmatrix} I_h \\ R_h \\ I_v \end{pmatrix}, \quad (5)$$

and

$$\bar{A}_2 = \begin{pmatrix} -\epsilon_h & 0 & \beta_{vh} \\ \gamma_h & -\theta_h & 0 \\ \beta_{hv} \frac{V^*}{m} & \hat{\beta}_{hv} \frac{V^*}{m} & -\mu_v \end{pmatrix}.$$

**Remark 2.** The sufficient condition for global stability of the DFE given by Proposition 2 is exactly the same as the one given by Theorem 6 in [2].

When there is no disease induced mortality, Proposition 2 implies the following result.

**Corollary 1.** Suppose $\delta_h = 0$, then the DFE is globally asymptotically stable if and only if $\mathcal{R}_0 \le 1$.

## 2.2 The uniform persistence

**Theorem 1.** If $\mathcal{R}_0 > 1$, then System (I) is uniformly persistent.

*Proof.* Let $X$ be the $\omega$-limit set of $\Omega$. Since $\Omega$ is positively invariant, we have that $X \subset \Omega$. The largest invariant compact set $M$ in $\partial\Omega$ (the boundary of $\Omega$) is reduced to the DFE, that is $M = \{E_0\}$ and so it is isolated in $X$. On the other hand, when $\mathcal{R}_0 > 1$, $E_0$ is unstable and the stable manifold of $E_0$ is contained in $\partial\Omega$. Therefore, conditions (1) and (2) of Theorem 4.1 in [14] are satisfied and so uniform persistence holds. $\square$

It follows that, when $\mathcal{R}_0 > 1$, there exists $r > 0$ such that for all initial conditions in $\tilde{\Omega}$ the interior of $\Omega$, one has

$$\liminf_{t \to \infty} S_h(t) > r, \liminf_{t \to \infty} I_h(t) > r, \liminf_{t \to \infty} R_h(t) > r, \liminf_{t \to \infty} S_v(t) > r, \text{ and } \liminf_{t \to \infty} I_v(t) > r.$$
## 2.3 Existence and uniqueness of the endemic equilibrium when $R_0 > 1$

A state $(\bar{H}, \bar{I}_h, \bar{R}_h, \bar{I}_v)$ is an equilibrium for System (4) if and only if the following relations are satisfied:

$$ egin{cases} 
\Lambda_h - \mu_h \bar{H} - \delta_h \bar{I}_h = 0, \\
\beta_{vh} \frac{\bar{I}_v}{\bar{H}} (\bar{H} - \bar{I}_h - \bar{R}_h) - (\alpha_h + \gamma_h + \mu_h + \delta_h) \bar{I}_h = 0, \\
\gamma_h \bar{I}_h - (\rho_h + \mu_h) \bar{R}_h = 0, \\
\left( \beta_{hv} \frac{\bar{I}_h}{\bar{H}} + \hat{\beta}_{hv} \frac{\bar{R}_h}{\bar{H}} \right) \left( \frac{\Lambda_v}{\mu_v} - \bar{I}_v \right) - \mu_v \bar{I}_v = 0.
\end{cases} 
\quad (6) $$

Using the first equation of (6) and the fact that $\frac{\Lambda_h}{\mu_h + \delta_h} \le H \le \frac{\Lambda_h}{\mu_h}$, we have that $\bar{I}_h$ satisfies

$$ 0 < \bar{I}_h \le \frac{\Lambda_h}{\mu_h + \delta_h}. \quad (7) $$

Substituting $\bar{R}_h = \frac{\gamma_h}{\rho_h + \mu_h} \bar{I}_h$, and $\bar{H} = \frac{\Lambda_h - \delta_h \bar{I}_h}{\mu_h}$ in the second and in the fourth relations above, we obtain

$$ egin{cases} 
\bar{I}_h = \frac{\Lambda_h \mu_v \bar{I}_v}{\mu_h \left( \beta_{hv} + \frac{\gamma_h \hat{\beta}_{hv}}{\rho_h + \mu_h} \right) \left( \frac{\Lambda_v}{\mu_v} - \bar{I}_v \right) + \delta_h \mu_v \bar{I}_v}, \\
\bar{I}_v = \frac{(\Lambda_h - \delta_h \bar{I}_h)(\alpha_h + \gamma_h + \mu_h + \delta_h) \bar{I}_h}{\mu_h \beta_{vh} \left( \frac{\Lambda_h - \delta_h \bar{I}_h}{\mu_h} - \left(1 + \frac{\gamma_h}{\rho_h + \mu_h}\right) \bar{I}_h \right)}
\end{cases} $$

Therefore, $(\bar{I}_h, \bar{I}_v)$ is a fixed point of the function

$$ f(x, y) = \begin{pmatrix} \frac{\Lambda_h \mu_v y}{\mu_h \left( \beta_{hv} + \frac{\gamma_h \hat{\beta}_{hv}}{\rho_h + \mu_h} \right) \left( \frac{\Lambda_v}{\mu_v} - y \right) + \delta_h \mu_v y} \\ \frac{(\Lambda_h - \delta_h x)(\alpha_h + \gamma_h + \mu_h + \delta_h)x}{\mu_h \beta_{vh} \left( \frac{\Lambda_h}{\mu_h} - \left(1 + \frac{\delta_h}{\mu_h} + \frac{\gamma_h}{\rho_h + \mu_h}\right)x \right)} \end{pmatrix} = \begin{pmatrix} f_1(y) \\ f_2(x) \end{pmatrix}. $$

We can write $f_2(x) = \frac{(\Lambda_h - \delta_h x)(\alpha_h + \gamma_h + \mu_h + \delta_h)x}{\beta_{vh} \left( \Lambda_h - \left( \mu_h + \delta_h + \frac{\gamma_h \mu_h}{\theta_h} \right)x \right)} = \frac{\theta_h (\Lambda_h - \delta_h x)(\alpha_h + \gamma_h + \mu_h + \delta_h)x}{\beta_{vh} ((\Lambda_h - (\mu_h + \delta_h)x)\theta_h - \gamma_h \mu_h x)}$.

For $x \in [0, \frac{\Lambda_h}{\mu_h}]$, $f_2(x) \ge 0$ if and only if $x < \frac{\Lambda_h}{\mu_h + \delta_h + \frac{\gamma_h \mu_h}{\rho_h + \mu_h}}$. For $y \in [0, \frac{\Lambda_v}{\mu_v}]$, $f_1(y) \ge 0$. $f(0, 0) = (0, 0)^T$.

The Jacobian of $f$ is $Jf(x, y) = \begin{pmatrix} 0 & f'_1(y) \\ f'_2(x) & 0 \end{pmatrix}$ with

$$ f'_1(y) = \frac{\theta_h \Lambda_h \mu_v^2 \Lambda_v \mu_h (\hat{\beta}_{hv} \gamma_h + \theta_h \beta_{hv})}{\left( - (\hat{\beta}_{hv} \gamma_h + \theta_h \beta_{hv}) (\mu_v y - \Lambda_v) \mu_h + \delta_h \theta_h \mu_v^2 y \right)^2}, $$
$$f_2'(x) = \frac{\theta_h (\alpha_h + \gamma_h + \mu_h + \delta_h) \left( (\delta_h (\delta_h + \mu_h) x^2 - 2 \delta_h \Lambda_h x + \Lambda_h^2) \theta_h + \delta_h \gamma_h \mu_h x^2 \right)}{\beta_{vh} (((\delta_h + \mu_h) x - \Lambda_h) \theta_h + \gamma_h \mu_h x)^2}.$$

$f_1'$ and $f_2'$ are positive on the sets where they are defined. So $f_2$ is increasing and tends to $+\infty$ when $x \to \frac{\Lambda_h}{\mu_h + \delta_h + \frac{\gamma_h \mu_h}{\rho_h + \mu_h}}$, hence there exists $x_1 < \frac{\Lambda_h}{\mu_h + \delta_h + \frac{\gamma_h \mu_h}{\rho_h + \mu_h}}$ such that $f_2(x_1) = \frac{\Lambda_v}{\mu_v}$ and $f_2(x) \le \frac{\Lambda_v}{\mu_v}$ for all $x \in [0, x_1]$. Therefore, $f_1 \circ f_2$ is well defined on $[0, x_1]$ and is increasing. Moreover $f_1 \circ f_2(0) = 0$ and $f_1 \circ f_2(x_1) = f_1(\frac{\Lambda_v}{\mu_v}) = \frac{\Lambda_h}{\delta_h} > x_1$. We also have that $(f_1 \circ f_2)'(0) = \frac{1}{\mathbb{R}_0^2} < 1$. Hence $f_1 \circ f_2$ has at least a fixed point $\bar{x} \in (0, x_1)$.

We shall prove the uniqueness of the fixed point in $(0, x_1)$ by proving that $f_1 \circ f_2$ is convex on $[0, x_1]$.

$$(f_1 \circ f_2)'(x) = f_1'(f_2(x)). f_2'(x) > 0.$$

$$(f_1 \circ f_2)''(x) = f_1''(f_2(x)). f_2''(x) + f_1'(f_2(x)). f_2'(x)$$

$$f_2''(x) = \frac{2(\rho_h + \mu_h)(\alpha_h + \gamma_h + \mu_h + \delta_h)\Lambda_h^2 \mu_h(\gamma_h + \rho_h + \mu_h)}{\beta_{vh}(\mu_h + \delta_h + \frac{\gamma_h \mu_h}{\rho_h + \mu_h}) \left( \frac{\Lambda_h}{\mu_h + \delta + \frac{\gamma_h \mu_h}{\rho_h + \mu_h}} - x \right)^3} > 0 \text{ for } x < \frac{\Lambda_h}{\mu_h + \delta + \frac{\gamma_h \mu_h}{\rho_h + \mu_h}}.$$

$$f_1''(y) = \frac{2\theta_h \Lambda_h \Lambda_v \mu_h \mu_v^3 (\gamma_h \hat{\beta}_{hv} + \theta_h \beta_{hv}) (\delta_h \theta_h \mu_v - (\gamma_h \hat{\beta}_{hv} + \theta_h \beta_{hv}) \mu_h)}{((\hat{\beta}_{hv} \gamma_h + \theta_h \beta_{hv}) (\Lambda_v - \mu_v y) \mu_h + \delta_h \theta_h \mu_v^2 y)^3} \text{ of fixed sign for } y \le \frac{\Lambda_v}{\mu_v}.$$

$$f_1'''(y) = \frac{6\theta_h \Lambda_h \Lambda_v \mu_h \mu_v^4 (\gamma_h \hat{\beta}_{hv} + \theta_h \beta_{hv}) (\delta_h \theta_h \mu_v - (\gamma_h \hat{\beta}_{hv} + \theta_h \beta_{hv}) \mu_h)^2}{((\hat{\beta}_{hv} \gamma_h + \theta_h \beta_{hv}) (\Lambda_v - \mu_v y) \mu_h + \delta_h \theta_h \mu_v^2 y)^4} > 0 \text{ for } y \le \frac{\Lambda_v}{\mu_v}.$$

* If $\delta_h \theta_h \mu_v > (\gamma_h \hat{\beta}_{hv} + \theta_h \beta_{hv}) \mu_h$, then $f_1'' > 0$ and so $(f_1 \circ f_2)'' > 0$. Therefore $(f_1 \circ f_2)$ is convex which ensures the uniqueness of the fixed point in $(0, x_1)$.

* If $\delta_h \theta_h \mu_v < (\gamma_h \hat{\beta}_{hv} + \theta_h \beta_{hv}) \mu_h$ then $f_1'' < 0$ and so $f_1'$ is decreasing. Using the facts that $f_2''$ is positive and increasing, $f_2$ is positive and increasing, $f_1'$ is decreasing, $f_2'$ is increasing, and $f_1''$ is increasing, we have that for all $x \in [0, x_1]$:

$$(f_1 \circ f_2)''(x) \ge f_1''(f_2(0)) f_2'(0)^2 + f_1'(f_2(x_1)) f_2''(0) = f_1''(f_2(0)) f_2'(0)^2 + f_1' \left(\frac{\Lambda_v}{\mu_v}\right) f_2''(0) := A$$

Some computations allow to obtain for A, using the notations $\epsilon_h = \alpha_h + \gamma_h + \mu_h + \delta_h$ and $\theta_h = \rho_h + \mu_h$,

$$A = \frac{2 \mu_v^4 \theta^2 \Lambda_h \epsilon_h^2 \delta_h}{\Lambda_v^2 \mu_h^2 (\gamma_h \hat{\beta}_{hv} + \theta_h \beta_{hv})^2 \beta_{vh}^2} \left( \frac{\mu_h^4 (\gamma_h \hat{\beta}_{hv} + \theta_h \beta_{hv})^3 (\gamma_h + \theta_h) \beta_{vh} \Lambda_v}{\mu_v^4 \theta_h^4 \Lambda_h \epsilon_h \delta_h^3} + \frac{\mu_h (\gamma_h \hat{\beta}_{hv} + \beta_{hv} \theta)}{\delta_h \mu_v \theta} - 1 \right).$$

Since $\delta_h \theta_h \mu_v < (\gamma_h \hat{\beta}_{hv} + \theta_h \beta_{hv}) \mu_h$, we have $A > 0$ and so $(f_1 \circ f_2)''(x) > 0$ for all $x \in [0, x_1]$. In both cases, $(f_1 \circ f_2)$ is convex on $[0, x_1]$ which ensures the uniqueness of the fixed point in $(0, x_1)$.
For $x_1 \le x < \frac{\Lambda_h}{\mu_h + \delta + \frac{\gamma_h \mu_h}{\rho_h + \mu_h}}$, we have $f_2(x) \in [\frac{\Lambda_v}{\mu_v}, \infty)$, and $f_1 \circ f_2(x) \ge \frac{\Lambda_h}{\delta_h} > \frac{\Lambda_h}{\mu_h + \delta}$. Therefore, $f_1 \circ f_2(x)$ cannot have a fixed point in $[x_1, \frac{\Lambda_h}{\mu_h + \delta + \frac{\gamma_h \mu_h}{\rho_h + \mu_h}})$ since a fixed point must satisfy relation (7).

So the function $f$ has a unique fixed point $(\bar{x}, f_2(\bar{x}))$ that belongs to $(0, x_1) \times (0, \frac{\Lambda_v}{\mu_v})$. We have then proved the following:

**Theorem 2.** If $\mathcal{R}_0 > 1$ then System (4) has a unique endemic equilibrium EE.

## 2.4 Stability of the EE

**Proposition 3.** Suppose $\mathcal{R}_0 > 1$. The endemic equilibrium EE is asymptotically stable if

$$ \frac{\delta_h}{\mu_h} \le \frac{\beta_{hv} + \frac{\hat{\beta}_{hv} \gamma_h}{\theta_h}}{\mu_v}. \qquad (8) $$

*Proof.* The Jacobian at the EE (with $\epsilon_h = \alpha_h + \gamma_h + \mu_h + \delta_h$ and $\theta_h = \rho_h + \mu_h$):

$$ J_{EE} = \begin{bmatrix} -\mu_h & -\delta & 0 & 0 \\ -\frac{\beta_{vh}\bar{I}_v(\bar{H}-\bar{I}_h-\bar{R}_h)}{\bar{H}^2} + \frac{\beta_{vh}\bar{I}_v}{\bar{H}} & -\frac{\beta_{vh}\bar{I}_v}{\bar{H}} - \epsilon_h & -\frac{\beta_{vh}\bar{I}_v}{\bar{H}} & \frac{\beta_{vh}(\bar{H}-\bar{I}_h-\bar{R}_h)}{\bar{H}} \\ 0 & \gamma_h & -\theta_h & 0 \\ \left(-\frac{\beta_{hv}\bar{I}_h}{\bar{H}^2} - \frac{\hat{\beta}_{hv}\bar{R}_h}{\bar{H}^2}\right)\left(\frac{\Lambda_v}{\mu_v} - \bar{I}_v\right) & \frac{\beta_{hv}}{\bar{H}}\left(\frac{\Lambda_v}{\mu_v} - \bar{I}_v\right) & \frac{\hat{\beta}_{hv}}{\bar{H}}\left(\frac{\Lambda_v}{\mu_v} - \bar{I}_v\right) & -\frac{\beta_{hv}\bar{I}_h}{\bar{H}} - \frac{\hat{\beta}_{hv}\bar{R}_h}{\bar{H}} - \mu_v \end{bmatrix} $$

The application of the Routh-Hurwitz criterion proves to be difficult and does not allow to give conditions that can easily be written in function of the system parameters. We will therefore proceed differently. A complex number $w \in \mathbb{C}$ is an eigenvalue of $J_{EE}$ if and only if there exists a vector $(Z_1, Z_2, Z_3, Z_4)^T \ne (0, 0, 0, 0)$ such that the following equations are satisfied:

$$ wZ_1 = -\delta Z_2 - Z_1 \mu_h \qquad (9a) $$

$$ wZ_2 = Z_1 \left( -\frac{\beta_{vh}\bar{I}_v \bar{S}_h}{\bar{H}^2} + \frac{\beta_{vh}\bar{I}_v}{\bar{H}} \right) + \left( -\frac{\beta_{vh}\bar{I}_v}{\bar{H}} - \epsilon \right) Z_2 - \frac{\beta_{vh}\bar{I}_v Z_3}{\bar{H}} + \frac{\beta_{vh}\bar{S}_h Z_4}{\bar{H}} \qquad (9b) $$

$$ wZ_3 = \gamma_h Z_2 + (-\rho - \mu_h) Z_3 \qquad (9c) $$

$$ wZ_4 = \left[ Z_1 \left( -\frac{\beta_{hv}\bar{I}_h}{\bar{H}^2} - \frac{\hat{\beta}_{hv}\bar{R}_h}{\bar{H}^2} \right) + \frac{\beta_{hv}Z_2}{\bar{H}} + \frac{\hat{\beta}_{hv}Z_3}{\bar{H}} \right] \left( \frac{\Lambda_v}{\mu_v} - \bar{I}_v \right) - \left( \frac{\beta_{hv}\bar{I}_h}{\bar{H}} + \frac{\hat{\beta}_{hv}\bar{R}_h}{\bar{H}} + \mu_v \right) Z_4 \qquad (9d) $$

Combining the first three equations we obtain

$$ wZ_2 = \left[ \frac{\delta_h}{w + \mu_h} \left( \frac{\beta_{vh}\bar{I}_v \bar{S}_h}{\bar{H}^2} - \frac{\beta_{vh}\bar{I}_v}{\bar{H}} \right) - \frac{\beta_{vh}\bar{I}_v}{\bar{H}} - \epsilon - \frac{\beta_{vh}\bar{I}_v}{\bar{H}} \frac{\gamma_h}{w + \theta_h} \right] Z_2 + \frac{\beta_{vh}\bar{S}_h}{\bar{H}} Z_4 \\ wZ_2 = - \left[ \frac{\beta_{vh}\bar{I}_v}{\bar{H}} \left( \frac{\delta_h}{w + \mu_h} \left( 1 - \frac{\bar{S}_h}{\bar{H}} \right) + 1 + \frac{\gamma_h}{w + \theta_h} \right) + \epsilon \right] Z_2 + \frac{\beta_{vh}\bar{S}_h}{\bar{H}} Z_4 $$
This can be written as follows

$$ (1 + G_2(w))Z_2 = \frac{\beta_{vh}\bar{S}}{\epsilon_h \bar{H}} Z_4 = \frac{\bar{I}_h}{\bar{I}_v} Z_4 \text{ with } G_2(w) = \frac{w + \frac{\beta_{vh}\bar{I}_v}{\bar{H}} \left( \frac{\delta_h}{w+\mu_h} \left(1 - \frac{\bar{S}}{\bar{H}}\right) + 1 + \frac{\gamma_h}{w+\theta_h} \right)}{\epsilon_h} \quad (10) $$

Using equilibrium relations, it is possible to write

$$ G_2(w) = \frac{w + \frac{\beta_{vh}\bar{I}_v}{\bar{H}} \left( \frac{\delta_h}{w+\mu_h} \left(1 + \frac{\gamma_h}{\theta_h}\right) \frac{\bar{I}_h}{\bar{H}} + 1 + \frac{\gamma_h}{w+\theta_h} \right)}{\epsilon_h} $$

On the other hand, using (9c) together with the expressions of $Z_1$ and $Z_2$ and the equilibrium relations, we have the following successive equalities:

$$ \begin{aligned} wZ_4 &= \left[ \frac{\delta_h}{w + \mu_h} Z_2 \left( \frac{\beta_{hv}\bar{I}_h}{\bar{H}^2} + \frac{\hat{\beta}_{hv}\bar{R}_h}{\bar{H}^2} \right) + \frac{\beta_{hv}Z_2}{\bar{H}} + \frac{\hat{\beta}_{hv}}{\bar{H}} \frac{\gamma_h}{w + \theta_h} Z_2 \right] \bar{S}_v - \left( \frac{\beta_{hv}\bar{I}_h}{\bar{H}} + \frac{\hat{\beta}_{hv}\bar{R}_h}{\bar{H}} + \mu_v \right) Z_4 \\ wZ_4 &= \left( \frac{\delta_h}{w + \mu_h} \left( \frac{\beta_{hv}\bar{I}_h}{\bar{H}} + \frac{\hat{\beta}_{hv}\bar{R}_h}{\bar{H}} \right) + \beta_{hv} + \frac{\hat{\beta}_{hv}\gamma_h}{w + \theta_h} \right) \frac{\bar{S}_v}{\bar{H}} Z_2 - \left( \frac{\beta_{hv}\bar{I}_h}{\bar{H}} + \frac{\hat{\beta}_{hv}\bar{R}_h}{\bar{H}} + \mu_v \right) Z_4 \\ wZ_4 &= \left( \frac{\delta_h}{w + \mu_h} \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \beta_{hv} + \frac{\hat{\beta}_{hv}\gamma_h}{w + \theta_h} \right) \frac{\bar{S}_v}{\bar{H}} Z_2 - \left( \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \mu_v \right) Z_4 \\ \left( w + \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \mu_v \right) Z_4 &= \left( \frac{\delta_h}{w + \mu_h} \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \beta_{hv} + \frac{\hat{\beta}_{hv}\gamma_h}{w + \theta_h} \right) \frac{\bar{S}_v}{\bar{H}} Z_2 \\ &\qquad \frac{w + \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \mu_v}{\frac{\delta_h}{w + \mu_h} \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \beta_{hv} + \frac{\hat{\beta}_{hv}\gamma_h}{w + \theta_h}} Z_4 = \frac{\bar{S}_v}{\bar{H}} Z_2 \end{aligned} $$

Therefore we obtain the following relation

$$ \underbrace{\frac{\left( w + \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \mu_v \right) \frac{\bar{I}_v \bar{H}}{\bar{S}_v \bar{I}_h}}{\frac{\delta_h}{w + \mu_h} \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \beta_{hv} + \frac{\hat{\beta}_{hv} \gamma_h}{w + \theta_h}}}_{f(w)} Z_4 = \frac{\bar{I}_v}{\bar{I}_h} Z_2 \quad (11) $$

We have

$$ \begin{aligned} & \left| \frac{\delta_h}{w + \mu_h} \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \beta_{hv} + \frac{\hat{\beta}_{hv} \gamma_h}{w + \theta_h} \right| \\ & \le \frac{\delta_h}{|w + \mu_h|} \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \beta_{hv} + \hat{\beta}_{hv} \gamma_h \frac{1}{|w + \theta_h|} \\ & \le \frac{\delta_h}{\mathrm{Re}(w) + \mu_h} \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \beta_{hv} + \hat{\beta}_{hv} \gamma_h \frac{1}{\mathrm{Re}(w) + \theta_h} \\ & \le \frac{\delta_h}{\mu_h} \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \beta_{hv} + \hat{\beta}_{hv} \gamma_h \frac{1}{\theta_h} \quad \text{if } \mathrm{Re}(w) \ge 0 \\ & = \left( \frac{\delta_h \bar{I}_h}{\mu_h \bar{H}} + 1 \right) \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \end{aligned} $$
So, if $Re(w) \ge 0$

$$ 
\begin{aligned}
|f(w)| & \ge \frac{\left( \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right) \frac{\bar{I}_h}{\bar{H}} + \mu_v \right)}{\left( \frac{\delta_h \bar{I}_h}{\mu_h \bar{H}} + 1 \right) \left( \beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h} \right)} \frac{\bar{I}_v \bar{H}}{\bar{S}_v \bar{I}_h} \\
& = \frac{\frac{\beta_{hv} + \frac{\hat{\beta}_{hv} \gamma_h}{\theta_h}}{\mu_v} \frac{\bar{I}_h}{\bar{H}} + 1}{\left( \frac{\delta_h \bar{I}_h}{\mu_h \bar{H}} + 1 \right) \frac{\beta_{hv} + \frac{\hat{\beta}_{hv} \gamma_h}{\theta_h}}{\mu_v}} \frac{\bar{I}_v \bar{H}}{\bar{S}_v \bar{I}_h}
\end{aligned}
$$

We have $\frac{\bar{I}_v}{\bar{S}_v} = \frac{\beta_{hv} + \frac{\hat{\beta}_{hv} \gamma_h}{\theta_h}}{\mu_v} \frac{\bar{I}_h}{\bar{H}}$. So, $|f(w)| \ge \frac{\frac{\beta_{hv} + \frac{\hat{\beta}_{hv} \gamma_h}{\theta_h}}{\mu_v} \frac{\bar{I}_h}{\bar{H}} + 1}{\frac{\delta_h \bar{I}_h}{\mu_h \bar{H}} + 1}$.

Therefore, we have that

If $\frac{\beta_{hv} + \frac{\hat{\beta}_{hv} \gamma_h}{\theta_h}}{\mu_v} > \frac{\delta_h}{\mu_h}$ then $Re(w) \ge 0 \Rightarrow |f(w)| > 1$.

We can write the relations (10-11) in a matrix form as follows

$$ \mathrm{diag}(1 + G_2(w), f(w)) Z = M Z \qquad (12) $$

The matrix $M$ is given by

$$ M = \begin{pmatrix} 0 & \frac{\bar{I}_h}{\bar{I}_v} \\ \frac{\bar{I}_v}{\bar{I}_h} & 0 \end{pmatrix} $$

It should be noted that the matrix $H$ has non-negative entries and that the equilibrium $\bar{x} = (\bar{I}_h, \bar{I}_v)$ satisfies $M\bar{x} = \bar{x}$.

Suppose that $Re(w) \ge 0$. We shall use a Krasnoselskii trick ([17], Proof of Theorem 6.1, see also [13]) to reach a contradiction. From (12), we get using modulus and denoting $|Z| = (|Z_2|, |Z_4|)^T$:

$$ \min(1 + \mathrm{Re}G_2(w), |f(w)|) |Z| \le M|Z|. \qquad (13) $$

Since the vector $\bar{x}$ is positive, there exists a minimal positive real number $r$ satisfying

$$ |Z| \le r\bar{x}. \qquad (14) $$

Let $\eta(w) = \min(1 + \mathrm{Re}G_2(w), |f(w)|)$. We have proved that $Re(w) \ge 0$ implies that $|f(w)| > 1$ and $\mathrm{Re}G_2(w) > 0$ and hence, $\eta(w) > 1$. Combining the relations (13) and (14), we obtain

$$ \eta(w) |Z| \le M|Z| \le M r \bar{x} = r \bar{x}. $$

The later contradict the minimality of $r$.

**Remark 3.** Numerical simulations show that the sufficient condition (8) is not a necessary one. We think that $R_0 > 1$ (without any supplementary condition) should be sufficient for the local stability of the endemic equilibrium but we could not prove it.
## 2.5 No disease induced death $\delta_h = 0$

When there is no disease induced death: $\delta_h = 0$, System (I) is written as follows

$$ 
\left\{
\begin{aligned}
\frac{dS_h}{dt} &= \Lambda_h + \rho_h R_h + \alpha_h I_h - \beta_{vh} \frac{I_v}{H} S_h - \mu_h S_h, \\
\frac{dI_h}{dt} &= \beta_{vh} \frac{I_v}{H} S_h - \tilde{\epsilon}_h I_h, \quad \text{with } \tilde{\epsilon}_h = \mu_h + \alpha_h + \gamma_h, \\
\frac{dR_h}{dt} &= \gamma_h I_h - \theta_h R_h, \\
\frac{dS_v}{dt} &= \Lambda_v - \left(\beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H}\right) S_v - \mu_v S_v, \\
\frac{dI_v}{dt} &= \left(\beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H}\right) S_v - \mu_v I_v.
\end{aligned}
\right.
\quad (15) $$

Since $\delta_h = 0$, the equations of the total population of mosquitoes and that of human are given as follows:

$$ 
\begin{aligned}
\frac{dH}{dt} &= \Lambda_h - \mu_h H, \\
\frac{dV}{dt} &= \Lambda_v - \mu_v V.
\end{aligned}
 $$

This last system admits a unique equilibrium point $(H^*, V^*) = (\frac{\Lambda_h}{\mu_h}, \frac{\Lambda_v}{\mu_v})$ which is GAS, in this case we can use the [26, Theorem 3.1], we can reduce the dimension of System (15) to a system of dimension 3 as follows:

$$ 
\left\{
\begin{aligned}
\frac{dI_h}{dt} &= \beta_{vh} \frac{I_v}{H^*} (H^* - I_h - R_h) - \tilde{\epsilon}_h I_h, \\
\frac{dR_h}{dt} &= \gamma_h I_h - \theta_h R_h, \\
\frac{dI_v}{dt} &= \left( \beta_{hv} \frac{I_h}{H^*} + \hat{\beta}_{hv} \frac{R_h}{H^*} \right) (V^* - I_v) - \mu_v I_v.
\end{aligned}
\right.
\quad (16) $$

### 2.5.1 The endemic equilibrium

The endemic equilibrium (EE) of the System (16) is solution of the following equations:

$$ 
\begin{aligned}
\beta_{vh} \frac{I_v}{H^*} (H^* - I_h - R_h) - \tilde{\epsilon}_h I_h &= 0 \\
\gamma_h I_h - \theta_h R_h &= 0 \\
\left( \beta_{hv} \frac{I_h}{H^*} + \hat{\beta}_{hv} \frac{R_h}{H^*} \right) (V^* - I_v) - \mu_v I_v &= 0
\end{aligned}
 $$

Then the EE is given by:

$$ 
\begin{aligned}
\bar{I}_h &= \frac{\tilde{\epsilon}_h \mu_v H^* (\mathcal{R}_0^2 - 1)}{\left(\beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h}\right) \left(\beta_{vh} \frac{V^*}{H^*} \left(1 + \frac{\theta_h}{\gamma_h}\right) + \tilde{\epsilon}_h\right)}, \\
\bar{R}_h &= \frac{\gamma_h}{\theta_h} \bar{I}_h, \\
\bar{I}_v &= \frac{\left(\beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h}\right) \bar{I}_h V^*}{\mu_v H^* + \left(\beta_{hv} + \hat{\beta}_{hv} \frac{\gamma_h}{\theta_h}\right) \bar{I}_h}.
\end{aligned}
 $$
The Jacobian matrix at the endemic equilibrium is given by:

$$ J_0(\bar{I}_h, \bar{R}_h, \bar{I}_v) = \begin{pmatrix} -\beta_{vh} \frac{\bar{I}_v}{H^*} - \tilde{\epsilon}_h & -\beta_{vh} \frac{\bar{I}_v}{H^*} & \beta_{vh} \frac{1}{H^*} (H^* - \bar{I}_h - \bar{R}_h) \\ \gamma_h & -\theta_h & 0 \\ \beta_{hv} \frac{1}{H^*} (V^* - \bar{I}_v) & \hat{\beta}_{hv} \frac{1}{H^*} (V^* - \bar{I}_v) & -\left(\beta_{hv} \frac{\bar{I}_h}{H^*} + \hat{\beta}_{hv} \frac{\bar{R}_h}{H^*}\right) - \mu_v \end{pmatrix}. $$

Using the equilibrium equations we can write the Jacobian as follow:

$$ J_0 = \begin{pmatrix} -\beta_{vh} \frac{\bar{I}_v}{H^*} - \tilde{\epsilon}_h & -\beta_{vh} \frac{\bar{I}_v}{H^*} & \tilde{\epsilon}_h \frac{\bar{I}_h}{\bar{I}_v} \\ \gamma_h & -\theta_h & 0 \\ \beta_{hv} \frac{1}{H^*} (V^* - \bar{I}_v) & \hat{\beta}_{hv} \frac{1}{H^*} (V^* - \bar{I}_v) & -\left(\beta_{hv} \frac{\bar{I}_h}{H^*} + \hat{\beta}_{hv} \frac{\bar{R}_h}{H^*}\right) \frac{V^*}{\bar{I}_v} \end{pmatrix}. $$

The second additive matrix at the endemic equilibrium is given as follow:

$$ J_0^{[2]} = \begin{pmatrix} -\beta_{vh} \frac{\bar{I}_v}{H^*} - \tilde{\epsilon}_h - \theta_h & 0 & -\tilde{\epsilon}_h \frac{\bar{I}_h}{\bar{I}_v} \\ \frac{\hat{\beta}_{hv}}{H^*} (V^* - \bar{I}_v) & -\beta_{vh} \frac{\bar{I}_v}{H^*} - \tilde{\epsilon}_h - \left(\beta_{hv} \frac{\bar{I}_h}{H^*} + \hat{\beta}_{hv} \frac{\bar{R}_h}{H^*}\right) \frac{V^*}{\bar{I}_v} & -\beta_{vh} \frac{\bar{I}_v}{H^*} \\ \frac{\beta_{hv}}{H^*} (V^* - \bar{I}_v) & \gamma_h & -\left(\beta_{hv} \frac{\bar{I}_h}{H^*} + \hat{\beta}_{hv} \frac{\bar{R}_h}{H^*}\right) \frac{V^*}{\bar{I}_v} - \theta_h \end{pmatrix} \quad (17) $$

**Proposition 4.** The endemic equilibrium exists and is locally asymptotically stable if and only if $R_0 > 1$.

*Proof.*

$$ \begin{aligned} \det(J_0) &= -\theta_h \frac{V^*}{\bar{I}_v} \left(\beta_{vh} \frac{\bar{I}_v}{H^*} + \tilde{\epsilon}_h\right) \left(\beta_{hv} \frac{\bar{I}_h}{H^*} + \hat{\beta}_{hv} \frac{\bar{R}_h}{H^*}\right) - \beta_{vh} \gamma_h \frac{V^*}{H^*} \left(\beta_{hv} \frac{\bar{I}_h}{H^*} + \hat{\beta}_{hv} \frac{\bar{R}_h}{H^*}\right) \\ &\quad + \tilde{\epsilon}_h \frac{\bar{I}_h}{\bar{I}_v} \left[\gamma_h \hat{\beta}_{hv} \frac{1}{H^*} (V^* - \bar{I}_v) + \theta_h \beta_{hv} \frac{1}{H^*} (V^* - \bar{I}_v)\right] \\ &= -(\theta_h + \gamma_h) \left[\beta_{vh} \frac{V^*}{H^*} \left(\beta_{hv} \frac{\bar{I}_h}{H^*} + \hat{\beta}_{hv} \frac{\bar{R}_h}{H^*}\right) + \tilde{\epsilon}_h \bar{I}_h \frac{\beta_{hv}}{H^*}\right] \end{aligned} $$

$$ \begin{aligned} \det(J_0^{[2]}) = & -\left(\beta_{vh} \frac{\bar{I}_v}{H^*} + \tilde{\epsilon}_h + \theta_h\right) \\ & \left[ \left(\beta_{vh} \frac{\bar{I}_v}{H^*} + \tilde{\epsilon}_h + \left(\beta_{hv} \frac{\bar{I}_h}{H^*} + \hat{\beta}_{hv} \frac{\bar{R}_h}{H^*}\right) \frac{V^*}{\bar{I}_v}\right) \left(\left(\beta_{hv} \frac{\bar{I}_h}{H^*} + \hat{\beta}_{hv} \frac{\bar{R}_h}{H^*}\right) \frac{V^*}{\bar{I}_v} + \theta_h\right) \right] \\ & - \left(\beta_{vh} \frac{\bar{I}_v}{H^*} + \tilde{\epsilon}_h + \theta_h\right) \beta_{vh} \frac{\bar{I}_v}{H^*} \gamma_h - \tilde{\epsilon}_h \frac{\bar{I}_h}{\bar{I}_v} \frac{\hat{\beta}_{hv}}{H^*} (V^* - \bar{I}_v) \gamma_h \\ & - \tilde{\epsilon}_h \frac{\bar{I}_h}{\bar{I}_v} \frac{\beta_{hv}}{H^*} (V^* - \bar{I}_v) \left(\beta_{vh} \frac{\bar{I}_v}{H^*} + \tilde{\epsilon}_h + \left(\beta_{hv} \frac{\bar{I}_h}{H^*} + \hat{\beta}_{hv} \frac{\bar{R}_h}{H^*}\right) \frac{V^*}{\bar{I}_v}\right) \end{aligned} $$
$\det(J_0)$, $tr(J_0)$ and $\det(J_0^{[2]})$ are all negative then using [21] (Lemma 3) we have that the endemic equilibrium is locally asymptotically stable. □

### 2.5.2 The global stability of the EE

For the global stability of the EE we have the following result.

**Proposition 5.** If $\rho_h - \alpha_h \le \frac{\gamma_h}{2}$ then the EE is globally asymptotically stable if $\mathcal{R}_0 > 1$.

*Proof.* To prove the global stability of the endemic equilibrium we shall use the geometric approach introduced by Li and Muldowney [18]. The second additive matrix at a point $(I_h, R_h, I_v)$ is given as follows:

$$J^{[2]} = \begin{pmatrix} -\beta_{vh} \frac{I_v}{H^*} - \tilde{\epsilon}_h - \theta_h & 0 & -\frac{\beta_{vh}}{H^*} (H^* - I_h - R_h) \\ \frac{\hat{\beta}_{hv}}{H^*} (V^* - I_v) & -\frac{\beta_{vh}}{H^*} (I_v + I_h) - \hat{\beta}_{hv} \frac{R_h}{H^*} - \tilde{\epsilon}_h - \mu_v & -\beta_{vh} \frac{I_v}{H^*} \\ -\frac{\beta_{hv}}{H^*} (V^* - I_v) & \gamma_h & -\left(\beta_{hv} \frac{I_h}{H^*} + \hat{\beta}_{hv} \frac{R_h}{H^*}\right) - \theta_h - \mu_v \end{pmatrix}$$

Consider the following matrix $P$ which is nonsingular in the interior of $\Omega$:

$$P = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{I_h}{I_v} & 0 \\ 0 & \frac{I_h}{I_v} & \frac{I_h}{I_v} \end{pmatrix}$$

With the same notations as in [18], let $P_f = (DP)(f)$, where $f$ is the vector field of (16), we have

$$P_f = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \frac{\dot{I}_h I_v - I_h \dot{I}_v}{I_v^2} & 0 \\ 0 & \frac{\dot{I}_h I_v - I_h \dot{I}_v}{I_v^2} & \frac{\dot{I}_h I_v - I_h \dot{I}_v}{I_v^2} \end{pmatrix}, \text{ and } P_f P^{-1} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \frac{\dot{I}_h}{I_h} - \frac{\dot{I}_v}{I_v} & 0 \\ 0 & 0 & \frac{\dot{I}_h}{I_h} - \frac{\dot{I}_v}{I_v} \end{pmatrix}$$

$$\text{Let } B = P_f P^{-1} + P J^{[2]} P^{-1} = \begin{pmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{pmatrix} \text{ where, } B_{11} = -\beta_{vh} \frac{I_v}{H^*} - \tilde{\epsilon}_h - \theta_h,$$

$$B_{12} = \left( \frac{\beta_{vh}}{H^*} S_h \frac{I_v}{I_h}, -\frac{\beta_{vh}}{H^*} S_h \frac{I_v}{I_h} \right), B_{21} = \left( \frac{\hat{\beta}_{hv}}{H^*} (V^* - I_v) \frac{I_h}{I_v} \frac{\beta_{hv} - \beta_{hv}}{H^*} (V^* - I_v) \frac{I_h}{I_v} \right)$$

$$B_{22} = \begin{pmatrix} \frac{\dot{I}_h}{I_h} - \frac{\dot{I}_v}{I_v} - \beta_{hv} \frac{I_h}{H^*} - \hat{\beta}_{hv} \frac{R_h}{H^*} - \tilde{\epsilon}_h - \mu_v & -\beta_{vh} \frac{I_v}{H^*} \\ \rho_h - \alpha_h & \frac{\dot{I}_h}{I_h} - \frac{\dot{I}_v}{I_v} - \beta_{vh} \frac{I_v}{H^*} - \beta_{hv} \frac{I_h}{H^*} - \hat{\beta}_{hv} \frac{R_h}{H^*} - \mu_v - \theta_h \end{pmatrix}$$
The vector norm |·| is chosen as

$$|(u, v, w)| = \sup \{|u|, |v| + |w|\}$$

The Lozinskii measure $\mu(B)$ with respect to |·| can be estimated as follows [18]

$$\mu(B) \leq \sup\{g_1, g_2\}.$$

Where,

$$g_1 = B_{11} + |B_{12}|$$

$$g_2 = |B_{21}| + \mu_1(B_{22})$$

We have: $B_{11} = -\beta_{vh} \frac{I_v}{H^*} - \tilde{\epsilon}_h - \theta_h$, and $|B_{12}| = \frac{\beta_{vh}}{H^*} S_h \frac{I_v}{I_h}$. Thus,

$$g_1 = B_{11} + |B_{12}| = -\beta_{vh} \frac{I_v}{H^*} - \tilde{\epsilon}_h - \theta_h + \frac{\beta_{vh}}{H^*} S_h \frac{I_v}{I_h}.$$

Using the fact that $\frac{\dot{I}_h}{I_h} = \frac{\beta_{vh}}{H^*} S_h \frac{I_v}{I_h} - \tilde{\epsilon}_h$, we obtain,

$$g_1 = -\beta_{vh} \frac{I_v}{H^*} - \theta_h + \frac{\dot{I}_h}{I_h} \leq \frac{\dot{I}_h}{I_h} - \theta_h$$

To estimate $g_2$, we have

$$|B_{21}| = \frac{(V^* - I_v)}{H^*} \frac{I_h}{I_v} (|\hat{\beta}_{hv}| + |\hat{\beta}_{hv} - \beta_{hv}|) = \beta_{hv} \frac{(V^* - I_v)}{H^*} \frac{I_h}{I_v} \text{ since } \hat{\beta}_{hv} < \beta_{hv}$$

Using the fact that $\frac{dI_v}{dt} = (\beta_{hv} \frac{I_h}{H^*} + \hat{\beta}_{hv} \frac{R_h}{H^*}) (V^* - I_v) - \mu_v I_v$, we can write,

$$|B_{21}| = \frac{\dot{I}_v}{I_v} - \hat{\beta}_{hv} \frac{R_h}{H^* I_v} (V^* - I_v) + \mu_v$$

$$\mu_1(B_{22}) = \max\{\frac{\dot{I}_h}{I_h} - \frac{\dot{I}_v}{I_v} - \beta_{hv} \frac{I_h}{H^*} - \hat{\beta}_{hv} \frac{R_h}{H^*} - \tilde{\epsilon}_h - \mu_v + |\rho_h - \alpha_h|, \frac{\dot{I}_h}{I_h} - \frac{\dot{I}_v}{I_v} - \beta_{hv} \frac{I_h}{H^*} - \hat{\beta}_{hv} \frac{R_h}{H^*} - \mu_v - \theta_h\}$$

Since $\tilde{\epsilon}_h = \alpha_h + \gamma_h + \mu_h$, and $\theta_h = \rho_h + \mu_h$, we have

$$\mu_1(B_{22}) \leq \frac{\dot{I}_h}{I_h} - \frac{\dot{I}_v}{I_v} - \beta_{hv} \frac{I_h}{H^*} - \hat{\beta}_{hv} \frac{R_h}{H^*} - \mu_v - \mu_h + \max\{-\alpha_h - \gamma_h + |\rho_h - \alpha_h|, -\rho_h\}$$

Let $M$ be the number defined by $M = \max\{-\alpha_h - \gamma_h + |\rho_h - \alpha_h|, -\rho_h\}$. Then we have

$$g_2 = |B_{21}| + \mu_1(B_{22}) \leq \frac{\dot{I}_h}{I_h} - \mu_h + M$$

1. If $\rho_h \leq \alpha_h$ then $|\rho_h - \alpha_h| = -\rho_h + \alpha_h$, and so $M = \max\{-\gamma_h - \rho_h, -\rho_h\} = -\rho_h$. Then,

$$g_2 \leq \frac{\dot{I}_h}{I_h} - \mu_h - \rho_h = \frac{\dot{I}_h}{I_h} - \theta_h$$

and therefore,

$$\mu(B) \leq \sup\{g_1, g_2\} \leq \sup\{\frac{\dot{I}_h}{I_h} - \theta_h, \frac{\dot{I}_h}{I_h} - \theta_h\} = \frac{\dot{I}_h}{I_h} - \theta_h$$
2. If $\rho_h > \alpha_h$ then $|\rho_h - \alpha_h| = \rho_h - \alpha_h$ and $M = \max\{\rho_h - 2\alpha_h - \gamma_h, -\rho_h\}$. Since, by assumption, $\rho_h - \alpha_h \le \frac{\gamma_h}{2}$ we get $M \le -\rho_h$, and hence

$$g_2 \le \frac{\dot{I}_h}{I_h} - \mu_h - \rho_h = \frac{\dot{I}_h}{I_h} - \theta_h.$$

Finally:

$$\mu(B) \le \sup\{g_1, g_2\} \le \sup\{\frac{\dot{I}_h}{I_h} - \theta_h, \frac{\dot{I}_h}{I_h} - \theta_h\} = \frac{\dot{I}_h}{I_h} - \theta_h$$

So both cases lead to the following inequality,

$$\mu(B) \le \frac{\dot{I}_h}{I_h} - \theta_h.$$

System (16) is uniformly persistence. So $\exists r \ge 0$ such that $I_h(t) > r$, $R_h(t) > r$ and $I_v(t) > r$. This also ensures the existence of a compact set $K$ which is absorbing in the interior of $\Omega$.

The uniform persistence constant $r$ can be adjusted so that there exists $\bar{t} > 0$ independent of the initial condition in $\Omega$, such that

$I_h(t) > r, R_h(t) > r$ and $I_v(t) > r$ for $t > \bar{t}$. This lead us to

$$\begin{aligned} \frac{1}{t} \int_{0}^{t} \mu(B) ds &= \frac{1}{t} \int_{0}^{\bar{t}} \mu(B) ds + \frac{1}{t} \int_{\bar{t}}^{t} \mu(B) ds \\ &\le \frac{1}{t} \int_{0}^{\bar{t}} \mu(B) ds + \frac{1}{t} \int_{\bar{t}}^{t} \left( \frac{\dot{I}_h}{I_h} - \theta_h \right) ds \end{aligned}$$

$$\frac{1}{t} \int_{0}^{t} \mu(B) ds \le \frac{1}{t} \int_{0}^{\bar{t}} \mu(B) ds + \frac{\ln(I_h(t))}{\ln(I_h(\bar{t}))} - \theta_h(1 - \frac{\bar{t}}{t}) \quad (18)$$

Define $\bar{q}_2 = \limsup_{t \to \infty} \sup_{(S(0),I(0)) \in K} \frac{1}{t} \int_0^t \mu(B) ds$. We have

$$\begin{aligned} \sup_{(S(0),I(0)) \in K} \frac{1}{t} \int_0^t \mu(B) ds &\le \sup_{(S(0),I(0)) \in K} \left( \frac{1}{t} \int_0^{\bar{t}} \mu(B) ds + \frac{1}{t} \frac{\ln(I_h(t))}{\ln(I_h(\bar{t}))} - \theta_h(1 - \frac{\bar{t}}{t}) \right) \\ &\le \sup_{(S(0),I(0)) \in K} \frac{1}{t} \int_0^{\bar{t}} \mu(B) ds + \frac{1}{t} \frac{\ln(H^*)}{\ln(I_h(\bar{t}))} - \theta_h(1 - \frac{\bar{t}}{t}). \end{aligned}$$

This implies

$$\bar{q}_2 \le -\theta_h.$$

Theorem 3.5 in [18] allows then to conclude that the endemic equilibrium EE is globally asymptotically stable in $\bar{\Omega}$. $\square$

## 3 State Estimation

We are interested here in estimating the size of the populations $S_h$, $I_h$ and $R_h$. Assuming that (I) (or, equivalently, System (4)) is a “quite good” model of the system under consideration. If it is possible to have the value of the state $(H, I_h, R_h, I_v)$ at some time $t_0$ then it is possible to compute $(H(t), I_h(t), R_h(t), I_v(t))$ for all $t \ge t_0$ by integrating the differential system (4) with the
initial condition $(H(t_0), I_h(t_0), R_h(t_0), I_v(t_0))$. Unfortunately, it is often not possible to measure the whole state at a given time and therefore it is not possible to integrate the differential equation because one does not know an initial condition. In practice, the only measurement available is the host incidence, that is the number of new infected humans per unit time. This information is usually accessible to the Public Health Services. This available measurement will be denoted $y(t)$. In the considered model, it corresponds to the term $\beta_{vh} \frac{I_v(t)}{H} S_h(t)$. Our goal is to use the information (that we assume to be continuously available) $y(t) = \beta_{vh} \frac{I_v(t)}{H} S_h(t) = \beta_{vh} \frac{I_v(t)}{H(t)} (H(t) - I_h(t) - R_h(t))$ together with the model (I) (or, equivalently, System (4)) in order to obtain dynamical estimates $\hat{H}(t)$, $\hat{I}_h(t)$ and $\hat{R}_h(t)$ of $H(t)$, $I_h(t)$ and $R_h(t)$. A solution to this estimation problem can be provided by a tool developed in automatic control theory: the use of so-called state observers or state estimators. An observer is an auxiliary dynamical system $\hat{\Sigma}$ designed to provide dynamical estimates of the complete state of another system $\Sigma$ —in this case the epidemiological model of interest—using the available information given by partial measurements $y(t)$ of the state of $\Sigma$. The solutions of this auxiliary dynamical system must converge (as fast as possible) towards the solutions of the original system. More precisely, an exponential observer (or state estimator) for (4)) is a dynamical system

$$ \left\{ \begin{aligned} \frac{d\hat{z}}{dt} &= \hat{F}(\hat{z}(t), y(t)), \\ \hat{x}(t) &= G(\hat{z}(t), y(t)), \end{aligned} \right. \qquad (19) $$

whose solutions $\hat{x}(t)$ converge exponentially to the solutions $x(t)$ of system (4), i.e., there exists $\lambda > 0$ such that, for all $t \ge 0$ and for all initial conditions $x(0), \hat{z}(0)$, the corresponding solutions of (4)–(19) satisfy

$$ \|\hat{x}(t) - x(t)\| \le \exp(-\lambda t) \|\hat{x}(0) - x(0)\|. $$

It must be pointed out that the initial state $x_0$ of (4) is unknown while the initial state $z_0$ of the observer (19) can be chosen arbitrarily. The problem to address is to construct the good vector field $\hat{F}$ and the map $G$ in such a way that the above condition is satisfied. Hereafter, we shall give two possible constructions.

## 3.1 A simple observer

We focus on the "human" part of model (I) and we write it introducing $y(t)$ and using the variable $H$ instead of $S$:

$$ \left\{ \begin{aligned} \frac{dH}{dt} &= \Lambda_h - \mu_h H - \delta_h I_h, \\ \frac{dI_h}{dt} &= \beta_{vh} \frac{I_v}{H} (H - I_h - R_h) - (\alpha_h + \gamma_h + \mu_h + \delta_h) I_h, \\ \frac{dR_h}{dt} &= \gamma_h I_h - (\rho_h + \mu_h) R_h, \\ y &= \beta_{vh} \frac{I_v}{H} (H - I_h - R_h) \end{aligned} \right. \qquad (20) $$
A simple observer for system (20) is given by:

$$ 
\left\{
\begin{aligned}
\frac{d\hat{H}}{dt} &= \Lambda_h - \mu_h \hat{H} - \delta_h \hat{I}_h, \\
\frac{d\hat{I}_h}{dt} &= y - (\alpha_h + \gamma_h + \mu_h + \delta_h) \hat{I}_h, \\
\frac{d\hat{R}_h}{dt} &= \gamma_h \hat{I}_h - (\rho_h + \mu_h) \hat{R}_h.
\end{aligned}
\right.
\qquad (21) $$

The estimation error $e(t) = \hat{x}(t) - x(t)$ satisfies differential equation $\dot{e} = Ae$ with

$$ A = \begin{pmatrix} -\mu_h & -\delta_h & 0 \\ 0 & -(\alpha_h + \gamma_h + \mu_h + \delta_h) & 0 \\ 0 & \gamma_h & -(\rho_h + \mu_h) \end{pmatrix}. $$

The matrix A is Hurwitz which implies that the estimation error $e(t)$ converges exponentially fast to zero.

It should be noticed that the observer (21) provides estimates $\hat{H}(t)$, $\hat{I}_h(t)$ and $\hat{R}_h(t)$ of the state without using the values of the various $\beta$ that are in general not well known. If $\beta_{vh}$ is known then we also have an estimate of $I_v(t)$ given by

$$ \hat{I}_v(t) = \frac{\hat{H}(t)}{\beta_{vh}(\hat{H}(t) - \hat{I}_h(t) - \hat{R}_h(t))} y(t). $$

The weakness of the state estimator (21) is that its convergence speed cannot be adjusted. Figures 2-3 show the convergence of the estimates delivered by the estimator (21) towards the unmeasured states solutions of (4). It must be noticed that the estimator is good for estimating the size of infected humans $I_h(t)$ and $R_h(t)$ but it is not that good for estimating the total human population $H(t)$: convergence of $\hat{H}(t)$ towards $H(t)$ is rather slow due to the fact that the natural mortality rate $\mu_h$ is very small.

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: The evolution of $R_h(t)$ (green curve), solution of (4), and its estimate $\hat{R}_h(t)$ (red curve) delivered by the observer (21).</figcaption></figure>

## 3.2 A faster estimator

If, in addition to the host incidence, the total population $H(t)$ can be measured then it is possible to construct another estimator to dynamically estimate $I_h(t)$. The convergence of this estimator is faster than the one of the estimator given by (21). Now System (20) has two outputs: $y(t) = \beta_{vh} \frac{I_v(t)}{H(t)}(H(t) - I_h(t) - R_h(t))$ and $y_1(t) = H(t)$. An estimate $\hat{x} = (\hat{H}, \hat{I}_h, \hat{R}_h)^T$ of the state $x = (H, I_h, R_h)^T$ can be computed thanks to the following exponential estimator given by the following Kalman deterministic observer ([6], page 16):

$$ \left\{ \begin{array}{l} \dot{\hat{z}} = F(\hat{x}, y) - \Sigma(t)C^T(C\hat{z} - y_1), \\ \dot{\Sigma} = Q_{\xi} + \Sigma A^T + A\Sigma - \Sigma C^T C\Sigma, \\ Q_{\xi} = \xi^2 \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1/\xi^2 & 0 \\ 0 & 0 & 1/\xi^4 \end{pmatrix}, C = \begin{pmatrix} 1 & 0 & 0 \end{pmatrix}, \\ A = \begin{pmatrix} -\mu_h & -\delta_h & 0 \\ 0 & -(\alpha_h + \gamma_h + \mu_h + \delta_h) & 0 \\ 0 & \gamma_h & -(\rho_h + \mu_h) \end{pmatrix}, F(\hat{x}, y) = A\hat{x} + \begin{pmatrix} \Lambda_h \\ y \\ 0 \end{pmatrix}, \end{array} \right. \quad (22) $$

The positive real number $\xi$ can be chosen to adjust the speed of the convergence of the estimate $\hat{I}_h(t)$ towards $I_h(t)$ but the convergence speed of $\hat{R}_h(t)$ towards $R_h(t)$ cannot be adjusted, it is given by $-(\rho_h + \mu_h)$.
Figure 4 show the convergence of the estimate $\hat{I}_h(t)$ delivered by the estimator (22) towards the unmeasured state $I_h(t)$ solution of (4). It can be noticed that the convergence is much faster than the one corresponding to the estimator (21).

<figure><img src="image_5.png" /><figcaption>Figure 4: The evolution of $I_h(t)$ (green curve), solution of (4), and its estimate $\hat{I}_h(t)$ (red curve) delivered by the observer (22).</figcaption></figure>

## 3.3 A high-gain estimator

If we suppose that all the parameters of model (I) are known then it is possible to built an observer whose convergence speed can be adjusted by the user in order to have a very fast convergence. This will be done using a "high-gain observer" whose construction has been developed in [12]. We consider the complete model (4). In the beginning of the epidemic, it is reasonable to assume that $S(t) = H(t) - I_h(t) - R_h(t)$ is close to $H(t)$, that is $\frac{H - I_h - R_h}{H} \approx 1$. Therefore, System (4) can be approximated by the following system:

$$ \left\{ \begin{aligned} \frac{dH}{dt} &= \Lambda_h - \mu_h H - \delta_h I_h, \\ \frac{dI_h}{dt} &= \beta_{vh} I_v - (\alpha_h + \gamma_h + \mu_h + \delta_h) I_h, \\ \frac{dR_h}{dt} &= \gamma_h I_h - (\rho_h + \mu_h) R_h, \\ \frac{dI_v}{dt} &= \left( \beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H} \right) \left( \frac{\Lambda_v}{\mu_v} - I_v \right) - \mu_v I_v. \end{aligned} \right. \quad (23) $$

The measurable output can be approximated by $y(t) = \beta_{vh}I_v(t)$. Since $\beta_{vh}$ is supposed to be known, we can suppose that the available measurable output is $y(t) = I_v(t)$. We denote $x =$
$(H, I_h, R_h, I_v)^T$ We perform a coordinates change $z = \Phi(x)$ as follows

$$z_1 = y, z_2 = \frac{dy}{dt}, z_3 = \frac{d^2y}{dt^2}, z_4 = \frac{d^3y}{dt^3}.$$

Hence, we have $z_1 = I_v$, $z_2 = (\beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H}) (\frac{\Lambda_v}{\mu_v} - I_v) - \mu_v I_v$. The explicit expressions of $z_3$ and $z_4$ are too long but are easily computable.

With the new coordinates, System (23) is given by the following simpler form:

$$\left\{
\begin{aligned}
\frac{dz_1}{dt} &= z_2, \\
\frac{dz_2}{dt} &= z_3, \\
\frac{dz_3}{dt} &= z_4, \\
\frac{dz_4}{dt} &= \psi(z_1, z_2, z_3, z_4) = \frac{d^4 y}{dt^4} (\Phi^{-1}(z)).
\end{aligned}
\right.
\qquad (24)$$

The function $\psi$ is smooth on the compact set $\Phi(\Omega)$. So it is globally Lipschitz on $\Phi(\Omega)$. We consider $\tilde{\psi}$ a Lipschitz extension of $\psi$ to the whole $\mathbb{R}^4$, i.e., $\tilde{\psi}$ is a globally Lipschitz function defined on $\mathbb{R}^4$ and satisfies $\tilde{\psi}(z) = \psi(z)$ for all $z \in \Phi(\Omega)$. Using the construction of [12], an exponential observer for System (24) is given by:

$$\left\{
\begin{aligned}
\frac{d\hat{z}_1}{dt} &= \hat{z}_2 - 4\xi(\hat{z}_1 - z_1), \\
\frac{d\hat{z}_2}{dt} &= \hat{z}_3 - 6\xi^2(\hat{z}_1 - z_1), \\
\frac{d\hat{z}_3}{dt} &= \hat{z}_4 - 4\xi^3(\hat{z}_1 - z_1), \\
\frac{d\hat{z}_4}{dt} &= \tilde{\psi}(\hat{z}_1, \hat{z}_2, \hat{z}_3, \hat{z}_4) - \xi^4(\hat{z}_1 - z_1).
\end{aligned}
\right.
\qquad (25)$$

If the positive real number $\xi$ is chosen sufficiently large than the solutions of (25) converge exponentially to the solutions of System (24). More precisely, the solutions $\hat{z}(t)$ of (25) and the solutions $z(t)$ of (24) satisfy for all initial conditions $(\hat{z}(0), z(0))$:

$$\|\hat{z}(t) - z(t)\| \le \exp(-\xi t/3) \|\hat{z}(0) - z(0)\|.$$

This shows that the user can adjust the convergence speed by choosing $\xi$. The dynamical estimates of $H(t), I_h(t), R_h(t)$ and $I_v(t)$ are given by:

$$(\hat{H}(t), \hat{I}_h(t), \hat{R}_h(t), \hat{I}_v(t)) = \tilde{\Phi}^{-1}(\hat{z}_2(t), \hat{z}_2(t), \hat{z}_3(t), \hat{z}_4(t)),$$

where $\tilde{\Phi}$ is a Lipschitz extension of $\Phi$. Figure 5 illustrate the performance of the high-gain estimator (25). One can remark that the convergence is much more faster than the convergence of the "simple" estimator (21).

# 4 Disease control tentative

## 4.1 Using constant pesticide induced death rate

This section is devoted to the study of the effect of pesticide for the control of malaria transmission, we introduce an other death rate for mosquitoes in order to eliminate malaria in the considered
<figure><img src="image_6.png" /><figcaption>Figure 5: The evolution of $I_h(t)$ (green curve), solution of (4), and its estimate $\hat{I}_h(t)$ (red curve) delivered by the observer (25).</figcaption></figure>

region. We denote the pesticide-induced death rate for mosquitoes by $\nu_v$, in this case our model is written as follow:

$$ 
\left\{ 
\begin{aligned}
\frac{dS_h}{dt} &= \Lambda_h + \rho_h R_h + \alpha_h I_h - \beta_{vh} \frac{I_h}{H} S_h - \mu_h S_h, \\
\frac{dI_h}{dt} &= \beta_{vh} \frac{I_h}{H} S_h - \alpha_h I_h - \gamma_h I_h - \mu_h I_h - \delta_h I_h, \\
\frac{dR_h}{dt} &= \gamma_h I_h - \rho_h R_h - \mu_h R_h \\
\frac{dS_v}{dt} &= \Lambda_v - \left( \beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H} \right) S_v - \mu_v S_v - \nu_v S_v, \\
\frac{dI_v}{dt} &= \left( \beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H} \right) S_v - \mu_v I_v - \nu_v I_v.
\end{aligned}
\right.
\qquad (26) $$

Let us denote by $\tilde{\mathcal{R}}_0$ the basic reproduction number of System (26). Using the same method described in [9], we get:

$$ \tilde{\mathcal{R}}_0 = \frac{1}{\mu_v + \nu_v} \sqrt{\beta_{vh} \frac{\mu_h}{\epsilon_h} \frac{\Lambda_v}{\Lambda_h} \left( \beta_{hv} + \frac{\gamma_h}{\theta_h} \hat{\beta}_{hv} \right)}. \qquad (27) $$

We aim to keep $\tilde{\mathcal{R}}_0 < \sqrt{\frac{\mu_h}{\mu_h + \delta_h}}$ because we have proved that the DFE is GAS in the general case when $\delta_h \neq 0$ if the latest condition is satisfied. Therefore the pesticide induce death rate $\nu_v$ has to be as follow:

$$ \nu_v > \mu_v \sqrt{1 + \frac{\delta_h}{\mu_h}} \left( \mathcal{R}_0 - \sqrt{\frac{\mu_h}{\mu_h + \delta_h}} \right) \qquad (28) $$
## 4.2 Treatment of infected individuals expressed as a feedback

Now we investigate the role of the treatment in controlling the disease. The treatment rate will be considered as a control $u$. We suppose that the recovered individuals due to the treatment control go back to the susceptible class at a rate $u$. The goal is to compute $u$ as a function of the state that makes the DFE globally asymptotically stable. Taking into account the treatment, the model becomes

$$ 
\left\{
\begin{aligned}
\frac{dS_h}{dt} &= \Lambda_h + \rho_h R_h + (\alpha_h + u)I_h - \beta_{vh} \frac{I_v}{H} S_h - \mu_h S_h, \\
\frac{dI_h}{dt} &= \beta_{vh} \frac{I_v}{H} S_h - (\alpha_h + \gamma_h + \mu_h + \delta_h)I_h - u I_h \\
\frac{dR_h}{dt} &= \gamma_h I_h - \rho_h R_h - \mu_h R_h \\
\frac{dI_v}{dt} &= \left(\beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H}\right) \left(\frac{\Lambda_v}{\mu_v} - I_v\right) - \mu_v I_v.
\end{aligned}
\right.
\qquad (29) $$

Using the variables $H = S_h + I_h + R_h$, $I_h$, $R_h$, and $I_v$, System (29) can be written:

$$ \dot{x} = X(x) + uY(x), \qquad (30) $$

where $x = (H, I_h, R_h, I_v)^T$, $u \in \mathbb{R}^+$, and

$$ X = \begin{pmatrix} X_1 \\ X_2 \\ X_3 \\ X_4 \end{pmatrix} = \begin{pmatrix} -\mu_h H - \delta I_h + \Lambda_h \\ \frac{\beta_{vh} I_v (H - I_h - R)}{H} - (\alpha_h + \gamma_h + \mu_h + \delta_h) I_h \\ \gamma I_h - (\rho + \mu_h) R \\ \left(\beta_{hv} \frac{I_h}{H} + \hat{\beta}_{hv} \frac{R_h}{H}\right) \left(\frac{\Lambda_v}{\mu_v} - I_v\right) - \mu_v I_v \end{pmatrix}, \quad Y = \begin{pmatrix} 0 \\ -I_h \\ 0 \\ 0 \end{pmatrix}. $$

The set $\Omega = \{ \frac{\Lambda_h}{\mu_h + \delta_h} \le H \le \frac{\Lambda_h}{\mu_h}, 0 \le I_v \le \frac{\Lambda_v}{\mu_v} \}$ is a positively invariant and attractive set.

Hence it is sufficient to consider system (29) on the set $\Omega$. The problem under consideration is to construct a feedback law $u(x)$ in such a way that the DFE is a globally asymptotically stable equilibrium for the closed-loop system $\dot{x} = X(x) + u(x)Y(x)$. To achieve this goal we shall construct a "Control Lyapunov Function" (CLF) for System (29) or System (30).

**Theorem 3.** The following function $W$ is a CLF for system (30):

$$ W(x) = H - \frac{\Lambda_h \ln(H)}{\mu_h} + \frac{\mu_v}{\beta_{vh}} I_h + \frac{\hat{\beta}_{hv} \Lambda_v (\mu_h + \delta_h)}{\mu_v^2 \mu_h \Lambda_h} (\mu_v + 1) R_h + \frac{\mu_v + 1}{\mu_v} I_v - \frac{\Lambda_h}{\mu_h} \left(1 - \ln\left(\frac{\Lambda_h}{\mu_h}\right)\right) $$

A stabilizing feedback control is then given by ([19]) :

$$ u(x) = \begin{cases} 0 & \text{if } I_h = 0 \\ -\frac{a(H, I_h, R_h, I_v) + \sqrt{(a(H, I_h, R_h, I_v))^2 + (b(H, I_h, R_h, I_v))^4}}{b(H, I_h, R_h, I_v)\left(1 + \sqrt{1 + b(H, I_h, R_h, I_v)^2}\right)} & \text{if } I_h \neq 0 \end{cases} \qquad (31) $$

where,

$$ a(H, I_h, R_h, I_v) = \left(1 - \frac{\Lambda_h}{\mu_h H}\right) X_1 + \frac{\mu_v}{\beta_{vh}} X_2 + \frac{\hat{\beta}_{hv} \Lambda_v (\mu_h + \delta_h)}{\mu_v^2 \mu_h \Lambda_h} (\mu_v + 1) X_3 + \frac{\mu_v + 1}{\mu_v} X_4 $$

$$ b(H, I_h, R_h, I_v) = \frac{\mu_v}{\beta_{vh}} I_h $$
This feedback control makes the DFE a globally asymptotically stable equilibrium for the closed-loop system (30-31), i.e., all the trajectories will converge to the DFE.

*Proof.* The function $W$ is well defined, differentiable on $\Omega$. For all $x \in \Omega$, $W(x) \ge 0$ et $W(x) = 0$ if and only if $x = \left(\frac{\Lambda_h}{\mu_h}, 0, 0, 0\right) = \text{DFE}$. Thus $W$ is definite positive on $\Omega$.

As in [25], we use the notations: $a(x) = \langle \nabla W(x), X(x) \rangle$ and $b(x) = \langle \nabla W(x), Y(x) \rangle$. Hence,

$$ \dot{W} = a(x) + u b(x). $$

To show that $W$ is a CLF for system (30), we have to show that (cf [25]):

$$ \forall x \in \Omega : b(x) = 0 \implies a(x) < 0. \qquad (32) $$

We have on one hand, $b(H, I_h, R_h, I_v) = \langle \nabla W(x), Y(x) \rangle = -\frac{\mu_v}{\beta_{vh}} I_h$.

On the other hand:

$$ \begin{aligned} a(H, I_h, R_h, I_v) &= \langle \nabla W(x), X(x) \rangle \\ &= \left(1 - \frac{\Lambda_h}{\mu_h H}\right) (-\mu_h H - I_h \delta_h + \Lambda_h) \\ &\quad + \frac{\mu_v}{\beta_{vh}} \left( \frac{\beta_{vh} I_v (H - I_h - R_h)}{H} - \epsilon_h I_h \right) \\ &\quad + \frac{\hat{\beta}_{hv} \Lambda_v (\mu_h + \delta_h)}{\mu_v^2 \mu_h \Lambda_h} (\mu_v + 1) (\gamma_h I_h - \theta_h R_h) \\ &\quad + \frac{\mu_v + 1}{\mu_v} \left( \left( \frac{\beta_{hv} I_h}{H} + \frac{\hat{\beta}_{hv} R_h}{H} \right) \left( \frac{\Lambda_v}{\mu_v} - I_v \right) - \mu_v I_v \right) \end{aligned} $$

Since $b(x) = 0 \implies I_h = 0$, we have to show that $a(H, 0, R_h, I_v) < 0$ for all $(H, R_h, I_v)$.

$$ \begin{aligned} a(H, 0, R_h, I_v) &= \left(1 - \frac{\Lambda_h}{\mu_h H}\right) (-\mu_h H + \Lambda_h) + \mu_v \frac{I_v (H - R_h)}{H} \\ &= -\frac{\hat{\beta}_{hv} \Lambda_v (\mu_h + \delta_h)}{\mu_v^2 \mu_h \Lambda_h} (\mu_v + 1) \theta_h R_h + \frac{\mu_v + 1}{\mu_v} \left( \frac{\hat{\beta}_{hv} R_h}{H} \left( \frac{\Lambda_v}{\mu_v} - I_v \right) - \mu_v I_v \right) \\ &= -\frac{(\mu_h H - \Lambda_h)^2}{\mu_h H} - I_v \\ &\quad + \left( -\left( \frac{\mu_v}{H} + \frac{\mu_v + 1}{\mu_v} \frac{\hat{\beta}_{hv}}{H} \right) I_v - \frac{\hat{\beta}_{hv} \Lambda_v (\mu_h + \delta_h)}{\mu_v^2 \mu_h \Lambda_h} (\mu_v + 1) \theta_h + \frac{\mu_v + 1}{\mu_v} \frac{\hat{\beta}_{hv}}{H} \frac{\Lambda_v}{\mu_v} \right) R_h \end{aligned} $$

Since $\frac{\Lambda_h}{\mu_h + \delta_h} \le H \le \frac{\Lambda_h}{\mu_h}$, we can write,
$$ 
\begin{aligned}
& \left( -\frac{\hat{\beta}_{hv}\Lambda_v(\mu_h + \delta_h)}{\mu_v^2 \mu_h \Lambda_h}(\mu_v + 1)\theta_h + \frac{\frac{\mu_v + 1}{\mu_v}\hat{\beta}_{hv}\Lambda_v}{H\mu_v} \right) R_h \le \\
& \left( -\frac{\hat{\beta}_{hv}\Lambda_v(\mu_h + \delta_h)}{\mu_v^2 \mu_h \Lambda_h}(\mu_v + 1)(\rho_h + \mu_h) + \frac{\frac{\mu_v + 1}{\mu_v}\hat{\beta}_{hv}\Lambda_v(\mu_h + \delta_h)}{\mu_v \Lambda_h} \right) R_h \\
& = \frac{\hat{\beta}_{hv}\Lambda_v(\mu_h + \delta_h)}{\mu_v^2 \mu_h \Lambda_h}(\mu_v + 1)(-(\rho_h + \mu_h) + \mu_h)R_h.
\end{aligned}
$$

We then obtain

$$ \left( -\frac{\hat{\beta}_{hv}\Lambda_v(\mu_h + \delta_h)}{\mu_v^2 \mu_h \Lambda_h}(\mu_v + 1)\theta_h + \frac{\frac{\mu_v + 1}{\mu_v}\hat{\beta}_{hv}\Lambda_v}{H\mu_v} \right) R_h \le -\frac{\hat{\beta}_{hv}\Lambda_v(\mu_h + \delta_h)}{\mu_v^2 \mu_h \Lambda_h}(\mu_v + 1)\rho_h R_h. $$

Thus,

$$ a(H, 0, R_h, I_v) \le -\frac{(\mu_h H - \Lambda_h)^2}{\mu_h H} - \frac{\hat{\beta}_{hv}\Lambda_v(\mu_h + \delta_h)}{\mu_v^2 \mu_h \Lambda_h}(\mu_v+1) \rho_h R_h - I_v - \left(\mu_v + \frac{\mu_v + 1}{\mu_v}\hat{\beta}_{hv}\right) \frac{R_h I_v}{H} < 0. $$

This proves that the function W is a CLF for System (30).

As in [19], we construct the stabilizing control as follows ($x = (H, I_h, R_h, I_v)$):

$$ u(x) = \begin{cases} 0 & \text{for } I_h = 0 \\ -\frac{a(x) + \sqrt{(a(x))^2 + (b(x))^4}}{b(x)(1 + \sqrt{1 + b(x)^2})} & \text{for } I_h \neq 0 \end{cases} \qquad (33) $$

<figure><img src="image_7.png" /></figure>

### 4.2.1 Simulations

We shall illustrate the effect of the treatment used as a feedback control and computed according to formula (31). We use the following parameters $\alpha_h = 0.0005$, $\delta_h = 5 \times 10^{-5}$, $\gamma_h = 0.00035$, $\rho_h = 0.00083$, $\Lambda_h = 1$, $\Lambda_v = 400$, $\mu_h = 5 \times 10^{-5}$, $\mu_v = 0.02\beta_{hv} = 0.027$, $\beta_{vh} = 0.029$, $\hat{\beta}_{hv} = 10^{-4}$. With these values, the basic reproduction number is equal to $R_0 = 1.91$. Figures 6-7 display the dynamics of infected populations ($I_h$ and $I_v$), given by Model (29) with no treatment ($u = 0$). These figures show that malaria remains persistent in the two populations. Figures 8-9 display the evolution of infected populations when a treatment is applied and its rate is given by the feedback control (31). It can be observed that the application of the feedback (31) make the disease dies out in both populations.

Figure 10 compare the evolution of human infected population ($I_h$) when no treatment is applied and when it is applied with a rate given by (31).
<figure><img src="image_8.png" /><figcaption>Figure 6: $I_h(t)$ without control</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: $I_v(t)$ without control</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: $I_h(t)$ with feedback control</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: $I_v(t)$ with feedback control</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10: Evolution of $I_h(t)$ solution of (29) with no control (red curve) and when feedback control (31) is applied (green curve)</figcaption></figure>

# 5 Conclusion

In this paper, we investigated the behavior and the control of a mathematical model for malaria transmission, based on the model proposed by Arino et al. [3].

When $\mathcal{R}_0 \le 1$, we proved the global stability of the disease free equilibrium (DFE) when there is no disease induced mortality. When the disease induced death rate $\delta_h$ is positive global asymptotic stability of the DFE cannot be, in general, expected with the only condition $\mathcal{R}_0 \le 1$ since it has been proved in the literature that the system may have endemic equilibria even if $\mathcal{R}_0 \le 1$. In this case we proved that global asymptotic stability of the DFE occurs if $\mathcal{R}_0 \le \sqrt{\frac{\mu_h}{\mu_h + \delta_h}}$.

When $\mathcal{R}_0 > 1$, we showed that the disease is uniformly persistent and proved the existence of a unique endemic equilibrium point whose local and global asymptotic stability are investigated.

We also investigated the state estimation problem for this model. We gave some tools that allow to estimate the various state variables using the only available data, namely the host incidence, i.e. the number of new infected humans per unit time, provided by Public Health Services. Our method consists of using elements of system theory in designing an auxiliary dynamical system, called observer, whose solutions converge exponentially to those of the original model.

Current strategies to control mosquito-transmitted infections use pesticides, we discussed in this work a possible way to control the impact of the disease using a pesticide death rate. We showed that the pesticide death rate have to be larger than a specific constant that depend basically on the basic reproduction number of the system $\mathcal{R}_0$ in order to control the situation and to eradicate the disease. We also addressed the problem of controlling the disease evolution using treatment of infected individuals. We gave a formula to compute the treatment rate as a feedback control. This feedback is based on the construction of a Control Lyapunov Function.
## Authors' contribution

All authors have participated in (a) conception, design, and analysis and interpretation of the research; (b) drafting the article or revising it critically for important intellectual content; and (c) approval of the final version.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

A. Iggidr and S. Yacheur acknowledge the financial support of the program STIC AmSud (project 20-STIC-05 - NEMBICA). A part of this work was completed when A. Iggidr was visiting ENSA (Universiy Cadi Ayyad, Marrakech).

## Data availability

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

## References

[1] R. M. Anderson and R. M. May. *Infectious diseases of humans: dynamics and control*. Oxford University Press, 1991.

[2] R. Anguelov, Y. Dumont, J. Lubuma, and E. Mureithi. Stability analysis and dynamics preserving nonstandard finite difference schemes for a malaria model. *Mathematical Population Studies*, 20(2):101–122, apr 2013. doi:10.1080/08898480.2013.777240.

[3] J. Arino, A. Ducrot, and P. Zongo. A metapopulation model for malaria with transmission-blocking partial immunity in hosts. *J Math Biol*, 64(3):423–48, Feb 2012. doi:10.1007/s00285-011-0418-4.

[4] J. Aron and R. May. *The population dynamics of malaria*, pages 139–179. Springer, Boston, MA, 1982. doi:10.1007/978-1-4899-2901-3_5.

[5] R. C. I. Barradas. Vector preference annihilates backward bifurcation and reduces endemicity. *Bulletin of Mathematical Biology*, 81:4447–4469, 2018. doi:10.1007/s11538-018-00561-1.

[6] G. Besançon, editor. *Nonlinear observers and applications. Papers based on the presentations at the 28th Grenoble international summer school on control, Grenoble, France, September 2007*. Lecture Notes in Control and Information Sciences 363. Berlin: Springer. xiii, 224 p. , 2007. doi:10.1007/978-3-540-73503-8.
[7] N. Chitnis, J. M. Cushing, and J. M. Hyman. Bifurcation analysis of a mathematical model for malaria transmission. *SIAM Journal on Applied Mathematics*, 67(1):24–45, 2006. doi: 10.1137/050638941.

[8] C. Chiyakaa, W. Gariraa, and S. Dubeb. Transmission model of endemic human malaria in a partially immune population. *Mathematical and Computer Modelling*, 46:806–822, September 2007. doi:10.1016/j.mcm.2006.12.010.

[9] P. V. den Driessche and J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Mathematical biosciences*, 180(1-2):29–48, 2002. doi:10.1016/S0025-5564(02)00108-6.

[10] A. Ducrot, S. B. Sirima, B. Somé, and P. Zongo. A mathematical model for malaria involving differential susceptibility, exposedness and infectivity of human host. *Journal of Biological Dynamics*, 3(6):574–598, 2009. doi:10.1080/17513750902829393.

[11] X. Feng, S. Ruan, Z. Teng, and K. Wang. Stability and backward bifurcation in a malaria transmission model with applications to the control of malaria in china. *Mathematical Biosciences*, 266:52–64, 2015. doi:10.1016/j.mbs.2015.05.005.

[12] J.-P. Gauthier, H. Hammouri, and S. Othman. A simple observer for nonlinear systems applications to bioreactors. *IEEE Trans. Autom. Control*, 37(6):875–880, 1992. doi:10.1109/9.256352.

[13] H. W. Hethcote and H. R. Thieme. Stability of the endemic equilibrium in epidemic models with subpopulations. *Math. Biosci.*, 75:205–227, 1985. doi:10.1016/0025-5564(85)90038-0.

[14] J. Hofbauer and J. W.-H. So. Uniform persistence and repellors for maps. *Proc. Am. Math. Soc.*, 107(4):1137–1142, 1989. doi:10.2307/2047679.

[15] J. C. Kamgang and G. Sallet. Computation of threshold conditions for epidemiological models and global stability of the disease-free equilibrium (DFE). *Math. Biosci.*, 213(1):1–12, 2008. doi:10.1016/j.mbs.2008.02.005.

[16] A. A. Kareem and A. R. Kimbir. Modeling the epidemiology of malaria and control with estimate of the basic repro-duction number. *Pure and Applied Mathematics Journal*, 2(1):42, 2013. doi:10.11648/j.pamj.20130201.17.

[17] M. Krasnosel'skij. *Positive solutions of operator equations. Translated from the Russian by Richard E. Flaherty. Edited by Leo F. Boron*. Groningen: The Netherlands: P.Noordhoff Ltd. 381 p. , 1964.

[18] M. Y. Li and J. S. Muldowney. A geometric approach to global-stability problems. *SIAM J. Math. Anal.*, 27(4):1070–1083, 1996. doi:10.1137/S0036141094266449.

[19] Y. Lin and E. D. Sontag. A universal formula for stabilization with bounded controls. *Systems & Control Letters*, 16(6):393 – 397, 1991. doi:10.1016/0167-6911(91)90111-Q.

[20] G. Macdonald. *The epidemiology and control of malaria*. Oxford University Press, 1957.

[21] C. McCluskey and P. Driessche. Global analysis of two tuberculosis models. *Journal of Dynamics and Differential Equations*, 16(1):139–166, 2004. doi:10.1023/B:JODY.0000041283.66784.3e.
[22] J. L. Mojeeb AL-Rahman EL-Nor Osman. Analysis of a vector-bias malaria transmission model with application to mexico, sudan and democratic republic of the congo. *Journal of Theoretical Biology*, 464:72–84, (2018). doi:10.1016/j.jtbi.2018.12.033.

[23] S. B. R. Senthamarai. Analytical solutions of sirssi malaria disease model using hpm. *Journal of Chemical and Pharmaceutical Research*, 8(7):651–666, 2016. URL: http://www.jocpr.com/articles/analytical-solutions-of-sirssi-malaria-disease-model-using-hpm.pdf.

[24] R. Ross. *The prevention of malaria*. John Murray, London, 1911.

[25] E. D. Sontag. A “universal” construction of Artstein’s theorem on nonlinear stabilization. *Systems Control Lett.*, 13(2):117–123, 1989. doi:10.1016/0167-6911(89)90028-5.

[26] M. Vidyasagar. Decomposition techniques for large-scale systems with nonadditive interactions: Stability and stabilizability. *IEEE Trans. Autom. Control*, 25(4):773–779, 1980. doi:10.1109/TAC.1980.1102422.

[27] WHO. World health organization. http://www.who.int/malaria/en/.

[28] WHO. World health organization travellers. https://www.who.int/malaria/travellers/en/.

[29] H. Yang. Malaria transmission model for different levels of acquired immunity and temperature-dependent parameters (vector). *Revista de saúde pública*, 34:223–31, 07 2000. doi:10.1590/S0034-89102000000300003.