1. Introduction

In many infectious diseases, such as HIV, schistosomiasis, tuberculosis, the infectiousness of an infected individual can be very different at various stages of infection. Hence, the age of infection may be an important factor to consider in modeling transmission dynamics of infectious diseases. In the epidemic model of Kermack and Mckendrick [9], infectivity is allowed to depend on the age of infection. Because the age-structured epidemic model is described by first order PDEs, it is more difficult to theoretically analyze the dynamical behavior of the PDE models, particularly the global stability. Several recent studies [10,11,18] have focused on age structured models, and the results show that age of infection may play an important role in the transmission dynamics of infectious diseases.

In our pervious work [4], we formulated an infection-age structured epidemic model to describe the transmission dynamics of vector-borne diseases. The model includes both incubation age of the exposed hosts and infection age of the infectious hosts, both of which describe the different removal rates in the latent period and the variable infectiousness in the infectious period, respectively. The results in [4] show that the basic reproduction number determines transmission dynamics of vector-borne diseases: the disease-free equilibrium is globally asymptotically stable if the basic reproduction number is less than 1, and the endemic equilibrium is globally asymptotically stable if the basic reproduction number is greater than 1. However, the vector-borne epidemic model formulated in [4] only incorporates a single strain. In reality, many diseases are caused by more than one antigenically different strains of the causative agent [15]. For instance, the dengue virus has 4 different serotypes [6], and bacterial pneumonia is caused by more than ninety different serotypes of Streptoccus pneumoniae. Therefore, it is necessary to study infection-age structured epidemic models with multiple strains.

In this paper, we will extend the model with a single strain to the model with multiple strains, and obtain the following infection-age-structured vector-borne epidemic model with multiple strains:

$ \left\{ $$\begin{array}{l} \displaystyle \dfrac{dS_v}{dt}= \Lambda_v- \sum\limits_{j=1}^nS_v \int_0^\infty\beta_v^j(a) I_h^j(a, t)da-\mu_v S_v,\\ \displaystyle \dfrac{dI_v^j}{dt}=S_v \int_0^\infty\beta_v^j(a) I_h^j(a, t)da- (\mu_v+\alpha_v^j)I_v^j,\\ \displaystyle \dfrac{dR_v}{dt}=\sum\limits_{j=1}^n\alpha_v^jI_v^j- \mu_vR_v,\\[2ex] \displaystyle \dfrac{dS_h}{dt}=\Lambda_h-\sum\limits_{j=1}^n\beta_h^jS_hI_v^j - \mu_hS_h,\\ \displaystyle \frac{\partial E_h^j(\tau, t)}{\partial\tau}+\frac{\partial E_h^j(\tau, t)}{\partial t}=-(\mu_h+m_h^j(\tau))E_h^j(\tau, t),\\[2ex] \displaystyle E_h^j(0, t)=\beta_h^j S_hI_v^j ,\\[2ex] \displaystyle \frac{\partial I_h^j(a, t)}{\partial a}+\frac{\partial I_h^j(a, t)}{\partial t}=-(\mu_h+\alpha_h^j(a)+r_h^j(a))I_h^j(a, t),\\[2ex] \displaystyle I_h^j(0, t)=\int^\infty_0 m_h^j(\tau)E_h^j(\tau, t)d\tau,\\[2ex] \displaystyle \dfrac{dR_h}{dt}= \sum\limits_{j=1}^n\int^\infty_0r_h^j(a) I_h^j(a, t)da- \mu_h R_h. \end{array}$$ \right. $

(1)

In the model (1), $S_h(t)$, $E^j_h(\tau, t)$, $I^j_h(a, t)$, $R_h(t)$ represent the number/density of the susceptible hosts, infected hosts with strain $j$ but not infectious, infected hosts with strain $j$ and infectious, and recovered hosts at time $t$, respectively. $S_v(t)$, $I^j_v(t)$ and $R_v(t)$ denote the number of the susceptible vectors, infected vectors with strain $j$ and infectious, and recovered vectors at time $t$, respectively. $\Lambda_v,\Lambda_h$ are the birth /recruitment rates of the vectors and hosts, respectively; $\mu_v, \mu_h$ are the natural death rates of the vectors and hosts, respectively. The parameter $m^j_h(\tau)$ denotes the removal rate of the infected hosts with strain $j$ of incubation age $\tau$ from the latent period; $\alpha^j_h(a)$ is the additional disease induced death rate due to the strain $j$ at age of infection $a$; $\alpha^j_v$ denotes the recovery rate of the infected vectors with strain $j$; $ r^j_h(a)$ denotes the recovery rate of the infected hosts of infection age $a$ with strain $j$; $\beta^j_v(a)$ is the transmission coefficient of the infected host individuals with strain $j$ at age of infection $a$, and $\beta^j_h$ is the transmission coefficient of strain $j$ from infected vectors to healthy host individuals.

The dynamics of the epidemic model involving multiple strains has fascinated researchers for a long time (see [3,5,6,7,17] and the references therein), and one of the important results is the competitive exclusion principle. In epidemiology, the competitive exclusion principle states that if multiple strains circulate in the population, only the strain with the largest reproduction number persists and the strains with suboptimal reproduction numbers are eliminated [13]. Using a multiple-strain ODE model Bremermann and Thieme [2] first proved that the principle of competitive exclusion is valid under the assumption that infection with one strain precludes additional infections with other strains. In 2013, Maracheva and Li [13] extended the competitive exclusion principle to a multi-stain age-since-infection structured model of SIR/SI-type. The goal of this paper is to extend this principle to model (1).

As we all know, the proof of competitive exclusion principle is based on the global stability of the single-strain equilibrium. The stability analysis of nonlinear dynamical systems has always been an important topic theoretically and practically since global stability is one of the most important issues related to their dynamic behaviors. Due to the lack of generically applicable tools proving the global stability is very challenging, especially for the continuous age-structured models which are described by first order PDEs. Although there are various approaches for some general nonlinear systems, the method of Lyapunov functions is the most common tool used to prove the global stability. In this paper, we will apply a class of Lyaponuv functions to study the global dynamics of system (1) and draw on the results to derive the competitive exclusion principle for infinite dimensional systems.

This paper is organized as follows. In the next section we derive an explicit formula for the basic reproduction number $\mathcal R^j_0$ of strain $j$ for $j=1,\cdots,n$, and then we will show that strain $j$ will die out if its basic reproduction number is less than one. In section 3, we will define the disease reproduction number $\mathcal R_0$, and then prove that the disease-free equilibrium (DFE) of the system is globally asymptotically stable if $\mathcal R_0< 1$. In Section 4, we will investigate the existence of single-strain equilibria and their local stabilities. In section 5 we will devote to prove the principle of competitive exclusion. Without loss of generality, we assume that strain one has the maximal reproduction number and $\mathcal R^1_0>1$. Under the assumption, we will show that strain one is uniformly strong persistent while the remaining strains become extinct. In Section 6, we use a class of Lyapunov functions to derive the global stability of the strain one equilibrium under the condition that $\mathcal R^i_0/\mathcal R^1_0<b_i/b_1<1,i\neq 1$, where $b_j$ denotes the probability of a given susceptible vector being transmitted by an infected host with strain $j$, which implies that complete competitive exclusion holds for the system. Finally, a brief discussion is given in Section 7.

---

2. The reproduction numbers and threshold dynamics

In this section, we mainly derive the reproduction numbers for each strain, and show that the stain will die out if its basic reproduction number is less than one.

Since the equations for the recovered individuals and the recovered vectors are decoupled from the system, it follows that the dynamical behavior of system (1) is equivalent to the dynamical behavior of the following system:

$ \left\{ $$\begin{array}{l} \displaystyle \dfrac{dS_v}{dt}= \Lambda_v- \sum\limits_{j=1}^nS_v \int_0^\infty\beta_v^j(a) I_h^j(a, t)da-\mu_v S_v,\\[2ex] \displaystyle \dfrac{dI_v^j}{dt}=S_v \int_0^\infty\beta_v^j(a) I_h^j(a, t)da- (\mu_v+\alpha_v^j)I_v^j,\\[2ex] \displaystyle \dfrac{dS_h}{dt}=\Lambda_h-\sum\limits_{j=1}^n\beta_h^jS_hI_v^j - \mu_hS_h,\\[2ex] \displaystyle \frac{\partial E_h^j(\tau, t)}{\partial\tau}+\frac{\partial E_h^j(\tau, t)}{\partial t}=-(\mu_h+m_h^j(\tau))E_h^j(\tau, t),\\[2ex] \displaystyle E_h^j(0, t)=\beta_h^j S_hI_v^j ,\\[2ex] \displaystyle \frac{\partial I_h^j(a, t)}{\partial a}+\frac{\partial I_h^j(a, t)}{\partial t}=-(\mu_h+\alpha_h^j(a)+r_h^j(a))I_h^j(a, t),\\[2ex] \displaystyle I_h^j(0, t)=\int^\infty_0 m_h^j(\tau)E_h^j(\tau, t)d\tau. \end{array}$$ \right. $

(2)

Model (2) is equipped with the following initial conditions:

$S_v(0)=S_{v_0}, I_v^j(0)=I_{v_0}^j, S_h(0)=S_{h_0} , E_h^j(\tau, 0)=\varphi_j(\tau), I_h^j(a, 0)=\psi_j(a).$

All parameters are nonnegative, $\Lambda_v > 0, ~\Lambda_h > 0$, and $\mu_v > 0, ~\mu_h > 0$. We make the following assumptions on the parameter-functions.

Assumption 2.1.

1. The function $\beta_v^j(a)$ is bounded and uniformly continuous for every $j$. When $\beta_v^j(a)$ is of compact support, the support has non-zero Lebesgue measure;

2. The functions $m_h^j(\tau),\ \alpha_h^j(a),\ r_h^j(a)$ belong to $L^\infty(0, \infty)$;

3. The functions $ \varphi_j(\tau),\ \psi_j(a)$ are integrable.

Let us define

$ $$\begin{align*} X=\mathbb{R}\times\prod^{n}_{j=1}\mathbb{R}\times\mathbb{R}\times\prod^{n}_{j=1}(L^1(0, \infty)\times L^1(0, \infty)). \end{align*}$$ $

It is easily verified that solutions of (2) with nonnegative initial conditions belong to the positive cone for $t \geq 0$. Adding the first and all equations for $I_v^j$ yields that

$ $$\begin{align*} \displaystyle \frac{d}{dt}\bigg(S_v(t)+\sum^n_{j=1}I_v^j( t) \bigg) \leq \displaystyle \Lambda_v-\mu_v\bigg(S_v(t)+\sum^n_{j=1}I_v^j( t) \bigg). \end{align*}$$ $

Hence,

$ $$\begin{align*} \limsup\limits_{t\rightarrow +\infty}\bigg(S_v(t)+\sum^n_{j=1}I_v^j( t)\bigg) \leq \frac{\Lambda_v}{\mu_v}. \end{align*}$$ $

Similarly, adding the equation for $S_h$ and all equations for $E_h^j,I_h^j$, we have

$ $$\begin{align*} \displaystyle \frac{d}{dt}\bigg(\displaystyle S_h(t) & +\displaystyle \sum^n_{j=1}\int^\infty_0E_h^j(\tau, t)d\tau+\sum^n_{j=1}\int^\infty_0I_h^j(a, t)da\bigg) \displaystyle \end{align*}$$ $

$ $$\begin{align*} & \displaystyle \leq\Lambda_h-\mu_h\bigg(S_h(t)+\sum^n_{j=1}\int^\infty_0E_h^j(\tau, t)d\tau+\sum^n_{j=1}\int^\infty_0I_h^j(a, t)da\bigg), \end{align*}$$ $

and it then follows that

$ $$\begin{align*} \displaystyle \limsup\limits_{t\rightarrow +\infty} \bigg(S_h(t)+\sum^n_{j=1}\int^\infty_0E_h^j(\tau, t)d\tau+\sum^n_{j=1}\int^\infty_0I_h^j(a, t)da\bigg) \leq\frac{\Lambda_h}{\mu_h}. \end{align*}$$ $

Therefore, the following set is positively invariant for system (2)

$ $$\begin{split} \displaystyle \Omega =\bigg\{(S_v, ~I_v^1, & ~\cdots,~I_v^n,~ S_h, ~E_h^1, ~I_h^1,~\cdots,~E_h^n, ~I_h^n )\in X_+\bigg|\\ & \displaystyle\bigg(S_v(t)+\sum^n_{j=1}I_v^j( t) \bigg) \leq \frac{\Lambda_v}{\mu_v},\\ & \displaystyle \bigg(S_h(t)+\sum^n_{j=1}\int^\infty_0E_h^j(\tau, t)d\tau+\sum^n_{j=1}\int^\infty_0I_h^j(a, t)da\bigg) \leq\frac{\Lambda_h}{\mu_h}\bigg\}. \end{split}$$ $

(3)

In what follows, we only consider the solutions of the system (2) with initial conditions which lie in the region $\Omega$. As we all know, the reproduction number is one of most important concepts in epidemiological model. Next, we will express the basic reproduction numbers for each strain. To simplify expression, let us introduce two notations.

Definition 2.1. The exit rate of exposed host individuals with strain $j$ from the incubation compartment is given by $\mu_h+m_h^j(\tau)$, the probability of still being latent after $\tau$ time units, denoted by $\pi^j_1(\tau)$, is given by

$ $$\begin{array}{l} \displaystyle \pi^j_1(\tau)= e^{-\mu_h\tau}e^{-\int^\tau_0m_h^j(\sigma))d\sigma}. \end{array}$$ $

(4)

Definition 2.2. The exit rate of infected individuals with strain $j$ from the infective compartment is given by $\mu_h+\alpha_h^j(a)+r_h^j(a)$, and it then follows that the probability of still being infectious after $a$ time units, denoted by $\pi^j_2(a)$, is given by

$ \label{eq:4} $$\begin{array}{l} \displaystyle\pi^j_2(a)= e^{-\mu_ha}e^{-\int^a_0 (\alpha_h^j(\sigma)+r_h^j(\sigma))d\sigma}. \end{array}$$ $

(5)

Then we can give the expression for the basic reproduction number of strain $j$ which can be expressed as

$ $$\begin{split} \displaystyle \mathcal R^j_0=\displaystyle \frac{\beta_h^j\Lambda_v\Lambda_h}{\mu_v\mu_h(\mu_v+\alpha_v^j)}\int^\infty_0m_h^j(\tau)\pi^j_1(\tau)d\tau\int^\infty_0\beta_v^j(a)\pi^j_2(a)da. \end{split}$$ $

(6)

The reproduction number of strain $j$ gives the number of secondary infections produced in an entirely susceptible population by a typical infected individual with strain $j$ during its entire infectious period. $\mathcal{R}^j_0$ gives the strength of strain $j$ to invade into the system when rare and alone. The reproduction number of strain $j$ consists of two terms:

$ $$\begin{align*}\displaystyle \mathcal R^j_h=\frac{\Lambda_v}{\mu_v}\int^\infty_0\beta_v^j(a)\pi^j_2(a)da, \mathcal R^j_v=\frac{\beta_h^j\Lambda_h}{\mu_h(\mu_v+\alpha_v^j)}\int^\infty_0m_h^j(\tau)\pi^j_1(\tau)d\tau. \end{align*}$$ $

The first term $\mathcal R^j_h$ represents the reproduction number of human-to-vector transmission of strain $j$, and the second term $\mathcal R^j_v$ is the reproduction number of vector-to-human transmission of strain $j$.

Now we are able to state the results on threshold dynamics of strain $j$:

Theorem 2.3. If $ \mathcal R^j_0<1$, strain $j$ will die out.

Proof. Let

$ $$\begin{align*} B_E^j(t)=E_h^j(0,t), B_I^j(t)=I_h^j(0, t). \end{align*}$$ $

Integrating along the characteristic lines of system (2) yields

$ $$\begin{array}{ll} \displaystyle E_h^j(\tau, t)= \left\{\begin{array} {ll} \displaystyle B_E^j(t-\tau)\pi^j_1(\tau), \mbox{ } & t>\tau,\\[1ex] \displaystyle \varphi_j(\tau-t)\frac{\pi^j_1(\tau)}{\pi^j_1(\tau-t)}, \mbox{ } & t<\tau, \end{array}$$ \right.\\[3ex] \displaystyle I_h^j (a, t)= \left\{\begin{array} {ll} \displaystyle B_I^j(t-a)\pi^j_2(a), \mbox{ } & t>a,\\[1ex] \displaystyle \psi_j(a-t)\frac{\pi^j_2(a)}{\pi^j_2(a-t)}, \mbox{ } & t<a. \end{array}\right. \end{array} $

(7)

From the first and the third equations of system (2), we obtain

$ \displaystyle \limsup\limits_{t\rightarrow +\infty} S_v(t)\leq\frac{\Lambda_v}{\mu_v}, \limsup\limits_{t\rightarrow +\infty} S_h(t)\leq\frac{\Lambda_h}{\mu_h}. $

(8)

Thus, from system (2) and inequalities (8), we have

$ \left\{ $$\begin{array}{ll} \displaystyle \frac{dI_v^j(t)}{dt}\leq \frac{\Lambda_v}{\mu_v}\int^\infty_0\beta_v^j(a)I_h^j(a,t)da-(\mu_v+\alpha_v^j)I_v^j,\\[2ex] \displaystyle E_h^j(\tau,t)=E_h^j(0,t-\tau)\pi^j_1(\tau), t>\tau,\\[2ex] \displaystyle I_h^j(a,t)=I_h^j(0,t-a)\pi^j_2(a), t>a. \end{array}$$ \right. $

(9)

From the first inequality of (9), we obtain that

$ $$\begin{array}{ll} \displaystyle I_v^j(t) & \leq\displaystyle I_v^j(0)e^{-(\mu_v+\alpha_v^j)t}+ \frac{\Lambda_v}{\mu_v}\int^t_0e^{-(\mu_v+\alpha_v^j)(t-s)}\int^\infty_0\beta_v^j(a)I_h^j(a,s)dads\\[2ex] & \leq\displaystyle I_v^j(0)e^{-(\mu_v+\alpha_v^j)t}+\frac{\Lambda_v}{\mu_v}\int^t_0e^{-(\mu_v+\alpha_v^j)(t-s)} \bigg (\int^s_0\beta_v^j(a)I_h^j(0,s-a)\pi^j_2(a)da\\[2ex] & + \displaystyle\int^t_s\beta_v^j(a)\psi_j(a-s)\frac{\pi^j_2(a)}{\pi^j_2(a-s)}da+\int^\infty_t\beta_v^j(a)I_h^j(a,s)da\bigg )ds. \end{array}$$ $

(10)

Notice that

$ $$\begin{array}{ll} & \displaystyle \limsup\limits_{t\rightarrow +\infty} \int^t_0e^{-(\mu_v+\alpha_v^j)(t-s)} \int^s_0\beta_v^j(a)I_h^j(0,s-a)\pi^j_2(a)dads\\[2ex] \leq & \displaystyle \bigg(\limsup\limits_{t\rightarrow +\infty}\int^t_0e^{-(\mu_v+\alpha_v^j)(t-s)}ds\bigg) \int^\infty_0\beta_v^j(a)\pi^j_2(a)da\bigg(\limsup\limits_{t\rightarrow +\infty}I_h^j(0,t)\bigg)\\[2ex] = & \displaystyle \frac{1}{\mu_v+\alpha_v^j} \int^\infty_0\beta_v^j(a)\pi^j_2(a)da\bigg(\limsup\limits_{t\rightarrow +\infty}I_h^j(0,t)\bigg),\\[2ex] \end{array}$$ $

(11)

$ $$\begin{array}{ll} & \displaystyle \limsup\limits_{t\rightarrow +\infty} \int^t_0e^{-(\mu_v+\alpha_v^j)(t-s)} \int^t_s\beta_v^j(a)\psi_j(a-s)\frac{\pi^j_2(a)}{\pi^j_2(a-s)}dads\\[2ex] \leq & \displaystyle \bar{\beta}\limsup\limits_{t\rightarrow +\infty} \int^t_0e^{-(\mu_v+\alpha_v^j)(t-s)} \int^t_s\psi_j(a-s)e^{-\int^a_{a-s}(\mu_h+\alpha_h^j(\sigma)+r_h^j(\sigma))d\sigma}dads\\[2ex] \leq & \displaystyle \bar{\beta}\limsup\limits_{t\rightarrow +\infty} \int^t_0e^{-(\mu_v+\alpha_v^j)(t-s)} \int^t_s\psi_j(a-s)e^{-\mu_hs}dads \end{array}$$ $

$ $$\begin{array}{ll}= & \displaystyle \bar{\beta}\limsup\limits_{t\rightarrow +\infty} \bigg(e^{-(\mu_v+\alpha_v^j)t}\int^t_0e^{(\mu_v+\alpha_v^j-\mu_h)s} \int^{t-s}_0\psi_j(a)dads\bigg)\\[2ex] = & \displaystyle \bar{\beta}\int^\infty_0\psi_j(a)da\limsup\limits_{t\rightarrow +\infty} \bigg(e^{-(\mu_v+\alpha_v^j)t}\frac{e^{(\mu_v+\alpha_v^j-\mu_h)t}-1}{\mu_v+\alpha_v^j-\mu_h}\bigg)\\[2ex] = & 0, \end{array}$$ $

(12)

and

$ $$\begin{array}{ll} & \displaystyle \limsup\limits_{t\rightarrow +\infty} \int^t_0e^{-(\mu_v+\alpha_v^j)(t-s)} \int^\infty_t\beta_v^j(a)I_h^j(a,s)dads=0. \end{array}$$ $

(13)

It then follows from (11), (12) and (13) that

$ $$\begin{array}{ll} & \displaystyle \limsup\limits_{t\rightarrow +\infty} I_v^j(t)\\[3ex] \leq & \displaystyle \frac{\Lambda_v}{\mu_v(\mu_v+\alpha_v^j)} \int^\infty_0\beta_v^j(a)\pi^j_2(a)da\bigg(\limsup\limits_{t\rightarrow +\infty}I_h^j(0,t)\bigg)\\[3ex] \leq & \displaystyle \frac{\Lambda_v}{\mu_v(\mu_v+\alpha_v^j)} \int^\infty_0\beta_v^j(a)\pi^j_2(a)da\bigg(\limsup\limits_{t\rightarrow +\infty}\int^\infty_0m_h^j(\tau)E_h^j(\tau,t)d\tau\bigg)\\[3ex] \leq & \displaystyle \frac{\Lambda_v}{\mu_v(\mu_v+\alpha_v^j)} \int^\infty_0\beta_v^j(a)\pi^j_2(a)da \limsup\limits_{t\rightarrow +\infty}\bigg(\int^t_0m_h^j(\tau)E_h^j(\tau,t)d\tau\\[3ex] & \displaystyle+\int^\infty_tm_h^j(\tau)E_h^j(\tau,t)d\tau\bigg)\\[3ex] = & \displaystyle \frac{\Lambda_v}{\mu_v(\mu_v+\alpha_v^j)} \int^\infty_0\beta_v^j(a)\pi^j_2(a)da \bigg(\limsup\limits_{t\rightarrow +\infty}\int^t_0m_h^j(\tau)E_h^j(0,t-\tau)\pi^j_1(\tau)d\tau\bigg)\\[3ex] \leq & \displaystyle \frac{\Lambda_v}{\mu_v(\mu_v+\alpha_v^j)} \int^\infty_0\beta_v^j(a)\pi^j_2(a)da\int^\infty_0m_h^j(\tau)\pi^j_1(\tau)d\tau \bigg(\limsup\limits_{t\rightarrow +\infty}E_h^j(0,t)\bigg)\\[3ex] \leq & \displaystyle \beta_h^j\frac{\Lambda_v\Lambda_h}{\mu_v\mu_h(\mu_v+\alpha_v^j)} \int^\infty_0\beta_v^j(a)\pi^j_2(a)da\int^\infty_0m_h^j(\tau)\pi^j_1(\tau)d\tau \limsup\limits_{t\rightarrow +\infty} I_v^j(t)\\[3ex] \leq & \displaystyle \mathcal R^j_0\limsup\limits_{t\rightarrow +\infty} I_v^j(t). \end{array}$$ $

(14)

Since $\mathcal R^j_0<1$ and $I_v^j(t),~j=1,\cdots,n$, are all bounded, the above expression implies that

$ \displaystyle \limsup\limits_{t\rightarrow +\infty} I_v^j(t)=0, j=1,\cdots,n. $

(15)

Hence, we have

$ \displaystyle \limsup\limits_{t\rightarrow +\infty}E_h^j(0,t)=0, \limsup\limits_{t\rightarrow +\infty}E_h^j(\tau,t)=\limsup\limits_{t\rightarrow +\infty} E_h^j(0,t-\tau)\pi^j_1(\tau)=0. $

(16)

By using the same argument, we have

$ \displaystyle \limsup\limits_{t\rightarrow +\infty}I_h^j(0,t)=0, \limsup\limits_{t\rightarrow +\infty}I_h^j(a,t)=0. $

(17)

Therefore, $(I_v^j(t),E_h^j(\tau,t),I_h^j(a,t))\rightarrow 0$ as $t\rightarrow \infty$. This means that strain $j$ will die out. The proof of Theorem 2.3 is completed.

---

3. Global stability of the disease-free equilibrium

In this section, we mainly define the disease reproduction number and show that the disease free equilibrium is globally asymptotically stable if the disease reproduction number $\mathcal R_0$ is less than one, where

$ $$\begin{align*} \mathcal R_0=\max \{\mathcal R^1_0, \cdots, \mathcal R^n_0 \}. \end{align*}$$ $

System (2) always has a unique disease-free equilibrium $\mathcal{E}_0$, which is given by

${{\cal E}_0} = \left( {S_{{v_0}}^*,{\bf{0}},S_{{h_0}}^*,{\bf{0}},{\bf{0}}} \right),$

where

$ $$\begin{align*} \displaystyle S_{v_0}^*=\frac{\Lambda_v}{\mu_v}, S_{h_0}^*=\frac{\Lambda_h}{\mu_h}, \end{align*}$$ $

and $\bf{0}=(0,\cdots,0)$ is an $n$-dimensional zero vector.

Now let us establish the local stability of the disease-free equilibrium. Let

$ $$\begin{align*} S_v(t) = & \displaystyle S_{v_0}^* +x_v(t), I_v^j{( t)} = y_v^j{( t)}, S_h(t) =S_{h_0}^* +x_h(t), \\[2ex] & \displaystyle E_h^j(\tau, t)=z_h^j(\tau, t), I_h^j(a, t)=y_h^j(a, t). \end{align*}$$ $

Then the linearized system of system (2) at the disease-free equilibrium $\mathcal{E}_0$ can be expressed as

$\left\{ $$\begin{array}{ll} \displaystyle \frac{dx_v(t)}{dt}=-\sum^n_{j=1}S_{v_0}^*\int^\infty_0\beta_v^j(a)y_h^j(a, t)da-\mu_v x_v(t),\\[4ex] \displaystyle \frac{d y_v^j( t)}{dt}=S_{v_0}^*\int^\infty_0\beta_v^j(a)y_h^j(a, t)da-(\mu_v+\alpha_v^j )y_v^j( t),\\[4ex] \displaystyle \frac{dx_h(t)}{dt}=-\sum^n_{j=1}\beta_h^jS_{h_0}^*y_v^j( t)-\mu_h x_h(t),\\[4ex] \displaystyle \frac{\partial z_h^j(\tau, t)}{\partial \tau}+\frac{\partial z_h^j(\tau, t)}{\partial t}=-(\mu_h+m_h^j(\tau))z_h^j(\tau, t),\\[4ex] \displaystyle z_h^j(0, t)=\beta_h^jS_{h_0}^*y_v^j( t),\\[4ex] \displaystyle \frac{\partial y_h^j(a, t)}{\partial a}+\frac{\partial y_h^j(a, t)}{\partial t}=-(\mu_h+\alpha_h^j(a)+r_h^j(a))y_h^j(a, t),\\[4ex] \displaystyle y_h^j(0, t)=\int^\infty_0 m_h^j(\tau)z_h^j(\tau, t)d\tau. \end{array}$$ \right. $

(18)

Let

$ \label{eq:char} y_v^j{(t)}=\bar{y}_v^je^{\lambda t},\ z_h^j{(\tau, t)}=\bar{z}_h^j(\tau)e^{\lambda t}, \ y_h^j{(a, t)}=\bar{y}_h^j(a)e^{\lambda t}, $

(19)

where $\bar{y}_v^j, \bar{z}_h^j(\tau)$ and $\bar{y}_h^j(a)$ are to be determined. Substituting (19) into (18), we obtain

$ \left\{ $$\begin{split} & \displaystyle \lambda\bar{y}_v^j=S_{v_0}^*\int^\infty_0\beta_v^j(a)\bar{y}_h^j(a)da-(\mu_v+\alpha_v^j )\bar{y}_v^j,\\ & \displaystyle \frac{d \bar{z}_h^j(\tau)}{d \tau}=-(\lambda+\mu_h+m_h^j(\tau))\bar{z}_h^j(\tau),\\ & \displaystyle \bar{z}_h^j(0)=\beta_h^jS_{h_0}^*\bar{y}_v^j,\\ & \displaystyle \frac{d \bar{y}_h^j(a)}{d a}=-(\lambda+\mu_h+\alpha_h^j(a)+r_h^j(a))\bar{y}_h^j(a),\\ & \displaystyle \bar{y}_h^j(0)=\int^\infty_0m_h^j(\tau)\bar{z}_h^j(\tau)d\tau.\\ \end{split}$$ \right. $

(20)

Solving the differential equation, we obtain

$ $$\begin{align*} \bar{z}_h^j(\tau)=\bar{z}_h^j(0)~e^{-\lambda \tau}\pi^j_1(\tau)=\beta_h^jS_{h_0}^*\bar{y}_v^j ~e^{-\lambda \tau}\pi^j_1(\tau). \end{align*}$$ $

Substituting the expression for $\bar{z}_h^j(\tau)$ into the equation for $\bar{y}_h^j(0)$, expressing $\bar{y}_h^j(0)$ in term of $\bar{z}_h^j(0)$, and replacing $\bar{y}_h^j(0)$ in the equation for $\bar{y}_h^j(a)$, we obtain

$ $$\begin{align*} \bar{y}_h^j(a)=\bar{y}_h^j(0)~e^{-\lambda a}\pi^j_2(a)=\beta_h^jS_{h_0}^*\bar{y}_v^j~e^{-\lambda a}\pi^j_2(a) \int^\infty_0m_h^j(\tau)~e^{-\lambda \tau}\pi^j_1(\tau)d\tau. \end{align*}$$ $

Substituting the above expression for $\bar{y}_h^j(a)$ into the first equation of (20), we can obtain

$ $$\begin{array}{ll} \displaystyle \lambda+\mu_v+\alpha_v^j =\beta_h^jS_{v_0}^*S_{h_0}^*\int^\infty_0m_h^j(\tau)e^{-\lambda \tau}\pi^j_1(\tau)d\tau\int^\infty_0\beta_v^j(a)e^{-\lambda a}\pi^j_2(a)da. \end{array}$$ $

(21)

Now we are able to state the following result.

Theorem 3.1. If

$ \mathcal R_0=\max \{\mathcal R^1_0,\cdots,\mathcal R^n_0\}< 1, $

then the disease-free equilibrium is locally asymptotically stable. If $\mathcal R_0>1$, it is unstable.

Proof. We first prove the first result. Let us assume $\mathcal R_0 <1$. For ease of notation, set

$ $$\begin{split} \displaystyle LHS\stackrel{\it def}{=} & \lambda+\mu_v+\alpha_v^j,\\ RHS\stackrel{\it def}{=} & \mathcal{G}_1(\lambda)=\beta_h^jS_{v_0}^*S_{h_0}^*\int^\infty_0m_h^j(\tau)e^{-\lambda \tau}\pi^j_1(\tau)d\tau\int^\infty_0\beta_v^j(a)e^{-\lambda a}\pi^j_2(a)da. \end{split}$$ $

(22)

We can easily verify that

$ $$\begin{align*} \displaystyle |LHS| & \geq\displaystyle\mu_v+\alpha_v^j,\\ \displaystyle|RHS| & \displaystyle\leq\mathcal{G}_1(\Re\lambda)\leq\mathcal{G}_1(0)=\beta_h^jS_{v_0}^*S_{h_0}^*\int^\infty_0m_h^j(\tau)\pi^j_1(\tau)d\tau\int^\infty_0\beta_v^j(a)\pi^j_2(a)da\\ & =\displaystyle\frac{\beta_h^j\Lambda_v\Lambda_h}{\mu_v\mu_h}\int^\infty_0m_h^j(\tau)\pi^j_1(\tau)d\tau\int^\infty_0\beta_v^j(a)\pi^j_2(a)da\\[2ex] & \displaystyle=\mathcal R^j_0(\mu_v+\alpha_v^j)<|LHS|, \end{align*}$$ $

for any $\lambda, \Re \lambda\geq 0$. There is a contradiction. The contradiction implies that the equation (21) cannot have any roots with non-negative real parts. Hence, the disease-free equilibrium is locally asymptotically stable.

Next, let us assume $ \max \{\mathcal R^1_0,\cdots,\mathcal R^n_0\}=\mathcal R^{j_0}_0>1. $ We rewrite the characteristic equation (21) in the form

$ \displaystyle \mathcal{G}_2(\lambda)=0, $

(23)

where

$ $$\begin{array}{ll} & \displaystyle \mathcal{G}_2(\lambda)\\ = & \displaystyle(\lambda+\mu_v+\alpha_v^{j_0}) -\beta_h^{j_0}S_{v_0}^*S_{h_0}^*\int^\infty_0m_h^{j_0}(\tau)e^{-\lambda \tau}\pi^{j_0}_1(\tau)d\tau\int^\infty_0\beta_v^{j_0}(a)e^{-\lambda a}\pi^{j_0}_2(a)da. \end{array}$$ $

It is easily verified that

$ $$\begin{array}{ll} \displaystyle \mathcal{G}_2(0) & \displaystyle=(\mu_v+\alpha_v^{j_0}) -\beta_h^{j_0}S_{v_0}^*S_{h_0}^*\int^\infty_0m_h^{j_0}(\tau)\pi^{j_0}_1(\tau)d\tau\int^\infty_0\beta_v^{j_0}(a)\pi^{j_0}_2(a)da\\[2ex] & =\displaystyle (\mu_v+\alpha_v^{j_0})(1-\mathcal R^{j_0}_0)<0, \end{array}$$ $

and

$\mathop {\lim }\limits_{\lambda \to + \infty } {{\cal G}_2}(\lambda ) = + \infty .$

Hence, the characteristic equation (23) has a real positive root. Therefore, the disease free equilibrium $\mathcal{E}_0$ is unstable. This concludes the proof.

We have proved that the disease-free equilibrium is locally stable if $\mathcal R_0<1$. It also follows from Theorem 2.3 that strain $j$ will die out if $\mathcal R^j_0<1$. Therefore we have the following result.

Theorem 3.2. If

$\mathcal R_0=\max \{\mathcal R^1_0,\cdots,\mathcal R^n_0\}< 1,$

then the disease-free equilibrium $\mathcal{E}_0$ is globally asymptotically stable.

---

4. Existence and stability of boundary equilibria

In this section, we mainly investigate the existence and stability of the boundary equilibria. For ease of notation, let

$ $$\begin{array}{ll} \displaystyle \Delta_j=\frac{\beta_h^j\Lambda_h\Lambda_v}{\mu_h\mu_v(\mu_v+\alpha_v^j)},\\[2ex] \displaystyle b_j=\int^\infty_0m_h^j(\tau)\pi^j_1(\tau)d\tau\int^\infty_0\beta_v^j(a)\pi^j_2(a)da,\\[2ex] \displaystyle b_j(\lambda)=\int^\infty_0m_h^j(\tau)e^{-\lambda \tau}\pi^j_1(\tau)d\tau\int^\infty_0\beta_v^j(a)e^{-\lambda a}\pi^j_2(a)da. \end{array}$$ $

(24)

From Theorem 2.3, it follows that strain $j$ will die out if $\mathcal R_0^j<1$. Thus in later sections we always assume that $\mathcal R_0^j>1$ for all $j,j=1,2,\cdots,n$. If $\mathcal R_0^j>1$, straightforward computation yields that system (2) has a corresponding single-strain equilibrium $\mathcal{E}_j$ which is given by

$ \mathcal{E}_j=(S_v^{j*},0,\cdots,0, I_v^{j*}, 0,\cdots,0,S_h^{j*}, 0,\cdots,0,E_h^{j*}(\tau), I_h^{j*}(a),0,\cdots,0). $

The non-zero components $I_v^{j*},E_h^{j*}$ and $I_h^{j*}$ are in positions $j+1,n+2j+1$ and $n+2j+2$, respectively, and

$ $$\begin{split} & \displaystyle I_v^{j*}= \displaystyle \frac{\mu_v\mu_h(\mathcal R^j_0-1)}{\beta_h^j(\Lambda_hb_j+\mu_v)},\\ & \displaystyle S_v^{j*}= \displaystyle\frac{\Lambda_v-(\mu_v+\alpha_v^j)I_v^{j*}}{\mu_v}=\frac{\beta_h^j\Lambda_v(\mu_v+\Lambda_hb_j)-\mu_v\mu_h(\mu_v+\alpha_v^j)(\mathcal R^j_0-1)}{\beta_h^j\mu_v(\mu_v+\Lambda_hb_j)}, \end{split}$$ $

$ $$\begin{split} & \displaystyle S_h^{j*}= \displaystyle\frac{\Lambda_h}{\beta_h^jI_v^{j*}+\mu_h}=\frac{\Lambda_h(\mu_v+\Lambda_hb_j)}{\mu_h(\mu_v\mathcal R^j_0+\Lambda_hb_j)},\\[2ex] & \displaystyle E_h^{j*}(\tau)= \displaystyle E_h^{j*}(0)\pi^j_1(\tau), E_h^{j*}(0)= \beta_h^jS_h^{j*}I_v^{j*}, \\ & \displaystyle I_h^{j*}(a)= \displaystyle I_h^{j*}(0)\pi^j_2(a), I_h^{j*}(0)=E_h^{j*}(0) \int^\infty_0m_h^j(\tau)\pi^j_1(\tau)d\tau. \end{split}$$ $

(25)

The results on the local stability of single-strain equilibrium $\mathcal{E}_{j_0}$ are summarized below:

Theorem 4.1. Assume $\mathcal{R}^{j_0}_0 > 1$ for a fixed $j_0$ and

$ \mathcal{R}^j_0 < \mathcal{R}^{j_0}_0 \;for\;all \; j\neq j_0. $

Then single-strain equilibrium $\mathcal{E}_{j_0}$ is locally asymptotically stable. If there exists $i_0$ such that

$ \mathcal{R}^{i_0}_0 > \mathcal{R}^{j_0}_0, $

then the single-strain equilibrium $\mathcal{E}_{j_0}$ is unstable.

Proof. Without loss of generality, we assume that $\mathcal{R}^1_0 > 1$ and $ \mathcal{R}^i_0 < \mathcal{R}^1_0$ for $i= 2, \cdots, n$. Let

$ $$\begin{array}{ll} \displaystyle S_v(t)=S_v^{1^*}+x_v(t),\ S_h(t) =S_h^{1^*} +x_h(t),\\[2ex] \displaystyle I_v^1(t) =I_v^{1^*} +y_v^1(t), \ E_h^1{(\tau, t)} =E_h^{1^*}(\tau)+ z_h^1{(\tau, t)}, \ I_h^1{(a,t)} =I_h^{1^*}(a)+ y_h^1{(a, t)},\\[2ex] \displaystyle I_v^i(t)=y_v^i(t), E_h^i{(\tau, t)} =z_h^i{(\tau, t)}, \ I_h^i{(a, t)}= y_h^i{(a, t)}, \end{array}$$ $

where $ i=2,\cdots,n.$ Then the linearization system of system (2) at the equilibrium $\mathcal{E}_{1}$ can be expressed as

$\left\{ $$\begin{split} & \displaystyle \frac{dx_v(t)}{dt}=-S_v^{1^*}\int^\infty_0\beta_v^1(a)y_h^1(a, t)da-x_v(t)\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da\\ & \displaystyle -\sum^n_{i=2}S_v^{1^*}\int^\infty_0\beta_v^i(a)y_h^i(a, t)da-\mu_v x_v(t),\\ & \displaystyle \frac{dy_v^1(t)}{dt}=S_v^{1*}\int^\infty_0\beta_v^1(a)y_h^1(a, t)da+x_v(t)\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da-(\mu_v+\alpha_v^1) y_v^1(t),\\ & \displaystyle \frac{dy_v^i(t)}{dt}=S_v^{1^*}\int^\infty_0\beta_v^i(a)y_h^i(a, t)da-(\mu_v+\alpha_v^i) y_v^i(t),\\ & \displaystyle \frac{dx_h(t)}{dt}=-\beta_h^1S_h^{1^*}y_v^1(t)-\beta_h^1x_h(t)I_v^{1^*}-\sum^n_{i=2}\beta_h^iS_h^{1^*}y_v^i(t)-\mu_h x_h(t),\\ & \displaystyle \frac{\partial z_h^j(\tau, t)}{\partial \tau}+\frac{\partial z_h^j(\tau, t)}{\partial t}=-(\mu_h+m_h^j (\tau))z_h^j(\tau, t),\\ & \displaystyle z_h^1(0, t)=\beta_h^1S_h^{1^*}y_v^1(t)+\beta_h^1x_h(t)I_v^{1^*},\\ & \displaystyle z_h^i(0, t)=\beta_h^iS_h^{1^*}y_v^i(t),\\ & \displaystyle \frac{\partial y_h^j(a, t)}{\partial a}+\frac{\partial y_h^j(a, t)}{\partial t}=-(\mu_h+\alpha_h^j(a)+r_h^j(a))y_h^j(a, t),\\ & \displaystyle y_h^j(0, t)=\int^\infty_0m_h^j(\tau)z_h^j(\tau,t)d\tau.\\ \end{split}$$ \right. $

(26)

An approach similar to [14] (see Appendix B in [14]) can show that the linear stability of the system is determined by the eigenvalues of the linearized system (26). In order to investigate the linear stability of the linearized system (26), we consider exponential solutions (see the case of the disease-free equilibrium) and obtain a linear eigenvalue problem. For the whole system, we only consider the equations for strains $i, i=2, \cdots, n$, and obtain the following eigenvalue problem:

$\left\{ $$\begin{array}{ll} \displaystyle \frac{dy_v^i(t)}{dt}=S_v^{1^*}\int^\infty_0\beta_v^i(a)y_h^i(a, t)da-(\mu_v+\alpha_v^i) y_v^i(t),\\[2ex] \displaystyle \frac{\partial z_h^i(\tau, t)}{\partial \tau}+\frac{\partial z_h^i(\tau, t)}{\partial t}=-(\mu_h+m_h^i (\tau))z_h^i(\tau, t),\\[2ex] \displaystyle z_h^i(0, t)=\beta_h^iS_h^{1^*}y_v^i(t),\\[2ex] \displaystyle \frac{\partial y_h^i(a, t)}{\partial a}+\frac{\partial y_h^i(a, t)}{\partial t}=-(\mu_h+\alpha_h^i(a)+r_h^i(a))y_h^i(a, t),\\[2ex] \displaystyle y_h^i(0, t)=\int^\infty_0m_h^i(\tau)z_h^i(\tau,t)d\tau.\\[2ex] \end{array}$$ \right. $

(27)

For each $i, i\neq 1$, by using the same argument to equation (21), we obtain the following characteristic equation

$ $$\begin{array}{ll} \displaystyle \lambda+\mu_v+\alpha_v^i =\beta_h^iS_v^{1*}S_h^{1^*}\int^\infty_0m_h^i(\tau)e^{-\lambda \tau}\pi^i_1(\tau)d\tau\int^\infty_0\beta_v^i(a)e^{-\lambda a}\pi^i_2(a)da. \end{array}$$ $

(28)

Notice that $S_v^{j*}$ and $S_h^{j*}$ satisfy

$ \displaystyle \beta_h^jS_v^{j*}S_h^{j*}\int^\infty_0m_h^j(\tau)\pi^j_1(\tau)d\tau\int^\infty_0\beta_v^j(a)\pi^j_2(a)da=\mu_v+\alpha_v^j, $

(29)

for $j=1,\cdots, n.$ It then follows from (6) and (24) that we have

$ \displaystyle S_v^{1^*}S_h^{1^*}=\frac{\mu_v+\alpha_v^1}{\beta_h^1b_1}=\frac{\Lambda_v\Lambda_h}{\mu_v\mu_h\mathcal R^1_0}. $

(30)

Substituting (30) into the equation (28), we get the following characteristic equation

$ $$\begin{array}{ll} \displaystyle \lambda+\mu_v+\alpha_v^i =\beta_h^i\frac{\Lambda_v\Lambda_h}{\mu_v\mu_h\mathcal R^1_0}b_i(\lambda), \end{array}$$ $

(31)

where $b_i(\lambda)$ is defined in (24).

First, assume that $ \mathcal{R}^{i_0}_0 > \mathcal{R}^1_0$ for some $i_0$, and set

$\displaystyle \mathcal{G}_{i_0}(\lambda)\stackrel{\it def}{=}(\lambda+\mu_v+\alpha_v^{i_0}) -\beta_h^{i_0}\frac{\Lambda_v\Lambda_h}{\mu_v\mu_h\mathcal R^1_0}b_{i_0}(\lambda).$

Straightforward computation yields that

$ \displaystyle \mathcal{G}_{i_0}(0)=(\mu_v+\alpha_v^{i_0}) -\beta_h^{i_0}\frac{\Lambda_v\Lambda_h}{\mu_v\mu_h\mathcal R^1_0}b_{i_0}=(\mu_v+\alpha_v^{i_0})(1-\frac{R^{i_0}_0}{R^1_0})<0. $

Furthermore, for $\lambda$ real, $\mathcal{G}_{i_0}(\lambda)$ is an increasing function of $\lambda$ such that $\lim \mathcal{G}_{i_0}(\lambda)\rightarrow +\infty$ as $\lambda\rightarrow+\infty$. Hence Intermediate Value Theorem implies that the equation (31) has a unique real positive solution. We conclude that in that case $\mathcal{E}_1$ is unstable.

Next, assume $\mathcal{R}^i_0<\mathcal{R}^1_0$ for all $i=2,\cdots,n$, and set

$ \displaystyle \mathcal{G}_3(\lambda)\stackrel{\it def}{=}\lambda+\mu_v+\alpha_v^i, \mathcal{G}_4(\lambda)\stackrel{\it def}{=}\displaystyle \beta_h^i\frac{\Lambda_v\Lambda_h}{\mu_v\mu_h\mathcal R^1_0}b_i(\lambda). $

(32)

Consider $\lambda$ with $\Re\lambda \geq0$. For such $\lambda$, following from (32), we have

$ $$\begin{array}{ll} \displaystyle |\mathcal{G}_3(\lambda)| & \geq\displaystyle\mu_v+\alpha_v^i,\\ \displaystyle|\mathcal{G}_4(\lambda)| & \displaystyle\leq\mathcal{G}_4(\Re\lambda)\leq\mathcal{G}_4(0) =\frac{1}{\mathcal R^1_0}\beta_h^i\frac{\Lambda_v\Lambda_h}{\mu_v\mu_h}\int^\infty_0m_h^i(\tau)\pi^i_1(\tau)d\tau\int^\infty_0\beta_v^i(a)\pi^i_2(a)da\\[2ex] & \displaystyle=\frac{\mathcal R^i_0}{\mathcal R^1_0}(\mu_v+\alpha_v^i)<|\mathcal{G}_3(\lambda)|. \end{array}$$ $

This gives a contradiction. Hence, the equation (31) have no solutions with positive real part and all eigenvalues of these equations have negative real parts. Therefore, the stability of $\mathcal{E}_1$ depends on the eigenvalues of the following system

$\left\{ $$\begin{array}{ll} \displaystyle \lambda x_v=-S_v^{1^*}\int^\infty_0\beta_v^1(a)y_h^1(a)da-x_v\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da-\mu_v x_v,\\[2ex] \displaystyle \lambda y_v^1=S_v^{1^*}\int^\infty_0\beta_v^1(a)y_h^1(a)da+x_v\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da-(\mu_v+\alpha_v^1) y_v^1,\\[2ex] \displaystyle \lambda x_h=- z_h^1(0)-\mu_h x_h,\\[2ex] \displaystyle \frac{d z_h^1(\tau)}{d \tau}=-(\lambda+\mu_h+m_h^1(\tau))z_h^1(\tau),\\[2ex] \displaystyle z_h^1(0)=\beta_h^1S_h^{1^*}y_v^1+\beta_h^1I_v^{1*}x_h,\\[2ex] \displaystyle \frac{d y_h^1(a)}{d a}=-(\lambda+\mu_h+\alpha_h^1(a)+r_h^1(a))y_h^1(a),\\[2ex] \displaystyle y_h^1(0)=\int^\infty_0m_h^1(\tau)z_h^1(\tau)d\tau.\\[2ex] \end{array}$$ \right. $

(33)

Solving the differential equation, we have

$ $$\begin{array}{ll} \displaystyle z_h^1(\tau)=z_h^1(0)~e^{-\lambda \tau}\pi^1_1(\tau),\\ \displaystyle y_h^1(a)=y_h^1(0)~e^{-\lambda a}\pi^1_2(a)=z_h^1(0)~e^{-\lambda a}\pi^1_2(a)\int^\infty_0m_h^1(\tau)~e^{-\lambda \tau}\pi^1_1(\tau)d\tau. \end{array}$$ $

Substituting the above expression for $y_h^1(a)$ into the first and the second equations of (33) yileds that

$\left\{ $$\begin{array}{ll} \displaystyle (\lambda+\mu_v +\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da)x_v+S_v^{1^*}b_1(\lambda)z_h^1(0)=0,\\[2ex] \displaystyle -x_v\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da+(\lambda+\mu_v+\alpha_v^1) y_v^1-S_v^{1^*}b_1(\lambda)z_h^1(0)=0,\\[2ex] \displaystyle (\lambda+\mu_h)x_h+z_h^1(0)=0,\\[2ex] \displaystyle -\beta_h^1I_v^{1^*}x_h-\beta_h^1S_h^{1^*}y_v^1+z_h^1(0)=0.\\[2ex] \end{array}$$ \right. $

(34)

Direct calculation yields the following characteristic equation

$ $$\begin{split} \displaystyle (\lambda+\mu_v +\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da)(\lambda+\mu_v+\alpha_v^1) & \displaystyle(\lambda+\mu_h+\beta_h^1I_v^{1^*})\\ = & \displaystyle\beta_h^1S_h^{1^*}S_v^{1^*}b_1(\lambda)(\lambda+\mu_v)(\lambda+\mu_h). \end{split}$$ $

(35)

Dividing both sides by $(\lambda+\mu_v)(\lambda+\mu_h)$ gives

$ \displaystyle \mathcal{G}_5(\lambda)=\mathcal{G}_6(\lambda), $

(36)

where

$ $$\begin{array}{ll} \displaystyle \mathcal{G}_5(\lambda) & \displaystyle=\frac{(\lambda+\mu_v +\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da)(\lambda+\mu_v+\alpha_v^1) (\lambda+\mu_h+\beta_h^1I_v^{1^*})}{(\lambda+\mu_v)(\lambda+\mu_h)},\\[2ex] \displaystyle \mathcal{G}_6(\lambda) & \displaystyle=\beta_h^1S_h^{1^*}S_v^{1^*}b_1(\lambda)\\[2ex] & \displaystyle=\beta_h^1S_h^{1^*}S_v^{1^*}\int^\infty_0m_h^1(\tau)e^{-\lambda \tau}\pi^1_1(\tau)d\tau\int^\infty_0\beta_v^1(a)e^{-\lambda a}\pi^1_2(a)da. \end{array}$$ $

(37)

If $\lambda$ is a root with $\Re\lambda\geq0$, it follows from equation (37) that

$ $$\begin{array}{ll} \displaystyle |\mathcal{G}_5(\lambda)|>|\lambda+\mu_v+\alpha_v^1|\geq\mu_v+\alpha_v^1. \end{array}$$ $

(38)

From (29), we have

$ $$\begin{split} \displaystyle |\mathcal{G}_6(\lambda)|\leq|\mathcal{G}_6(\Re\lambda)|\leq\mathcal{G}_6(0) & \displaystyle=\beta_h^1S_h^{1^*}S_v^{1^*}\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau\int^\infty_0\beta_v^1(a)\pi^1_2(a)da\\[2ex] & =\displaystyle \mu_v+\alpha_v^1<|\mathcal{G}_5(\lambda)|. \end{split}$$ $

(39)

This leads to a contradiction. The contradiction implies that (36) has no roots such that $\Re\lambda\geq0$. Thus, the characteristic equation for strain one has only roots with negative real parts. Thus, the single strain equilibrium $\mathcal{E}_1$ is locally asymptotically stable if $\mathcal{R}^1_0>1$ and $\mathcal{R}^i_0<\mathcal{R}^1_0,\ i=2,\cdots,n$. This concludes the proof.

---

5. Preliminary results and uniform persistence

In the previous section, we proved that if the disease reproduction number is less than one, all strains are eliminated and the disease dies out. Our next step is to show that the competitive exclusion principle holds for system (2). In the later sections, we always assume that $\mathcal R_0 >1$. Without loss of generality, we assume that

$\displaystyle \mathcal R^1_0=\max \{\mathcal R^1_0,\cdots,\mathcal R^n_0\}>1.$

In the following we will show that strain $1$ persists while the other strains die out if $\mathcal R^i_0/\mathcal R^1_0 < b_i/b_1 < 1,i\neq1$, where $b_j$ denotes the probability of a given susceptible vector being transmitted by an infected host with strain $j$. Hence, the strain with the maximal reproduction number eliminates all the rest and the competitive exclusion principle will be established for system (2).

Mathematically speaking, establishing the competitive exclusion principle means establishing the global stability of the single-strain equilibrium $\mathcal{E}_1$. From Theorem 4.1 we know that if $\mathcal R^i_0/\mathcal R^1_0 < 1,i\neq 1,$ the equilibrium $\mathcal{E}_1$ is locally asymptotically stable. In the following we only need to show that $\mathcal{E}_1$ is a global attractor. The method used here to show this result is similar to the one used in [1,12,13,20].

Set

$f(x)=x-1-\ln x.$

It is easy to check that $f(x)\geq0$ for all $x>0$ and $f(x)$ reaches its global minimum value $f(1)=0$ when $x=1$. Next, let us define the following Lyapunov function

$ $$\begin{array}{ll} \displaystyle U(t)=\displaystyle U_1(t)+U^1_2(t)+\sum^n_{i=2}U^i_2(t)+U_3(t)+U^1_4(t)+\sum^n_{i=2}U^i_4(t)+U^1_5(t)+\sum^n_{i=2}U^i_5(t), \end{array}$$ $

(40)

where

$ \left. $$\begin{array}{ll} \displaystyle U_1(t)=\frac{1}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} f\bigg(\frac{S_v}{S_v^{1^*}}\bigg),\\ \displaystyle U^1_2(t)=\frac{1}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} I_v^{1^*}f\bigg(\frac{I_v^1}{I_v^{1^*}}\bigg) ,\\[2ex] \displaystyle U^i_2(t)=\frac{1}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}I_v^i ,\\[2ex] \displaystyle U_3(t)=S_h^{1^*} f\bigg(\frac{S_h}{S_h^{1^*}}\bigg) ,\\[2ex] \displaystyle U^1_4(t)=\frac{1}{\mathcal R^1_0}\int^\infty_0p_1(\tau)E_h^{1^*}(\tau)f\bigg(\frac{E_h^1(\tau,t)}{E_h^{1^*}(\tau)}\bigg)d\tau,\\[2ex] \displaystyle U^i_4(t)=\frac{1}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \int^\infty_0p_i(\tau)E_h^i(\tau,t)d\tau,\\[2ex] \displaystyle U^1_5(t)=\frac{1}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \int^\infty_0 q_1(a)I_h^{1^*}(a)f\bigg(\frac{I_h^1(a,t)}{I_h^{1^*}(a)}\bigg)da.\\[2ex] \displaystyle U^i_5(t)=\frac{1}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \int^\infty_0 q_i(a)I_h^i(a,t)da,\\[2ex] \end{array}$$ \right. $

(41)

and

$ $$\begin{array}{ll} \displaystyle q_j(a)=\int^\infty_a\beta_v^j(s)e^{-\int^s_a(\mu_h+\alpha_h^j(\sigma)+r_h^j(\sigma))d\sigma}ds,\\[2ex] \displaystyle p_j(\tau)=\Delta_jq_j(0)\int^\infty_{\tau}m_h^j(s)e^{-\int^s_{\tau}(\mu_h+m_h^j(\sigma))d\sigma}ds. \end{array}$$ $

(42)

Direct computation gives

$ \displaystyle p_j(0)=\mathcal R^j_0, $

and

$ $$\begin{array}{ll} \displaystyle q'_j(a)=-\beta_v^j(a)+ (\mu_h+\alpha_h^j(a)+r_h^j(a))q_j(a),\\[2ex] \displaystyle p'_j(\tau)=-\Delta_jq_j(0)m_h^j(\tau)+(\mu_h+m_h^j(\tau))p_j(\tau). \end{array}$$ $

(43)

The main difficulty with the Lyapunov function $U$ above is that the Lyapunov function $U$ is well defined. Thus in the following we first show that strain one persists both in the hosts and in the vectors as the other strains die out. Let

$ \displaystyle \hat{X }_1=\bigg\{\varphi_1\in L^1_+(0, \infty)\bigg|\exists s\geq0:\ \int^\infty_0 m_h^1(\tau+s)\varphi_1(\tau)d\tau>0\bigg\}, $

$ \displaystyle \hat{X }_2=\bigg\{\psi_1\in L^1_+(0, \infty)\bigg|\exists s\geq0:\ \int^\infty_0 \beta_v^1(a+s)\psi_1(a)da>0\bigg\}, $

and define

$ \displaystyle X_0=\mathbb{R}_+\times \prod^n_{j=1}\mathbb{R}_+\times \mathbb{R}_+\times \hat{X }_1\times \hat{X }_2\times \prod^n_{i=2}(L^1(0,\infty)\times L^1(0,\infty)), $

$ \Omega_0=\Omega\cap X_0.$

Note that $\Omega_0$ is forward invariant. This is because (3) show that $\Omega$ is forward invariant. To see $X_0$ is forward invariant, we firstly demonstrate that $\hat{X }_2$ is forward invariant. Let us assume that the inequality holds for the initial condition. The inequality says that the support of $\beta_v^1(a)$ will intersect the support of the initial condition if it is transferred $s$ units to the right. Since the support of the initial condition only moves to the right, the intersection will take place for any other time if that happens for the initial time. Similarly, $\hat{X }_1$ is also forward invariant. Therefore, $\Omega_0$ is forward invariant.

Now let us recall two important definitions.

Definition 5.1. Strain one is called uniformly weakly persistence if there exists some $\gamma>0$ independent of the initial conditions such that

$ \displaystyle \limsup\limits_{t\rightarrow\infty}\int^\infty_0E_h^1(\tau, t)d\tau > \gamma \text{whenever}\int^\infty_0 \varphi_1(\tau)d\tau > 0, $

$ \displaystyle \limsup\limits_{t\rightarrow\infty}\int^\infty_0I_h^1(a,t)da > \gamma \text{whenever} \int^\infty_0\psi_1(a)da> 0, $

and

$ \displaystyle \limsup\limits_{t\rightarrow\infty}I_v^1(t) > \gamma \text{whenever} I_{v_0}^1> 0, $

for all solutions of system (2).

One of the important implications of uniform weak persistence of the disease is that the disease-free equilibrium is unstable.

Definition 5.2. Strain one is uniformly strongly persistence if there exists some $\gamma>0$ independent of the initial conditions such that

$ \displaystyle \liminf\limits_{t\rightarrow\infty}\int^\infty_0E_h^1(\tau, t)d\tau > \gamma \text{whenever}\int^\infty_0\varphi_1(\tau)d\tau> 0, $

$ \displaystyle \liminf\limits_{t\rightarrow\infty}\int^\infty_0I_h^1(a,t)da > \gamma \text{whenever} \int^\infty_0\psi_1(a)da > 0, $

and

$ \displaystyle \liminf\limits_{t\rightarrow\infty}I_v^1(t) > \gamma \text{whenever} I_{v_0}^1 > 0, $

for all solutions of model (2).

It is evident from the definitions that, if the disease is uniformly strongly persistent, it is also uniformly weakly persistent.

Now we are able to state the main results in this section.

Theorem 5.3. Assume $\mathcal R^1_0 >1$ and $\mathcal R^i_0 <\mathcal R^1_0$ for $i=2, \cdots, n$. Furthermore, assume that the other strains except stain 1 will die out, i.e.,

$ \displaystyle\limsup\limits_{t\rightarrow +\infty}I_v^i(t)=0,\ \limsup\limits_{t\rightarrow +\infty}\int^\infty_0 E_h^i(\tau, t)d\tau=0 \ \;and\;\ \limsup\limits_{t\rightarrow +\infty}\int^\infty_0I_h^i(a, t)da=0, $

for $ i=2, \cdots, n$. Then strain $1$ is uniformly weakly persistent for the initial conditions that belong to $\Omega_0$, i.e., there exists $\gamma>0$ such that

$ \displaystyle \limsup\limits_{t\rightarrow+\infty}\beta_h^1I_v^1(t)\geq\gamma, \limsup\limits_{t\rightarrow+\infty}\int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau\geq\gamma, \limsup\limits_{t\rightarrow+\infty}\int^\infty_0 \beta_v^1(a)I_h^1(a, t)da\geq\gamma. $ }

Proof. We argue by contradiction. Assume that strain $1$ also dies out. For any $\varepsilon>0$ and every initial condition in $\Omega_0$ such that

$ \displaystyle \limsup\limits_{t\rightarrow+\infty}\beta_h^1I_v^1(t) < \varepsilon, \limsup\limits_{t\rightarrow+\infty}\int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau < \varepsilon, \limsup\limits_{t\rightarrow+\infty}\int^\infty_0 \beta_v^1(a)I_h^1(a, t)da<\varepsilon. $

Following that there exist $T>0$ such that for all $t>T$ we have

$ \displaystyle \beta_h^jI_v^j(t)<\varepsilon,\ \int^\infty_0 m_h^j(\tau)E_h^j(\tau, t)d\tau<\varepsilon,\ \int^\infty_0 \beta_v^j(a)I_h^j(a, t)da<\varepsilon,\ j=1,\cdots,n. $

We may assume that the above inequality holds for all $t\geq0$ by shifting the dynamical system. From the first equation in (2) we have

$ \displaystyle S_v'(t) \geq\Lambda_v-n\varepsilon S_v-\mu_v S_v, S_h'(t) \geq\Lambda_h-n\varepsilon S_h-\mu_h S_h. $

Exploiting the comparison principle, we have

$ \displaystyle \limsup\limits_{t\rightarrow+\infty}S_v(t)\geq \liminf\limits_{t\rightarrow+\infty}S_v(t)\geq\frac{\Lambda_v}{n\varepsilon+\mu_v},\ \limsup\limits_{t\rightarrow+\infty}S_h(t)\geq \liminf\limits_{t\rightarrow+\infty}S_h(t)\geq\frac{\Lambda_h}{n\varepsilon+\mu_h}. $

Since $B_ E^1(t)=E_h^1(0, t),\ B_ I^1(t)=I_h^1(0, t) $, it then follows from system (2) that

$ \left\{ $$\begin{array}{ll} \displaystyle B_ E^1(t)=E_h^1(0, t)=\beta_h^1S_hI_v^1(t)\geq\beta_h^1\frac{\Lambda_h}{n\varepsilon+\mu_h}I_v^1(t),\\[2ex] \displaystyle \frac{dI_v^1(t)}{dt}\geq \frac{\Lambda_v}{n\varepsilon+\mu_v}\int^\infty_0 \beta_v^1(a)I_h^1(a, t)da-(\mu_v+\alpha_v^1)I_v^1(t). \end{array}$$ \right. $

(44)

By using the equations in (7), we can easily obtain the following inequalities on $B_ E^1(t),\ B_ I^1(t)$ and $I_v^1(t)$:

$ \left\{ $$\begin{array}{ll} \displaystyle B_ E^1(t)\geq\beta_h^1\frac{\Lambda_h}{n\varepsilon+\mu_h}I_v^1(t),\\[2ex] \displaystyle B_ I^1(t)=\int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau\geq\int^t_0 m_h^1(\tau)B_ E^1( t-\tau)\pi^1_1(\tau)d\tau,\\[2ex] \displaystyle \frac{dI_v^1(t)}{dt}\geq \frac{\Lambda_v}{n\varepsilon+\mu_v}\int^t_0 \beta_v^1(a)B_ I^1( t-a)\pi^1_2(a)da-(\mu_v+\alpha_v^1)I_v^1(t). \end{array}$$ \right. $

(45)

Let us take the Laplace transform of both sides of inequalities (45). Since all functions above are bounded, the Laplace transforms of the functions exist for $\lambda>0$. Denote the Laplace transforms of the functions $ B_E^1(t)$, $ B_I^1(t)$ and $ I_v^1(t)$ by $\hat{B}_E^1(\lambda)$, $\hat{B}_I^1(\lambda)$ and $\hat{I}_v^1(\lambda)$, respectively. Furthermore, set

$ \displaystyle \hat{K}_1(\lambda)=\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)e^{-\lambda\tau}d\tau, \displaystyle \hat{K}_2(\lambda)=\int^\infty_0\beta_v^1(a)\pi^1_2(a)e^{-\lambda a}da. $

(46)

Using the convolution property of the Laplace transform, we obtain the following inequalities for $\hat{B}_E^1(\lambda),\ \hat{B}_I^1(\lambda)$ and $\hat{I}_v^1(\lambda)$:

$ \left\{ $$\begin{array}{ll} \displaystyle\hat{ B}_E^1(\lambda)\geq\beta_h^1\frac{\Lambda_h}{n\varepsilon+\mu_h}\hat{I}_v^1(\lambda),\\[2ex] \displaystyle \hat{B}_I^1(\lambda)\geq\hat{K}_1(\lambda)\hat{ B}_E^1(\lambda),\\[2ex] \displaystyle \lambda \hat{I}_v^1(\lambda)-I_v^1(0)\geq \frac{\Lambda_v}{n\varepsilon+\mu_v}\hat{K}_2(\lambda)\hat{B}_I^1(\lambda)-(\mu_v+\alpha_v^1)\hat{I}_v^1(\lambda). \end{array}$$ \right. $

(47)

Eliminating $\hat{B}_I^1(\lambda)$ and $ \hat{I}_v^1(\lambda)$ yields

$ $$\begin{array}{ll} \displaystyle \hat{B}_E^1(\lambda)\geq\frac{\beta_h^1\Lambda_v\Lambda_h\hat{K}_1(\lambda)\hat{K}_2(\lambda)} {(n\varepsilon+\mu_v)(n\varepsilon+\mu_h)(\lambda+\mu_v+\alpha_v^1)}\hat{B}_E^1(\lambda)\displaystyle+\frac{\beta_h^1\Lambda_h} {(n\varepsilon+\mu_h)(\lambda+\mu_v+\alpha_v^1)}I_v^1(0). \end{array}$$ $

(48)

This is impossible since

$ \displaystyle \frac{\beta_h^1\Lambda_v\Lambda_h\hat{K}_1(0)\hat{K}_2(0)} {\mu_v\mu_h(\mu_v+\alpha_v^1)}:=\mathcal R^1_0>1, $

we can choose $\varepsilon$ and $\lambda$ small enough such that

$ \displaystyle \frac{\beta_h^1\Lambda_v\Lambda_h\hat{K}_1(\lambda)\hat{K}_2(\lambda)} {(n\varepsilon+\mu_v)(n\varepsilon+\mu_h)(\lambda+\mu_v+\alpha_v^1)}>1. $

The contradiction implies that there exists $\gamma>0$ such that for any initial condition in $\Omega_0$, we have

$ \displaystyle \limsup\limits_{t\rightarrow+\infty}\beta_h^1I_v^1(t)\geq\gamma, \limsup\limits_{t\rightarrow+\infty}\int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau\geq\gamma, \limsup\limits_{t\rightarrow+\infty}\int^\infty_0 \beta_v^1(a)I_h^1(a, t)da\geq\gamma. $

In addition, the equation for $I_v^1$ can be rewritten in the form

$ \displaystyle \frac{dI_v^1}{dt}\geq\frac{\Lambda_v\gamma}{n\gamma+\mu_v}-(\mu_v+\alpha_v^1)I_v^1,$

which implies a lower bound for $I_v^1$. This concludes the proof.

Next, we claim that system (2) has a global compact attractor $\mathfrak{T}$. Firstly, define the semiflow $\Psi: [0, \infty)\times \Omega_0\rightarrow \Omega_0$ generated by the solutions of system (2)

$ $$\begin{array}{ll} \displaystyle \Psi\bigg(t; S_{v_0}, I_{v_0}^1,\cdots,I_{v_0}^n, S_{h_0}, \varphi_1(\cdot), \psi_1(\cdot),\cdots,\varphi_n(\cdot), \psi_n(\cdot)\bigg)\\[2ex] =\displaystyle \bigg( S_v(t), I_v^1(t),\cdots,I_v^n(t),S_h(t), E_h^1(\tau, t),I_h^1(a,t), \cdots, E_h^n(\tau, t),I_h^n(a,t)\bigg). \end{array}$$ $

Definition 5.4. A set $\mathfrak{T}$ in $\Omega_0$ is called a global compact attractor for $\Psi$ if $\mathfrak{T}$ is a maximal compact invariant set and for all open sets $\mathfrak{U}$ containing $\mathfrak{T}$ and all bounded sets $\mathcal{B}$ of $\Omega_0$ there exists some $T>0$ such that $\Psi(t, \mathcal{B})\subseteq\mathfrak{U}$ holds for $t>T$.

Theorem 5.5. Under the hypothesis of Theorem 5.3, there exists $\mathfrak{T}$, a compact subset of $\Omega_0$, which is a global attractor for the semiflow $\Psi$ on $\Omega_0$. Moreover, we have

$ \displaystyle \Psi(t, x^0)\subseteq\mathfrak{T} \text{for every } x^0\in\mathfrak{T},\ \forall t\geq0. $

Proof. We split the solution semiflow into two components. For an initial condition $x^0\in \Omega_0$, }let $\Psi(t, x^0)=\hat{\Psi}(t, x^0)+\tilde{\Psi}(t, x^0)$, where

$ $$\begin{array}{ll} \displaystyle \hat{\Psi}\bigg(t; S_{v_0}, I_{v_0}^1,\cdots,I_{v_0}^n, S_{h_0}, \varphi_1(\cdot), \psi_1(\cdot),\cdots,\varphi_n(\cdot), \psi_n(\cdot)\bigg)\\[2ex] =\displaystyle \bigg( 0, 0,\cdots,0,0, \hat{E}_h^1(\tau, t),\hat{I}_h^1(a,t), \cdots, \hat{E}_h^n(\tau, t),\hat{I}_h^n(a,t)\bigg), \end{array}$$ $

(49)

$ $$\begin{array}{ll} \displaystyle \tilde{\Psi}\bigg(t; S_{v_0}, I_{v_0}^1,\cdots,I_{v_0}^n, S_{h_0}, \varphi_1(\cdot), \psi_1(\cdot),\cdots,\varphi_n(\cdot), \psi_n(\cdot)\bigg)\\[2ex] =\displaystyle \bigg( S_v(t), I_v^1(t),\cdots,I_v^n(t),S_h(t), \tilde{E}_h^1(\tau, t),\tilde{I}_h^1(a,t), \cdots, \tilde{E}_h^n(\tau, t),\tilde{I}_h^n(a,t)\bigg), \end{array}$$ $

(50)

and $\displaystyle E_h^j(\tau, t)=\hat{E}_h^j(\tau, t)+\tilde{E}_h^j(\tau, t),\ I_h^j(a, t)=\hat{I}_h^j(a, t)+\tilde{I}_h^j(a, t)$ for $j=1,\cdots,n$. $\hat{E}_h^j(\tau, t),$ $\hat{I}_h^j(a, t)$, $\tilde{E}_h^j(\tau, t)$, $\tilde{I}_h^j (a, t)$ are the solutions of the following equations

$ \left\{ $$\begin{array}{ll} \displaystyle \frac{\partial \hat{E}_h^j}{\partial t}+\frac{\partial \hat{E}_h^j}{\partial \tau}=-(\mu_h+m_h^j(\tau))\hat{E}_h^j(\tau, t),\\[2ex] \displaystyle \hat{E}_h^j(0, t)=0,\\[2ex] \displaystyle \hat{E}_h^j(\tau, 0)=\varphi_j(\tau), \end{array}$$ \right. $

(51)

$ \left\{ $$\begin{array}{ll} \displaystyle \frac{\partial \hat{I}_h^j}{\partial t}+\frac{\partial \hat{I}_h^j}{\partial a}=-(\mu_h+\alpha_h^j(a)+r_h^j(a))\hat{I}_h^j(a, t),\\[2ex] \displaystyle \hat{I}_h^j(0, t)=0,\\[2ex] \displaystyle \hat{I}_h^j(a, 0)=\psi_j(a), \end{array}$$ \right. $

(52)

and

$ \left\{ $$\begin{array}{ll} \displaystyle \frac{\partial \tilde{E}_h^j}{\partial t}+\frac{\partial \tilde{E}_h^j}{\partial \tau}=-(\mu_h+m_h^j(\tau))\tilde{E}_h^j(\tau, t),\\[2ex] \displaystyle \tilde{E}_h^j(0, t)=\beta_h^jS_hI_v^j,\\ \displaystyle \tilde{E}_h^j(\tau, 0)=0, \end{array}$$ \right. $

(53)

$ \left\{ $$\begin{array}{ll} \displaystyle \frac{\partial \tilde{I}_h^j}{\partial t}+\frac{\partial \tilde{I}_h^j}{\partial a}=-(\mu_h+\alpha_h^j(a)+r_h^j(a))\tilde{I}_h^j(a, t),\\[2ex] \displaystyle \tilde{I}_h^j(0, t)=\int^\infty_0m_h^j(\tau)\tilde{E}_h^j(\tau, t)d\tau,\\[2ex] \displaystyle \tilde{I}_h^j(a, 0)=0. \end{array}$$ \right. $

(54)

We can easily see that system (51) and (52) are decoupled from the remaining equations. Using the formula (7) to integrate along the characteristic lines, we obtain

$ $$\begin{array}{ll} \displaystyle \hat{E}_h^j(\tau, t)=\left\{\begin{array}{ll} \displaystyle 0,\mbox{ } & t>\tau,\\[2ex] \displaystyle \varphi_j(\tau-t)\frac{\pi^j_1(\tau)}{\pi^j_1(\tau-t)},\mbox{ } & t<\tau, \end{array}$$ \right. \end{array} $

(55)

$ $$\begin{array}{ll} \displaystyle \hat{I}_h^j(a, t)=\left\{\begin{array}{ll} \displaystyle 0,\mbox{ } & t>a,\\[2ex] \displaystyle \psi_j(a-t)\frac{\pi^j_2(a)}{\pi^j_2(a-t)},\mbox{ } & t<a. \end{array}$$ \right.\\[2ex] \end{array} $

(56)

Integrating $\hat{E}_h^j$ with respect to $\tau$ yields

$ \displaystyle \int^\infty_t \varphi_j(\tau-t)\frac{\pi^j_1(\tau)}{\pi^j_1(\tau-t)}d\tau=\int^\infty_0 \varphi_j(\tau)\frac{\pi^j_1(t+\tau)}{\pi^j_1(\tau)}d\tau\leq e^{-\mu_h t}\int^\infty_0 \varphi_j(\tau)d\tau\rightarrow 0 $

as $t\rightarrow\infty$. Integrating $\hat{I}_h^j$ with respect to $a$, we have

$ \displaystyle \int^\infty_t \psi_j(a-t)\frac{\pi^j_2(a)}{\pi^j_2(a-t)}da=\int^\infty_0 \psi_j(a)\frac{\pi^j_2(t+a)}{\pi^j_2(a)}da\leq e^{-\mu_h t}\int^\infty_0 \psi_j(a)da \rightarrow 0 $

as $t\rightarrow\infty$. This implies that $\hat{\Psi}(t, x^0)\rightarrow0$ as $t\rightarrow\infty$ uniformly for every $x^0\in\mathcal{B}\subseteq \Omega_0$, where $\mathcal{B}$ is a ball of a given radius.

In the following we need to show $\tilde{\Psi}(t, x)$ is completely continuous. We fix $t$ and let $x^0\in \Omega_0$. Note that $\Omega_0$ is bounded. We have to show that the family of functions defined by

$ $$\begin{array}{ll} \displaystyle \tilde{\Psi}(t, x^0)=\displaystyle \bigg( S_v(t), I_v^1(t),\cdots,I_v^n(t),S_h(t), \tilde{E}_h^1(\tau, t),\tilde{I}_h^1(a,t), \cdots, \tilde{E}_h^n(\tau, t),\tilde{I}_h^n(a,t)\bigg) \end{array}$$ $

is a compact family of functions for that fixed $t$, which are obtained by taking different initial conditions in $\Omega_0$. The family

$ \{\tilde{\Psi}(t, x^0)|x^0\in \Omega_0, t-\text{fixed}\}\subseteq \Omega_0, $

and, therefore, it is bounded. Thus, we have established the boundedness of the set. To show that $\tilde{\Psi}(t, x)$ is precompact, we first see the third condition of $\lim_{t\rightarrow\infty}\int^\infty_t\tilde{E}_h^j(\tau, t)d\tau=0$ and $\lim_{t\rightarrow\infty}\int^\infty_t\tilde{I}_h^j(a, t)da=0$ in the Frechet-Kolmogorov Theorem of [21]. The third condition in [21] is trivially satisfied since $ \tilde{E}_h^j(\tau, t)=0$ for $ \tau>t$ and $ \tilde{I}_h^j(a, t)=0$ for $ a>t$. To use the second condition of the Frechet-Kolmogorov Theorem in [21], we must bound by two constants the L$^1$-norms of $\partial E_h^j/\partial\tau$ and $\partial I_h^j/\partial a$. Notice that

$ $$\begin{array}{ll} \displaystyle \tilde{E}_h^j(\tau, t)=\left\{\begin{array}{ll} \displaystyle \tilde{B}_E^j(t-\tau)\pi^j_1(\tau),\mbox{ } & t>\tau,\\[2ex] \displaystyle 0,\mbox{ } & t<\tau,\\[2ex] \end{array}$$ \right.\\[2ex] \displaystyle \tilde{I}_h^j(a, t)=\left\{\begin{array}{ll} \displaystyle \tilde{B}_I^j(t-a)\pi^j_2(a),\mbox{ } & t>a,\\[2ex] \displaystyle 0,\mbox{ } & t<a, \end{array}\right. \end{array} $

(57)

where

$ $$\begin{array}{ll} \displaystyle \tilde{B}_E^j(t)=\displaystyle \beta_h^jS_h(t)I_v^j(t),\\[2ex] \displaystyle \tilde{B}_I^j(t)=\displaystyle \int^\infty_0m_h^j(\tau)\tilde{E}_h^j(\tau, t)d\tau=\int^t_0m_h^j(\tau)\tilde{B}_E^j(t-\tau)\pi^j_1(\tau)d\tau.\\[2ex] \end{array}$$ $

(58)

$\tilde{B}_E^j(t)$ is bounded because of the boundedness of $S_h$ and $I_v^j$. Hence, the $\tilde{B}_E^j(t)$ satisfies

$ $$\begin{array}{ll} \displaystyle \tilde{B}_E^j(t)\leq k_1. \end{array}$$ $

Therefore, we obtain

$ $$\begin{array}{ll} \displaystyle \tilde{B}_I^j(t) & \displaystyle=\int^t_0m_h^j(\tau)\tilde{B}_E^j(t-\tau)\pi^j_1(\tau)d\tau\\[2ex] & \displaystyle\leq k_2\int^t_0\tilde{B}_E^j(t-\tau)d\tau= k_2\int^t_0\tilde{B}_E^j(\tau)d\tau\\[2ex] & \displaystyle \leq k_1k_2t. \end{array}$$ $

Next, we differentiate (57) with respect to $\tau $ and $a$:

$ $$\begin{array}{ll} \displaystyle \bigg|\frac{\partial\tilde{E}_h^j(\tau, t)}{\partial\tau}\bigg|\leq\left\{\begin{array}{ll} |(\tilde{B}_E^j(t-\tau))'|\pi^j_1(\tau)+\tilde{B}_E^j(t-\tau)|(\pi^j_1(\tau))'|,\mbox{} & t>\tau,\\[2ex] \displaystyle 0,\mbox{ } & t<\tau, \end{array}$$ \right.\\[2ex] \displaystyle \bigg|\frac{\partial\tilde{I}_h^j(a, t)}{\partial a}\bigg|\leq\left\{\begin{array}{ll} |(\tilde{B}_I^j(t-a))'|\pi^j_2(a)+\tilde{B}_I^j(t-a)|(\pi^j_2(a))'|,\mbox{} & t>a,\\[2ex] \displaystyle 0,\mbox{ } & t<a. \end{array}\right.\\[2ex] \end{array} $

We see that $|(\tilde{B}_E^j(t-\tau))'|,\ |(\tilde{B}_I^j(t-a))'|$ are bounded. Differentiating (58), we obtain

$ $$\begin{array}{ll} \displaystyle (\tilde{B}_E^j(t))'=\displaystyle \beta_h^j\bigg(S'_h(t)I_v^j(t)+S_h(t)(I_v^j(t))'\bigg),\\[2ex] (\tilde{B}_I^j(t))'=\displaystyle m_h^j(t)\tilde{B}_E^j(0)\pi^j_1(t)+\int^t_0m_h^j(\tau)(\tilde{B}_E^j(t-\tau))'\pi^j_1(\tau)d\tau. \end{array}$$ $

(59)

Taking an absolute value and bounding all terms, we can rewrite the above equality as the following inequality:

$ $$\begin{array}{ll} \displaystyle |(\tilde{B}_E^j(t))'|\leq k_3, |(\tilde{B}_I^j(t))'|\leq k_4. \end{array}$$ $

Putting all these bounds together, we have

$ $$\begin{array}{ll} \displaystyle \parallel\partial_\tau\tilde{E}_h^j\parallel & \displaystyle\leq k_3\int^\infty_0\pi^j_1(\tau)d\tau+k_1(\mu_h+\bar{m}_h)\int^\infty_0\pi^j_1(\tau)d\tau<\displaystyle \mathfrak{b}_1,\\[2ex] \displaystyle \parallel\partial_a\tilde{I}_h^j\parallel & \displaystyle\leq k_4\int^\infty_0\pi^j_2(a)da+k_1k_2(\mu_h+\bar{\alpha}_h+\bar{r}_h)t\int^\infty_0\pi^j_2(a)da<\displaystyle \mathfrak{b}_2, \end{array}$$ $

where $\bar{m}_h=\sup_{\tau,j}\{m_h^j(\tau)\},\ \bar{\alpha}_h=\sup_{a,j}\{\alpha_h^j(a)\},\ \bar{r}_h=\sup_{a,j}\{r_h^j(a)\}$. To complete the proof, we notice that

$ $$\begin{array}{ll} \displaystyle \int^\infty_0|\tilde{E}_h^j(\tau+h, t)-\tilde{E}_h^j(\tau, t)|d\tau\leq\parallel\partial_\tau\tilde{E}_h^j\parallel |h|\leq \mathfrak{b}_1|h|,\\[2ex] \displaystyle \int^\infty_0|\tilde{I}_h^j(a+h, t)-\tilde{I}_h^j(a, t)|da\leq\parallel\partial_a\tilde{I}_h^j\parallel |h|\leq \mathfrak{b}_2|h|.\\[2ex] \end{array}$$ $

Thus, the integral can be made arbitrary small uniformly in the family of functions. That establishes the second condition of the Frechet-Kolmogorov Theorem. We conclude that the family is asymptotically smooth.

(3) means that the semigroup $\Psi(t)$ is point dissipative and the forward orbit of boundedness sets is bounded in $\Omega_0$. Thus, we prove Theorem 5.5 in accordance with Lemma 3.1.3 and Theorem 3.4.6 in [8].

Now we have all components to establish the uniform strong persistence. The next proposition states the uniform strong persistence of $I_v^1,E_h^1$ and $I_h^1$.

Theorem 5.6. Under the hypothesis of Theorem 5.3 strain one is uniformly strongly persistent for all initial conditions that belong to $\Omega_0$, that is, there exists $\gamma>0$ such that

$ \displaystyle \liminf\limits_{t\rightarrow+\infty}\beta_h^1I_v^1(t)\geq\gamma, \liminf\limits_{t\rightarrow+\infty}\int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau\geq\gamma, \liminf\limits_{t\rightarrow+\infty}\int^\infty_0 \beta_v^1(a)I_h^1(a, t)da\geq\gamma. $ }

Proof. We apply Theorem 2.6 in [19]. We consider the solution semiflow $\Psi$ on $\Omega_0$. Let us define three functionals $\rho_l: \Omega_0\rightarrow$R$_+,\ l=1,2,3$ as follows:

$ \left. $$\begin{array}{ll} \displaystyle \rho_1(\Psi(t, x^0))=\beta_h^1I_v^1(t),\\[2ex] \displaystyle \rho_2(\Psi(t, x^0))=\int^\infty_0 m_h^1(\tau)\tilde{E}_h^1(\tau, t)d\tau,\\[2ex] \displaystyle \rho_3(\Psi(t, x^0))=\int^\infty_0 \beta_v^1(a)\tilde{I}_h^1(a, t)da. \end{array}$$ \right. $

Theorem 5.3 implies that the semiflow is uniformly weakly $\rho$-persistent. Theorem 5.5 shows that the solution semiflow has a global compact attractor $\mathfrak{T}$. Total orbits are solutions to the system (2) defined for all times $t\in\mathbb{R}$. Since the solution semiflow is nonnegative, we have

$\displaystyle \beta_h^1I_v^1(t)= \beta_h^1I_v^1(s)e^{-(\mu_v+\alpha_v^1)(t-s)},$

$ $$\begin{array}{rl} \displaystyle \int^\infty_0 m_h^1(\tau)\tilde{E}_h^1(\tau, t)d\tau & \displaystyle =\tilde{B}_I^1(t)=\int^t_0m_h^1(\tau)\tilde{B}_E^1(t-\tau)\pi^1_1(\tau)d\tau \\[2ex] & \displaystyle\geq k^1\int^t_0\tilde{B}_E^1(t-\tau)d\tau= k^1\int^t_0\tilde{B}_E^1(\tau)d\tau\\[2ex] & \displaystyle=k^1\int^t_0\beta_h^1S_h(\tau)I_v^1(\tau)d\tau\geq k^2\int^t_0I_v^1(\tau)d\tau \end{array}$$ $

$ $$\begin{array}{rl} & \displaystyle=k^2\int^t_0I_v^1(s)e^{-(\mu_v+\alpha_v^1)(\tau-s)}d\tau\\[2ex] & \displaystyle=\frac{k^2I_v^1(s)}{\mu_v+\alpha_v^1}e^{(\mu_v+\alpha_v^1)s}(1-e^{-(\mu_v+\alpha_v^1)t}),\\[2ex] \displaystyle \int^\infty_0 \beta_v^1(a)\tilde{I}_h^1(a, t)da & \displaystyle=\int^t_0 \beta_v^1(a)\tilde{B}_I^1(t-a)\pi^1_2(a) da\geq k^3\int^t_0\tilde{B}_I^1(t-a)da\\[2ex] & \displaystyle= k^3\int^t_0\tilde{B}_I^1(a)da\\[2ex] & \displaystyle\geq\frac{k^2k^3I_v^1(s)}{\mu_v+\alpha_v^1}e^{(\mu_v+\alpha_v^1)s}\int^t_0(1-e^{-(\mu_v+\alpha_v^1)a})da, \end{array}$$ $

for any $s$ and any $t>s$. Therefore,

$\beta_h^1I_v^1(t)>0, \int^\infty_0 m_h^1(\tau)\tilde{E}_h^1(\tau, t)d\tau>0, \int^\infty_0\beta_v^1(a)\tilde{I}_h^1(a, t)da>0$

for all $t>s$ provided $\tilde{I}_v^1(s)>0$. Theorem 2.6 in [19] now implies that the semiflow is uniformly strongly $\rho$-persistent. Hence, there exists $\gamma$ such that

$ \displaystyle \liminf\limits_{t \rightarrow +\infty}\beta_h^1I_v^1(t)\geq\gamma, \liminf\limits_{t \rightarrow +\infty}\int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau\geq\gamma, \liminf\limits_{t \rightarrow +\infty}\int^\infty_0 \beta_v^1(a)I_h^1(a, t)da\geq\gamma. $

According to Theorem 5.6, we obtain that for all initial conditions that belong to $\Omega_0$, strain 1 persists. Furthermore we had verified that the solutions of (2) with nonnegative initial conditions belong to the positive cone for all $t\geq0$. All the solutions are in a positively invariant set. Therefore we can obtain the following Theorem 5.7 from Theorem 5.6.

Theorem 5.7. Under the hypothesis of Theorem 5.3, $\forall t\in\mathrm{R}$, there exists constants $\vartheta>0$ and $M>0$ such that

$ \vartheta\leq S_v(t)\leq M, \vartheta\leq S_h(t)\leq M,$

and

$ \displaystyle\vartheta\leq \beta_h^1I_v^1(t)\leq M,\ \displaystyle\vartheta\leq \int^\infty_0 m_h^1(\tau)E_h^1(\tau, t)d\tau\leq M,\ \displaystyle\vartheta\leq \int^\infty_0 \beta_v^1(a)I_h^1(a, t)da\leq M, $

for each orbit $( S_v(t), I_v^1(t),\cdots, I_v^n(t), S_h(t), E_h^1(\tau, t), I_h^1(a,t),\cdots,E_h^n(\tau, t), I_h^n(a,t))$ of $\Psi$ in $\mathfrak{T}$.

---

6. Principle of competitive exclusion

In this section we mainly state the main result of the paper.

Theorem 6.1. Assume $\mathcal R^1_0>1,\mathcal R^i_0/\mathcal R^1_0<b_i/b_1<1,i=2,\cdots,n$. Then the equilibrium $\mathcal{E}_1$ is globally asymptotically stable.

Proof. From Theorem 4.1 we know that the endemic equilibrium $\mathcal{E}_1$ is locally asymptotically stable. In the following we only need to show that the endemic equilibrium $\mathcal{E}_1$ is global attractor. From Theorem 5.5 there exists an invariant compact set $\mathfrak{T}$ which is global attractor of system (2). Furthermore, it follows from Theorem 5.7 that there exist $\varepsilon_1>0$ and $M_1>0$ such that

$ \displaystyle \varepsilon_1\leq\frac{I_v^1}{I_v^{1^*}}\leq M_1, \varepsilon_1\leq \frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\leq M_1, \varepsilon_1\leq \frac{I_h^1(a, t)}{I_h^{1^*}(a)}\leq M_1 $

for any solution in $\Psi$. This makes the Lyapunov function defined in (40) well defined.

After extensive computation, differentiating $U(t)$ along the solution of system (2) yields that

$ $$\begin{split} & \displaystyle \frac{dU_1(t)}{dt}\\ = & \displaystyle \frac{1}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\bigg(1-\frac{S_v^{1^*}}{S_v }\bigg) \bigg[S_v^{1^*}\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da+\mu_vS_v^{1^*}\\ & \displaystyle-S_v\int^\infty_0\beta_v^1(a)I_h^1(a, t)da-\mu_vS_v-\sum^n_{i=2}S_v\int^\infty_0\beta_v^i(a)I_h^i(a, t)da\bigg]\\ = & \displaystyle -\frac{\mu_v(S_v-S_v^{1^*})^2}{S_v^{1^*}S_vq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\ & \displaystyle+\frac{1}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\int^\infty_0 \beta_v^1(a)I_h^{1^*}(a)\bigg(1-\frac{S_v^{1^*}}{S_v}-\frac{S_vI_h^1(a, t)} {S_v^{1^*}I_h^{1^*}(a)}+\frac{I_h^1(a, t)}{I_h^{1^*}(a)}\bigg)da\\ & \displaystyle-\sum^n_{i=2} \frac{S_v\int^\infty_0\beta_v^i(a)I_h^i(a, t)da-S_v^{1^*}\int^\infty_0\beta_v^i(a)I_h^i(a, t)da}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}; \end{split}$$ $

(60)

$ $$\begin{array}{ll} & \displaystyle \frac{dU^1_2(t)}{dt}\\[3ex] = & \displaystyle \frac{(1-\frac{I_v^{1^*}}{I_v^1})\bigg(S_v\int^\infty_0\beta_v^1(a)I_h^1(a, t)da-\frac{S_v^{1^*}\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da}{I_v^{1^*}}I_v^1\bigg)}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[3ex] = & \displaystyle \frac{(1-\frac{I_v^{1^*}}{I_v^1})S_v^{1^*}\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)\bigg(\frac{S_vI_h^1(a, t)} {S_v^{1^*}I_h^{1^*}(a)}-\frac{I_v^1}{I_v^{1^*}}\bigg)da}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[3ex] = & \displaystyle \frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)\bigg(\frac{S_vI_h^1(a, t)} {S_v^{1^*}I_h^{1^*}(a)}-\frac{I_v^1}{I_v^{1^*}}-\frac{S_vI_h^1(a, t)I_v^{1^*}} {S_v^{1^*}I_h^{1^*}(a)I_v^1}+1\bigg)da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}; \end{array}$$ $

(61)

$ $$\begin{array}{ll} \displaystyle \frac{dU^i_2(t)}{dt}=\frac{S_v\int^\infty_0\beta_v^i(a)I_h^i(a, t)da-(\mu_v+\alpha_v^i)I_v^i}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}; \end{array}$$ $

(62)

$ $$\begin{array}{ll} \displaystyle \frac{dU_3(t) }{dt} & =\displaystyle \bigg(1-\frac{S_h^{1^*}}{S_h}\bigg)\bigg(E_h^{1^*}(0)+\mu_hS_h^{1^*}-E_h^1(0,t)-\mu_hS_h-\sum^n_{i=2}\beta_h^iS_hI_v^i\bigg)\\[2ex] & =\displaystyle -\frac{\mu_h(S_h-S_h^{1^*})^2}{S_h} \displaystyle+\bigg(E_h^{1^*}(0)-E_h^1(0,t)-\frac{S_h^{1^*}}{S_h}E_h^{1^*}(0)+\frac{S_h^{1^*}}{S_h}E_h^1(0,t)\bigg)\\[2ex] & \displaystyle-\sum^n_{i=2}\bigg(E_h^i(0,t)-\beta_h^iS_h^{1^*}I_v^i\bigg),\\[2ex] \end{array}$$ $

(63)

and

$ $$\begin{array}{ll} & \displaystyle \frac{dU^1_4(t)}{dt}\\[2ex] = & \displaystyle \frac{1}{\mathcal R^1_0}\int^\infty_0p_1(\tau)E_h^{1^*}(\tau)f'\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)\frac{1}{E_h^{1^*}(\tau)}\frac{\partial E_h^1(\tau, t)}{\partial t}d\tau\\[2ex] = & \displaystyle -\frac{1}{\mathcal R^1_0}\int^\infty_0p_1(\tau)E_h^{1^*}(\tau)f'\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)\frac{1}{E_h^{1^*}(\tau)}\bigg(\frac{\partial E_h^1(\tau, t)}{\partial \tau}+(\mu_h+m_h^1(\tau))E_h^1(\tau, t)\bigg)d\tau\\[2ex] = & \displaystyle -\frac{1}{\mathcal R^1_0}\int^\infty_0p_1(\tau)E_h^{1^*}(\tau)df\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)\\[2ex] = & \displaystyle -\frac{1}{\mathcal R^1_0}\bigg[p_1(\tau)E_h^{1^*}(\tau)f\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)\bigg|^{\infty}_0-\int^\infty_0f\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)d\bigg(p_1(\tau)E_h^{1^*}(\tau)\bigg)\bigg]\\[2ex] = & \displaystyle \frac{1}{\mathcal R^1_0}\bigg[p_1(0)E_h^{1^*}(0)f\bigg(\frac{E_h^1(0, t)}{E_h^{1^*}(0)}\bigg)-\Delta_1 q_1(0)\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)d\tau\bigg]\\[2ex] = & \displaystyle E_h^{1^*}(0)f\bigg(\frac{E_h^1(0, t)}{E_h^{1^*}(0)}\bigg)-\frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)})d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] = & \displaystyle E_h^1(0, t)-E_h^{1^*}(0)-E_h^{1^*}(0)\ln\frac{E_h^1(0, t)}{E_h^{1^*}(0)}-\frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)})d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}. \end{array}$$ $

(64)

The above equality follows from (24) and the fact

$ $$\begin{array}{ll} & \displaystyle p_1'(\tau)E_h^{1^*}(\tau)+p_1(\tau)(E_h^{1^*}(\tau))'\\[2ex] = & \displaystyle \bigg[-\Delta_1 q_1(0)m_h^1(\tau)+(\mu_h+m_h^1(\tau))p_1(\tau)\bigg]E_h^{1^*}(\tau)-p_1(\tau)(\mu_h+m_h^1(\tau))E_h^{1^*}(\tau) \\[2ex] = & \displaystyle -\Delta_1 q_1(0)m_h^1(\tau)E_h^{1^*}(\tau). \end{array}$$ $

We also have

$ $$\begin{array}{ll} & \displaystyle q_1'(a)I_h^{1^*}(a)+q_1(a)(I_h^{1^*}(a))'\\[2ex] = & \displaystyle \bigg[-\beta_v^1(a)+(\mu_h+\alpha_h^1(a)+r_h^1(a))q_1(a)\bigg]I_h^{1^*}(a)-q_1(a)(\mu_h+\alpha_h^1(a)+r_h^1(a))I_h^{1^*}(a) \\[2ex] = & \displaystyle -\beta_v^1(a)I_h^{1^*}(a). \end{array}$$ $

Similar to the differentiation of $U^1_4(t)$, we have

$ $$\begin{array}{ll} & \displaystyle \frac{dU^1_5(t)}{dt}\\[2ex] = & \displaystyle\frac{1}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\displaystyle \int^\infty_0q_1(a)I_h^{1^*}(a)f'\bigg(\frac{I_h^1(a, t)}{I_h^{1^*}(a)}\bigg)\frac{1}{I_h^{1^*}(a)}\frac{\partial I_h^1(a, t)}{\partial t}da\\[2ex] = & \displaystyle \frac{1}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\displaystyle \int^\infty_0q_1(a)I_h^{1^*}(a)df\bigg(\frac{I_h^1(a, t)}{I_h^{1^*}(a)}\bigg)\\[2ex] = & \displaystyle \frac{q_1(0)I_h^{1^*}(0)f(\frac{I_h^1(0, t)}{I_h^{1^*}(0)})-\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)f(\frac{I_h^1(a, t)}{I_h^{1^*}(a)})da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] = & \displaystyle \frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)(\frac{I_h^1(0, t)}{I_h^{1^*}(0)}-1-\ln\frac{I_h^1(0, t)}{I_h^{1^*}(0)})d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}-\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)f(\frac{I_h^1(a, t)}{I_h^{1^*}(a)})da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}.\\[2ex] \end{array}$$ $

(65)

Noting that (43), we differentiate the last two terms with respect to $t$, and have

$ $$\begin{array}{ll} & \displaystyle \frac{dU^i_4(t)}{dt}\\[2ex] = & \displaystyle \frac{1}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\int^\infty_0p_i(\tau)\frac{\partial E_h^i(\tau,t)}{\partial t}d\tau \\[2ex] = & \displaystyle \displaystyle -\frac{1}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\int^\infty_0p_i(\tau)\bigg[\frac{\partial E_h^i(\tau,t)}{\partial \tau}+(\mu_h+m_h^i(\tau))E_h^i(\tau,t)\bigg]d\tau \\[2ex] = & \displaystyle \displaystyle -\frac{\int^\infty_0p_i(\tau)d E_h^i(\tau,t) +\int^\infty_0(\mu_h+m_h^i(\tau))p_i(\tau)E_h^i(\tau,t)d\tau}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] = & \displaystyle -\frac{p_i(\tau) E_h^i(\tau,t)|^\infty_0 -\int^\infty_0E_h^i(\tau,t)d p_i(\tau) +\int^\infty_0(\mu_h+m_h^i(\tau))p_i(\tau)E_h^i(\tau,t)d\tau}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] = & \displaystyle \frac{p_i(0) E_h^i(0,t) -\Delta_i q_i(0)\int^\infty_0m_h^i(\tau)E_h^i(\tau,t)d\tau}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] = & \displaystyle \frac{\mathcal R^i_0 E_h^i(0,t)}{\Delta_iq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} -\frac{ q_i(0)I_h^i(0,t)}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] = & \displaystyle\frac{b_i}{b_1}E_h^i(0,t) -\frac{ q_i(0)I_h^i(0,t)}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}. \end{array}$$ $

(66)

Similarly, we have

$ $$\begin{array}{ll} & \displaystyle \frac{dU^i_5(t)}{dt}\\[2ex] = & -\displaystyle \frac{\int^\infty_0q_i(a)\bigg[\frac{\partial I_h^i(a,t)}{\partial a}+(\mu_h+\alpha_h^i(a)+r_h^i(a))I_h^i(a,t)\bigg]da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] = & \displaystyle \frac{q_i(0) I_h^i(0,t)}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} - \frac{\int^\infty_0\beta_v^i(a)I_h^i(a,t)da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}.\\[2ex] \end{array}$$ $

(67)

Adding all five components of the Lyapunov function, we have

$ \displaystyle U'(t)=U^1+U^2, $

where

$ $$\begin{array}{ll} & \displaystyle U^1(t)\\[2ex] = & \displaystyle -\frac{\mu_v(S_v-S_v^{1^*})^2}{S_v^{1^*}S_vq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] & \displaystyle+\frac{\int^\infty_0 \beta_v^1(a)I_h^{1^*}(a)\bigg(1-\frac{S_v^{1^*}}{S_v}-\frac{S_vI_h^1(a, t)} {S_v^{1^*}I_h^{1^*}(a)}+\frac{I_h^1(a, t)}{I_h^{1^*}(a)}\bigg)da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] & +\displaystyle \frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)\bigg(\frac{S_vI_h^1(a, t)} {S_v^{1^*}I_h^{1^*}(a)}-\frac{I_v^1}{I_v^{1^*}}-\frac{S_vI_h^1(a, t)I_v^{1^*}} {S_v^{1^*}I_h^{1^*}(a)I_v^1}+1\bigg)da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] & \displaystyle -\frac{\mu_h(S_h-S_h^{1^*})^2}{S_h} \displaystyle+\bigg(E_h^{1^*}(0)-E_h^1(0,t)-\frac{S_h^{1^*}}{S_h}E_h^{1^*}(0)+\frac{S_h^{1^*}}{S_h}E_h^1(0,t)\bigg) \end{array}$$ $

$ $$\begin{array}{ll} & \displaystyle+ E_h^1(0, t)-E_h^{1^*}(0)-E_h^{1^*}(0)\ln\frac{E_h^1(0, t)}{E_h^{1^*}(0)}-\frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)})d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] & \displaystyle +\frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)(\frac{I_h^1(0, t)}{I_h^{1^*}(0)}-1-\ln\frac{I_h^1(0, t)}{I_h^{1^*}(0)})d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}-\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)f(\frac{I_h^1(a, t)}{I_h^{1^*}(a)})da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}, \end{array}$$ $

(68)

and

$ $$\begin{array}{ll} & \displaystyle U^2(t)\\[2ex] = & \displaystyle-\sum^n_{i=2} \frac{S_v\int^\infty_0\beta_v^i(a)I_h^i(a, t)da-S_v^{1^*} \int^\infty_0\beta_v^i(a)I_h^i(a, t)da}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] & +\displaystyle\sum^n_{i=2}\frac{S_v\int^\infty_0\beta_v^i(a)I_h^i(a, t)da-(\mu_v+\alpha_v^i)I_v^i}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}-\sum^n_{i=2}\bigg(E_h^i(0,t)-\beta_h^iS_h^{1^*}I_v^i\bigg)\\[2ex] & +\displaystyle\sum^n_{i=2} \bigg(\frac{b_i}{b_1}E_h^i(0,t) -\frac{q_i(0)I_h^i(0,t)}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\bigg)\\[2ex] & + \displaystyle \sum^n_{i=2}\bigg(\frac{q_i(0) I_h^i(0,t)}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} - \frac{\int^\infty_0\beta_v^i(a)I_h^i(a,t)da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\bigg).\\[2ex] \end{array}$$ $

(69)

Canceling terms, (68) can be simplified as

$ $$\begin{array}{ll} & \displaystyle U^1(t) \\[2ex] = & \displaystyle -\frac{\mu_v(S_v-S_v^{1^*})^2}{S_v^{1^*}S_vq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} -\frac{\mu_h(S_h-S_h^{1^*})^2}{S_h}\\[2ex] & \displaystyle+ \frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)(3-\frac{S_v^{1^*}}{S_v}-\frac{I_v^1}{I_v^{1^*}}-\frac{S_vI_h^1(a, t)I_v^{1^*}} {S_v^{1^*}I_h^{1^*}(a)I_v^1}+\ln\frac{I_h^1(a, t)} {I_h^{1^*}(a)})da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\\[2ex] & \displaystyle +E_h^{1^*}(0)\bigg(-\frac{S_h^{1^*}}{S_h}+\frac{S_h^{1^*}E_h^1(0, t)}{S_hE_h^{1^*}(0)}-\ln\frac{E_h^1(0, t)}{E_h^{1^*}(0)} \bigg)\\[2ex] & \displaystyle+ \frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)(\frac{I_h^1(0, t)}{I_h^{1^*}(0)}-\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}+\ln\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\frac{I_h^{1^*}(0)}{I_h^1(0, t)})d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}.\\[2ex] \end{array}$$ $

(70)

Direct computation yields that

$ $$\begin{array}{ll} & \displaystyle \int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)\bigg(\frac{I_h^1(0, t)}{I_h^{1^*}(0)}-\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\bigg)d\tau\\[2ex] = & \displaystyle \frac{I_h^1(0, t)}{I_h^{1^*}(0)}\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)d\tau-\int^\infty_0m_h^1(\tau)E_h^1(\tau, t)d\tau\\[2ex] = & \displaystyle \frac{I_h^1(0, t)}{I_h^{1^*}(0)}I_h^{1^*}(0)-I_h^1(0, t)=0,\\[2ex] & \displaystyle \int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)\bigg(\frac{E_h^1(\tau, t)}{E_h^{1^*}(\tau)}\frac{I_h^{1^*}(0)}{I_h^1 (0, t)}-1\bigg)\\[2ex] = & \displaystyle \frac{I_h^{1^*}(0)}{I_h^1 (0, t)}\int^\infty_0m_h^1(\tau)E_h^1(\tau, t)d\tau-\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)d\tau\\[2ex] = & \displaystyle \frac{I_h^{1^*}(0)}{I_h^1 (0, t)}I_h^1 (0, t)-I_h^{1^*}(0)=0. \end{array}$$ $

(71)

By using (71), (70) can be simplified as

$ $$\begin{array}{ll} \displaystyle U^1(t) & =\displaystyle -\frac{\mu_v(S_v-S_v^{1^*})^2}{S_v^{1^*}S_vq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} -\frac{\mu_h(S_h-S_h^{1^*})^2}{S_h}\\[2ex] & \displaystyle -\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)[f(\frac{S_v^{1^*}}{S_v}) +f(\frac{I_v^1}{I_v^{1^*}})+f(\frac{S_vI_h^1(a, t)I_v^{1^*}} {S_v^{1^*}I_h^{1^*}(a)I_v^1})]da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] & \displaystyle -\frac{1}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f\bigg(\frac{E_h^1(\tau, t)I_h^{1^*}(0)}{E_h^{1^*}(\tau)I_h^1(0, t)}\bigg)d\tau\\[2ex] & \displaystyle +E_h^{1^*}(0)\bigg[-f\bigg(\frac{S_h^{1^*}}{S_h}\bigg)+f\bigg(\frac{S_h^{1^*}E_h^1(0, t)}{S_hE_h^{1^*}(0)}\bigg)\bigg].\\[2ex] \end{array}$$ $

(72)

Noting that $E_h^1(0, t)=\beta_h^1 S_h I_v^1,\ E_h^{1^*}(0)=\beta_h^1 S_h^{1^*} I_v^{1^*}$, we get

$ \displaystyle \frac{S_h^{1^*}E_h^1(0,t)}{S_hE_h^{1^*}(0)}=\frac{S_h^{1^*}\beta_h^1 S_h I_v^1}{S_h\beta_h^1 S_h^{1^*} I_v^{1^*}}=\frac{I_v^1}{I_v^{1^*}}. $

(73)

Furthermore, from (25) and (42) we have

$ $$\begin{array}{rl} \displaystyle\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)f(\frac{I_v^1}{I_v^{1^*}})da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} & =\displaystyle\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}f(\frac{I_v^1}{I_v^{1^*}})\\[2ex] & =\displaystyle \frac{I_h^{1^*}(0)}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}f(\frac{I_v^1}{I_v^{1^*}})\\[2ex] & =\displaystyle E_h^{1^*}(0)f(\frac{I_v^1}{I_v^{1^*}}). \end{array}$$ $

(74)

Finally, simplifying (72) with (73) and (74), we obtain

$ $$\begin{array}{ll} \displaystyle U^1(t) & =\displaystyle -\frac{\mu_v(S_v-S_v^{1^*})^2}{S_v^{1^*}S_vq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} -\frac{\mu_h(S_h-S_h^{1^*})^2}{S_h}\\[2ex] & \displaystyle -\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)[f(\frac{S_v^{1^*}}{S_v}) +f(\frac{S_vI_h^1(a, t)I_v^{1^*}} {S_v^{1^*}I_h^{1^*}(a)I_v^1})]da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] & \displaystyle -\frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f\bigg(\frac{E_h^1(\tau, t)I_h^{1^*}(0)}{E_h^{1^*}(\tau)I_h^1(0, t)}\bigg)d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}-E_h^{1^*}(0)f\bigg(\frac{S_h^{1^*}}{S_h}\bigg) . \end{array}$$ $

(75)

Canceling terms, (69) can be simplified as

$ $$\begin{array}{ll} & \displaystyle U^2(t)\\[2ex] = & \displaystyle\sum^n_{i=2}\bigg[ \bigg(\frac{b_i}{b_1}-1\bigg)E_h^i(0,t)+ \bigg(\beta_h^iS_h^{1^*}-\frac{\mu_v+\alpha_v^i}{S_v^{1^*}q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}\bigg)I_v^i\bigg]. \end{array}$$ $

(76)

Simplifying (76) with (25), we get

$ $$\begin{array}{ll} \displaystyle U^2(t)= & \displaystyle \sum^n_{i=2}\bigg[ \bigg(\frac{b_i}{b_1}-1\bigg)E_h^i(0,t) +\frac{\beta_h^i\Lambda_h(\mu_v+\Lambda_hb_1)}{\mu_h(\mu_v\mathcal R^1_0+\Lambda_hb_1)} \bigg(1-\frac{\mathcal R^1_0b_i}{\mathcal R^i_0b_1}\bigg)I_v^i\bigg]. \end{array}$$ $

(77)

Hence, by using (75) and (77) we obtain

$ $$\begin{array}{ll} \displaystyle U'(t) & =\displaystyle -\frac{\mu_v(S_v-S_v^{1^*})^2}{S_v^{1^*}S_vq_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} -\frac{\mu_h(S_h-S_h^{1^*})^2}{S_h}\\[2ex] & \displaystyle -\frac{\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)[f(\frac{S_v^{1^*}}{S_v}) +f(\frac{S_vI_h^1(a, t)I_v^{1^*}} {S_v^{1^*}I_h^{1^*}(a)I_v^1})]da}{q_1(0)\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau} \\[2ex] & \displaystyle -\frac{\int^\infty_0m_h^1(\tau)E_h^{1^*}(\tau)f\bigg(\frac{E_h^1(\tau, t)I_h^{1^*}(0)}{E_h^{1^*}(\tau)I_h^1(0, t)}\bigg)d\tau}{\int^\infty_0m_h^1(\tau)\pi^1_1(\tau)d\tau}-E_h^{1^*}(0)f\bigg(\frac{S_h^{1^*}}{S_h}\bigg) \\[2ex] & \displaystyle+\sum^n_{i=2}\bigg[ \bigg(\frac{b_i}{b_1}-1\bigg)E_h^i(0,t) +\frac{\beta_h^i\Lambda_h(\mu_v+\Lambda_hb_1)}{\mu_h(\mu_v\mathcal R^1_0+\Lambda_hb_1)} \bigg(1-\frac{\mathcal R^1_0b_i}{\mathcal R^i_0b_1}\bigg)I_v^i\bigg]. \end{array}$$ $

(78)

Since $f(x)\geq0$ for $x>0$, $\mathcal R^i_0/\mathcal R^1_0<b_i/b_1<1,i\neq1$ we have $U'\leq0$. Define,

$ \displaystyle \Theta_2= \bigg\{(S_v, I_v^1,\cdots,I_v^n, S_h, E_h^1, I_h^1,\cdots,E_h^n, I_h^n)\in \Omega_0\bigg|U'(t)=0\bigg\}. $

We want to show that the largest invariant set in $\Theta_2$ is the singleton $\mathcal{E}_1$. First, we notice that equality in (78) occurs if and only if $S_v( t)=S_v^{1^*},\ S_h( t)=S_h^{1^*},\ E_h^i(0,t)=0,\ I_v^i=0$, and

$ \displaystyle \frac{I_h^1(a, t)I_v^{1^*}} {I_h^{1^*}(a)I_v^1}=1, \frac{E_h^1(\tau, t)I_h^{1^*}(0)}{E_h^{1^*}(\tau)I_h^1(0, t)}=1. $

(79)

Thus, we obtain

$ \displaystyle \frac{I_h^1(a, t)} {I_h^{1^*}(a)}=\frac{I_v^1(t)} {I_v^{1^*}}. $

(80)

It is obvious that the left term $\frac{I_h^1(a, t)} {I_h^{1^*}(a)}$ of (80) is a function with $a, t$, while the right term $\frac{I_v^1(t)} {I_v^{1^*}}$ is a function with $t$. So we can assume that $I_h^1(a, t)=I_h^{1^*}(a)g(t)$. Thus we have

$ \displaystyle I_v^1=I_v^{1^*}g(t). $

(81)

It follows from (2) we can also obtain

$ $$\begin{array}{ll} \displaystyle I_v^{1'}(t) & \displaystyle= S_v\int^\infty_0\beta_v^1(a)I_h^1(a,t)da-(\mu_v+\alpha_v^1)I_v^1,\\[2ex] & \displaystyle= S_v^{1^*}\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)g(t)da-(\mu_v+\alpha_v^1)I_v^1,\\[2ex] & \displaystyle= g(t)S_v^{1^*}\int^\infty_0\beta_v^1(a)I_h^{1^*}(a)da-(\mu_v+\alpha_v^1)I_v^1,\\[2ex] & \displaystyle= g(t)(\mu_v+\alpha_v^1)I_v^{1^*}-(\mu_v+\alpha_v^1)I_v^1,\\[2ex] & \displaystyle=(\mu_v+\alpha_v^1)(I_v^{1^*}g(t)-I_v^1)=0. \end{array}$$ $

(82)

Therefore, we can get

$I_v^1=I_v^{1^*}.$

Subsequently, it follows from (80) we have

$I_h^1(a,t)=I_h^{1^*}(a).$

Specially, when $a=0$, we have $I_h^1(0,t)=I_h^{1^*}(0).$ Thus from (79) we get

$ E_h^1(\tau, t)=E_h^{1^*}(\tau). $

Since $E_h^i(0,t)=0$, then $E_h^i(\tau,t)=E_h^i(0,t-\tau)\pi^i_1(\tau)=0$ for $t>\tau,i=2,\cdots,n$. Similarly, we also have $I_h^i(0,t)=\int^\infty_0m_h^i(\tau)E_h^i(\tau,t)d\tau=0,\ I_h^i(a,t)=I_h^i(0,t-a)\pi^i_2(a)=0$ for $t>a$. At last we conclude that the largest invariant set in $\Theta_2$ is the singleton $\mathcal{E}_1$. Since $\Psi(t)\Omega_0^+\subset \Omega^+_0$, the global attractor, $\mathfrak{T}$, is actually contained in $\Omega^+_0$. Furthermore, the interior global attractor $\mathfrak{T}$ is invariant. By using the above result, we show that the compact global attractor $\mathfrak{T}=\{\mathcal{E}_1\}$. This completes the proof of Theorem 6.1.

---

7. Discussion

In this paper, we formulate a multi-strain partial differential equation (PDE) model describing the transmission dynamics of a vector-borne disease that incorporates both incubation age of the exposed hosts and infection age of the infectious hosts, respectively. The formulas for the reproduction number $\mathcal R^j_0$ of strain $j, ~j=1,\cdots,n$ are obtained from the biological meanings of models. And we define the basic number of the disease as the maximum of the reproduction numbers of each strain. We show that if $\mathcal R_0<1$, the disease-free equilibrium is locally and globally asymptotically stable. That means the disease dies out and the number of infected with each strain goes to zero. If $\mathcal R_0>1$, without loss of generality, assuming $\mathcal R^1_0=\max\{\mathcal R^1_0,\cdots,\mathcal R^n_0\}>1$, we show that the single-strain equilibrium $\mathcal E_1$ corresponding to strain one exists. The single-strain equilibrium $\mathcal E_1$ is locally asymptotically stable when $\mathcal R^1_0>1$ and $\mathcal R^i_0<\mathcal R^1_0, ~i=2,\cdots,n$.

The main purpose in this article is to extend the competitive exclusion result established by Bremermann and Thieme in [2], who using a multiple-strain ODE model derives that if multiple strains circulate in the population only the strain with the largest reproduction number persists, the strains with suboptimal reproduction numbers are eliminated. The proof of the competitive exclusion principle is based on the proof of the global stability of the single-strain equilibrium $\mathcal E_1$. We approach the result by using a Lyapunov function under a stronger condition that

$ \displaystyle\frac{\mathcal R^i_0}{\mathcal R^1_0}<\displaystyle\frac{b_i}{b_1}<1, i\neq1. $

(83)

Our results do not include the case of

$ \max\{\mathcal R^1_0,\cdots,\mathcal R^n_0\}=\mathcal R^1_0=\mathcal R^2_0=\cdots=\mathcal R^m_0>1, m\leq n, m\geq 2. $

According to Proposition 3.3 in [16], where the authors proved and simulated by data that if there is no mutation between two strains and if the basic reproduction numbers corresponding to the two strains are the same, then for the two strain epidemic model there exist many coexistence equilibria, we guess that the coexistence of multi-strains may occur and it is impossible for competitive exclusion in this case.

From the expression (6) of the basic reproduction number $\mathcal R^j_0$ corresponding to strain $j$ and the inequality $ \mathcal R^i_0/\mathcal R^1_0<b_i/b_1,i\neq 1$, it follows that

${r_i}<{r_1},$

where

$r_j=\frac{\beta_h^j}{\mu_v+\alpha^j_v},~\text{for} ~~j=1,2, \cdots, n.$

$r_j$ represents the transmission rate of an infectious vector with strain $j$ during its entire infectious period. The condition (83) implies that the following three inequalities hold at the same time,

$ \mathcal R^i_0<\mathcal R^1_0, {r_i}<{r_1}, {b_i}<{b_1}, i\neq1. $

Recall that $b_j$ denotes the probability of a given susceptible vector being transmitted by an infected host with strain $j$. Then the condition (83) for the occurrence of competition exclusion of strain $1$ means that the basic reproduction number corresponding to strain $1$, the transmission rate of an infectious vector with strain $1$ during its entire infectious period, and the probability of a given susceptible vector being transmitted by an infected host with strain $1$ are all biggest comparing to three quantities of other strains.

---

Acknowledgments

Y. Dang is supported by NSF of Henan Province 142300410350, Z. Qiu's research is supported by NSFC grants No. 11671206 and No. 11271190 and X. Li is supported by NSF of China grant 11271314 and Plan For Scientific Innovation Talent of Henan Province 144200510021. We are very grateful to two anonymous referees for their careful reading, valuable comments and helpful suggestions, which help us to improve the presentation of this work significantly.

---