An immuno-eco-epidemiological model of competition

ABSTRACT

This paper introduces a novel immuno-eco-epidemiological model of competition in which one of the species is affected by a pathogen. The infected individuals from species one are structured by time-since-infection and the within-host dynamics of the pathogen and the immune response is also modelled. A novel feature of the model is the impact of the species two numbers on the ability of species one to mount an immune response. The within-host model has three equilibria: an extinction equilibrium, pathogen-only equilibrium and pathogen and immune response equilibrium which exists if the immune response reproduction number ${\mathcal{R}}_0>1$. The extinction equilibrium is always unstable, the pathogen-only equilibrium is stable if ${\mathcal{R}}_0<1$, and the coexistence equilibrium is stable whenever it exists. The between-host competition model has six equilibria: an extinction equilibrium, three disease-free equilibria: species one-only equilibrium, species two-only equilibrium and a disease-free species coexistence equilibrium. There are also two disease-present equilibria: species one-only disease equilibrium and disease coexistence equilibrium. The existence and stability of these equilibria are governed by six reproduction numbers. Results show that for a non-fatal disease, the disease coexistence equilibrium is stable whenever it exists.

KEYWORDS:

• Immunology
• epidemiology
• ecology
• competition
• multiscale modeling
• immuno-eco-epidemiology

AMS SUBJECT CLASSIFICATION:

• 92D30

1. Introduction

Lotka–Volterra models have been successfully used to describe ecological interactions, such as competition and predation, since the beginning of the twentieth century. Lotka–Volterra models by themselves, however, did not address questions related to the interplay between host-pathogen and ecological interactions.

Approximately 75% of the recently emerging pathogens affecting humans come from animal origin [15]. Understanding the distribution of pathogens in wild and domestic animal populations is a first step in studying dangerous zoonotic microorganisms. Wild animals, however, rarely exist in isolation. They are subjected to basic ecological interactions, such as predation and competition [39]. There is a rich and growing literature on the interplay of community interactions and host-pathogen dynamics (see reviews in [29,35,49]). Early papers that incorporated infectious disease into classic models of competition and predation include in particular [4,28]. Parasites with multiple hosts can also lead to indirect competitive interactions (apparent competition) both by direct transmission [32] and by transmission via vectors [12]; parasitism can also modify the outcome of preexisting competitive interactions [10].

The interest in the interplay of infectious disease and community interactions, an area now termed eco-epidemiology, was enhanced by Chattopadhyay and Arino [16]. The introduction of disease in the competition between two species has led to destabilization in the dynamics [20]. More recently, disease infecting a predator has been found to lead to complex dynamics and chaos [52]. Most of these models consist of ordinary differential equations and do not address how ecological interactions affect the host's immune response and how the within-host pathogen reproduction and immune response affect the ecological interaction (but see [33,51] for the impact of immunity).

Bridging the immunological scale and the epidemiological scale has been fascinating mathematical biology researchers for more than a decade [21,30]. The dependence of the epidemiological reproduction number ${\mathcal{R}}_0$ on within-host pathogen load and immune responses is of key interest as it explains how the within-host dynamics affects the population-level persistence of the disease [19,31,34,36,40,50,53]. Martcheva studied how the pathogen load affects the population-level prevalence of HIV and found that a medication-mediated decrease in viral load increases the population-level prevalence of HIV [43]. This rather perplexing finding is well known among the public health authorities and is attributed to increasing the lifespan of HIV-infected individuals.

The importance of multi-scale immuno-epidemiological modelling is best highlighted by its role in studying evolution [46]. Gilchrist and Sasaki [27] were the first to address co-evolution but since then evolution of virulence has been attracting significant attention. Because evolution typically involves multiple strains, a number of approaches have been developed to handle the emergent complexity [1,2,5,17,18,25,26,37,38]. One key limiting case (and the assumption we make in our model) is that each infected host harbors only a single strain of the pathogen.

Mathematically speaking, several different ways have been proposed for linking the within-host and between-host scales through (single-strain) differential equation models on both scales. The simplest approach seems to result in a system of ODE models, part of which describes the within-host dynamics while another part describes the between-host dynamics [14,22,23]. These models often represent an environmentally driven disease and the within-host and between-host system are linked through the contaminated environment. The simplicity of the approach has multiple advantages and allows for significant analysis, which in turn leads to insight into biological implications. The other two modelling approaches use ODEs for the immunological model and PDEs for the epidemiological. In one of these approaches, perhaps the most complex one, the epidemiological model consists of physiologically structured PDEs in which the structural independent variables are the dynamical variables of the ODE immune model. The first such model was proposed by Martcheva and Pilyugin [44] (see also [41]) and since then has been improved upon by several authors [6,7,24,47]. This approach presents interesting mathematical challenges (such as the potential for measure-valued solutions) but somewhat restricts the incorporation of significant realism into the immune system. The second approach involving PDEs, originally proposed by Gilchrist and Sasaki [27], links the within-host model and the between-host model through an age-since-infection variable and dependence of the epidemic parameters on the within-host variables. This approach allows for incorporation of the necessary realism in the immunological and the epidemiological systems and has been studied mathematically more often [19,40,42,43,45,48]. This approach leads to what is called nested immuno-epidemiological models and it is the approach we will take in this article.

Little work has been done at the interface of eco-epidemiological and immuno-epidemiological modelling, although incorporating immune responses in the Lotka– Volterra predation or competition model is a natural step to take the community interactions one step further and study how ecological interactions affect the host's immune response and how the within-host pathogen reproduction and immune response affect the ecological interaction (but see [33,51] for the impact of immunity). The only article we are aware of studies the interplay between predation and within-host dynamics of disease in the prey [11]. In this article we introduce and investigate an eco-epidemiological model of competition with disease in species one and a nested immunological model of the within-host dynamics of the pathogen with the immune response. We further incorporate one further level of competition – we allow the population size of species two to affect the ability of species one to mount an immune response to the pathogen. We term these novel type of models immuno-eco-epidemiological models.

In the next section we introduce the within-host and the between-host eco-epidemiological model. In Section 3 we analyse the within-host model. In Section 4 we study the equilibria of the nested immuno-eco-epidemiological model and we define six multi-scale reproduction numbers that determine the existence of the equilibria. In Section 5 we study the stability of the equilibria of the nested immuno-eco-epidemiological model. Section 6 contains summary of our results.

2. An immuno-eco-epidemiological model

In this section we introduce a two-species ecological competition model, in which one of the species is infected by a pathogen, while the other is not affected by the same pathogen. Species one, that is infected by the pathogen, is structured by age-since-infection and its immune status is taken into account. The resulting immuno-eco-epidemiological model portrays the interaction of the pathogen with the immune system and how this interaction affects the global dynamics among the species. We also assume that the competition between the species exercises pressure on species one that affects its immune response to the pathogen. This model is a combination of ODEs and PDEs. The PDE model represents the population-level dynamics between the species. The host system has been classified into two species, where the number of species one is given by $N_1\left(t\right)$. Species one is in competition with species two whose numbers are given by $N_2\left(t\right)$. The following sets of equations represent the PDE model. (1) $\begin{array}{rl}{S}^{\prime }& =r_1\left(1-\frac{N_1+{\alpha }_{12}{N}_2}{K_1}\right)N_1-S\left(t\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i(\tau ,t)\mathrm{d}\tau -\mu S\left(t\right),\\ i_{\tau }(\tau ,t)+i_t(\tau ,t)& =-\mu i(\tau ,t)-\alpha \left(\tau \right)i(\tau ,t),\\ i(0,t)& =S\left(t\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i(\tau ,t)\mathrm{d}\tau ,\\ N_2^{\prime }\left(t\right)& =r_2\left(1-\frac{N_2+{\alpha }_{21}{N}_1}{K_2}\right)N_2-{\mu }_2{N}_2,\end{array}$(1)

$S\left(t\right)$ represents the number hosts of species one susceptible to the pathogen and $i(\tau ,t)$ is the corresponding number of infected hosts of species one at time t with age of infection τ. This definition leads to the following equation: $N_1\left(t\right)=S\left(t\right)+{\int }_0^{\mathrm{\infty }}i(\tau ,t)\mathrm{d}\tau .$

We assume that the disease is not transmitted vertically and for species one only susceptible individuals reproduce. The growth rate of the number of susceptible hosts is described by a logistic function where $r_1$ is the intrinsic growth and $K_1$ is the carrying capacity of species one. ${\alpha }_{12}$ represents the interference that comes from species two in the growth of species one population. ${\mu }_i$ is the natural death rate of species i and ${\beta }_1\left(\tau \right)$ represents the rate of infection for species one. α represents the rate at which the infected hosts are removed. The dynamics of species two have been modelled in a similar way using a logistic function to represent the growth of the population. $r_2$ is the intrinsic growth and $K_2$ is the carrying capacity of species two. ${\alpha }_{21}$ represents the interference of species one in the growth of the species two population. The PDE model above combines a Lotka–Volterra competition model with and SI epidemic model.

The ODE model represents the multiplication of virus/parasite inside a host. The within-host model captures the interplay of the virus with the immune system. We focus on these two components only as the main components of what we want to model within-host – the pathogen and the immune response, which is affected by the competition. The model is general enough to capture the within-host processes in a variety of diseases. ‘Predator–prey style’ immune models have been used widely in the literature to understand a variety of within-host and within-host-between-host processes [8]. In this ODE system, $V\left(\tau \right)$ represents the number of virus particles and $z\left(\tau \right)$ represents the healthy immune cells with age of infection τ. The following ODE model gives the dynamics in the number of the healthy immune cells and the number of virus particles. (2) $\begin{array}{rl}{V}^{\prime }\left(\tau \right)& =qV\left(\tau \right)\left(1-\frac{V\left(\tau \right)}{K_v}\right)-aV\left(\tau \right)z\left(\tau \right),\\ z^{\prime }\left(\tau \right)& =\frac{b_{0}V\left(\tau \right)z\left(\tau \right)}{N_2\left(t\right)+D}-dz\left(\tau \right).\end{array}$(2)

We assume that the growth of the virus population occurs in a logistic fashion, where $K_v$ is the carrying capacity of the virus inside a host and q is the intrinsic growth rate. a is the elimination rate of the virus particles by the healthy immune system of the host. Healthy immune cells are represented by $z\left(\tau \right)$. We use immune cells as a proxy to antibodies which mark the pathogen for destruction. Furthermore, we assume that the immune-response cells are cleared at a rate d. The constant $b_0/D$ is the immune response activation rate in the absence of interference from species two. D represents a scaling factor. We note that the immune system contains as a parameter the number of species two $N_2\left(t\right)$ which affects negatively the immune response. We assume that being subjected to competition with species two, species one is stressed to an extent that modifies it's internal ability to mount an effective immune response to the pathogen. The impact of stress on immune system functioning is well recognized [3].

The parameters of the population-level model are linked to the dynamic variables of the within-host model in the following way. It has been assumed that the rate of infection ${\beta }_1\left(\tau \right)$ and the rate of removal of the infected hosts $\alpha \left(\tau \right)$ are proportional to the within-host viral load. (3) ${\beta }_1\left(\tau \right)=c_{1}V\left(\tau \right)\mathrm{and}\alpha \left(\tau \right)=c_{2}V\left(\tau \right),$(3) where $c_1$ and $c_2$ are constants of proportionality. The assumption for linking ${\beta }_1\left(\tau \right)$ linearly to the pathogen load is common [27,50]; however, newer results obtained from comparison with data suggest that a better linking is a Hill function [42]. Nonetheless, we stay with the simplest case. The linking for $\alpha \left(\tau \right)$ typically reflects also the immune response [27,50]; however, the presence of the immune response creates a trade-off for the impact of the competition of species two. In this article, we are only trying to investigate the negative effect due to stress. We notice that we have assumed that the disease has no recovery, which is characteristic of many diseases of wild animals (e.g. Feline Immunodeficiency Virus (FIV) in lions).

3. Analysis of the ODE model

The within-host ODE model describes the interplay between the virus and the host immune system. To understand better its dynamics, we begin by defining the equilibria of the model and then investigate the stability of each equilibrium. The equilibria are time-independent solutions of the ODE model, obtained by equating each of the equations in the ODE model to 0.

The question of equilibria in the immunological model is not completely trivial, as the model involves the time-dependent size of species two $N_2\left(t\right)$. We would investigate the equilibria at some fixed value of $N_2\left(t\right)=N_2^{\ast }$ which is to be determined later. Assume $V\left(\tau \right)=V^{\ast }$ and $z\left(\tau \right)=z^{\ast }$ be the equilibrium points. (4) $\begin{array}{rl}0& =q{V}^{\ast }\left(1-\frac{V^{\ast }}{K_v}\right)-a{V}^{\ast }{z}^{\ast },\\ 0& =\frac{b_0{V}^{\ast }{z}^{\ast }}{N_2^{\ast }+D}-d{z}^{\ast }.\end{array}$(4)

Solving the above set of equations leads to the following equilibria of this system, given as ordered pairs $(V^{\ast },z^{\ast })$. The immunological extinction equilibrium ${\mathcal{E}}_0=(0,0)$ and the endemic equilibrium ${\mathcal{E}}_1=(K_v,0)$ always exist. The co-existence equilibrium is given by ${\mathcal{E}}^{\ast }=(\hat{V},\hat{z})$ where $\hat{V}=\frac{d}{b_0}(N_2^{\ast }+D)\mathrm{and}\hat{z}=\frac{q}{a}\left(1-\frac{1}{{\mathcal{R}}_0}\right).$

We see that the equilibrial level of species two $N_2$ suppresses the immune response and supports higher level of the virus. We observe that this equilibrium exists if and only if ${\mathcal{R}}_0>1$. The immune-response reproduction number ${\mathcal{R}}_0$ is given by ${\mathcal{R}}_0=\frac{K_v{b}_0}{d(N_2^{\ast }+D)}.$ This number gives the number of secondary immune particles that one immune particle stimulates when the virus is at carrying capacity and species two is at equilibrial value $N_2^{\ast }$.

3.1. Local stability of the equilibria

We explore the stability of the equilibria defined in the previous subsection. The presence of time-dependent species two size $N_2\left(t\right)$ in system (2(2) $\begin{array}{rl}{V}^{\prime }\left(\tau \right)& =qV\left(\tau \right)\left(1-\frac{V\left(\tau \right)}{K_v}\right)-aV\left(\tau \right)z\left(\tau \right),\\ z^{\prime }\left(\tau \right)& =\frac{b_{0}V\left(\tau \right)z\left(\tau \right)}{N_2\left(t\right)+D}-dz\left(\tau \right).\end{array}$(2) ) is a significant complication. We assume that the equilibrium in the ODE immune system occurs after the epidemic system has stabilized at equilibrium. Therefore, $N_2\left(t\right)\approx N_2^{\ast }$ where $N_2^{\ast }$ is any equilibrial value of epidemic system. We believe the analysis of this section can be carried out with $N_2\left(t\right)$ rather than $N_2^{\ast }$; however, the dependence of $N_2$ on t may have implications for the conclusions. A more careful analysis of this scenario will be considered in a future work.

We consider the following asymptotically autonomous ODE model: (5) $\begin{array}{rl}{V}^{\prime }\left(\tau \right)& =qV\left(\tau \right)\left(1-\frac{V\left(\tau \right)}{K_v}\right)-aV\left(\tau \right)z\left(\tau \right),\\ z^{\prime }\left(\tau \right)& =\frac{b_{0}V\left(\tau \right)z\left(\tau \right)}{N_2^{\ast }+D}-dz\left(\tau \right).\end{array}$(5)

Theorem 3.1.

The extinction equilibrium ${\mathcal{E}}_0$ is always unstable. The pathogen-only equilibrium ${\mathcal{E}}_1$ is locally asymptotically stable if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$. The coexistence equilibrium ${\mathcal{E}}^{\ast }$ is locally asymptotically stable if ${\mathcal{R}}_0>1$.

Proof.

The local stability of an equilibrium is given by the Jacobian of the system (5(5) $\begin{array}{rl}{V}^{\prime }\left(\tau \right)& =qV\left(\tau \right)\left(1-\frac{V\left(\tau \right)}{K_v}\right)-aV\left(\tau \right)z\left(\tau \right),\\ z^{\prime }\left(\tau \right)& =\frac{b_{0}V\left(\tau \right)z\left(\tau \right)}{N_2^{\ast }+D}-dz\left(\tau \right).\end{array}$(5) ). The Jacobian of the above system (5(5) $\begin{array}{rl}{V}^{\prime }\left(\tau \right)& =qV\left(\tau \right)\left(1-\frac{V\left(\tau \right)}{K_v}\right)-aV\left(\tau \right)z\left(\tau \right),\\ z^{\prime }\left(\tau \right)& =\frac{b_{0}V\left(\tau \right)z\left(\tau \right)}{N_2^{\ast }+D}-dz\left(\tau \right).\end{array}$(5) ) is: $J=\left(\begin{array}{cc}q-\frac{2q{V}^{\ast }}{K_v}-a{z}^{\ast }& -a{V}^{\ast }\\ \frac{b_0{z}^{\ast }}{N_2^{\ast }+D}& \frac{b_0{V}^{\ast }}{N_2^{\ast }+D}-d\end{array}\right)$ At the trivial equilibrium ${\mathcal{E}}_0=(0,0)$ the Jacobian reduces to $J=\left(\begin{array}{cc}q& 0\\ 0& -d\end{array}\right).$ The two roots of this matrix are q and $-d$. Since q>0, the extinction equilibrium is always unstable. At the endemic equilibrium ${\mathcal{E}}_1$ the Jacobian reduces to the following matrix $J=\left(\begin{array}{cc}-q& -a{K}_v\\ 0& \frac{b_0{K}_v}{N_2^{\ast }+D}-d\end{array}\right)$ The roots of this Jacobian matrix are $-q$ and $b_0{K}_v/(N_2^{\ast }+D)-d$. The roots are both negative if and only if $b_0{K}_v/(N_2^{\ast }+D)-d<0$ which is negative if and only if ${\mathcal{R}}_0<1$. When ${\mathcal{R}}_0>1$, we have that $b_0{K}_v/(N_2^{\ast }+D)-d>0$ and the endemic, immune-response-free equilibrium ${\mathcal{E}}_1$ is unstable.

Now, we will explore the co-existence equilibrium ${\mathcal{E}}^{\ast }$. At this equilibrium, $V^{\ast }$ and $z^{\ast }$ satisfy the following set of equations: (6) $\begin{array}{rl}q\left(1-\frac{V^{\ast }}{K_v}\right)-a{z}^{\ast }& =0,\\ \frac{b_0{V}^{\ast }}{N_2^{\ast }+D}-d& =0.\end{array}$(6)

This simplifies the Jacobian to the following form: $J=\left(\begin{array}{cc}-\frac{q{V}^{\ast }}{K_v}& -a{V}^{\ast }\\ \frac{b_0{z}^{\ast }}{N_2^{\ast }+D}& 0\end{array}\right).$

Exploring the properties of this Jacobian matrix, we observe that $\mathrm{Trace}\left(J\right)<0$ and $\mathrm{Det}\left(J\right)=b_0{z}^{\ast }/(N_2^{\ast }+D)>0$. This shows that ${\mathcal{E}}^{\ast }$ is locally asymptotically stable.

3.2. Global stability of equilibria

In the case ${\mathcal{R}}_0<1$, the only locally stable equilibrium is ${\mathcal{E}}_1$. To see the global stability of ${\mathcal{E}}_1$, notice that from the first equation in Equation (2(2) $\begin{array}{rl}{V}^{\prime }\left(\tau \right)& =qV\left(\tau \right)\left(1-\frac{V\left(\tau \right)}{K_v}\right)-aV\left(\tau \right)z\left(\tau \right),\\ z^{\prime }\left(\tau \right)& =\frac{b_{0}V\left(\tau \right)z\left(\tau \right)}{N_2\left(t\right)+D}-dz\left(\tau \right).\end{array}$(2) ) we have that $\underset{\tau }{lim sup}V\left(\tau \right)\le K_v.$ Then from the second equation, we have $z^{\prime }\left(\tau \right)\le d({\mathcal{R}}_0-1)z\left(\tau \right).$ Hence $\underset{\tau \to \mathrm{\infty }}{lim}z\left(\tau \right)=0$. From here, it is not hard to prove that $V\left(\tau \right)\to K_v$ which establishes the global stability of ${\mathcal{E}}_1$.

In the case of ${\mathcal{R}}_0>1$, only equilibrium ${\mathcal{E}}^{\ast }$ is locally asymptotically stable. To its see global stability for system (2(2) $\begin{array}{rl}{V}^{\prime }\left(\tau \right)& =qV\left(\tau \right)\left(1-\frac{V\left(\tau \right)}{K_v}\right)-aV\left(\tau \right)z\left(\tau \right),\\ z^{\prime }\left(\tau \right)& =\frac{b_{0}V\left(\tau \right)z\left(\tau \right)}{N_2\left(t\right)+D}-dz\left(\tau \right).\end{array}$(2) ), we use the fact that the coexistence equilibrium is locally stable whenever it exists. The term $N_2\left(t\right)$ is treated as a parameter. We use Dulac's criterium on the open first quadrant ${\mathbf{R}}_{+}^2$ to rule out periodic orbits. Let $f(V,z)=qV\left(\tau \right)\left(1-\frac{V\left(\tau \right)}{K_v}\right)-aV\left(\tau \right)z\left(\tau \right)$ and $g(V,z)=\frac{b_{0}V\left(\tau \right)z\left(\tau \right)}{N_2\left(t\right)+D}-dz\left(\tau \right).$

Define Dulac's auxiliary function $\psi =1/Vz$. Then we find that $\frac{\mathrm{\partial }}{\mathrm{\partial }V}\left(f\psi \right)+\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(g\psi \right)=-\frac{q}{K_{v}z}<0.$ By Dulac's criterion, there are no periodic solutions. By Poincare Bendixon theorem, the solutions converges to an equilibrium. When ${\mathcal{R}}_0>1$, both ${\mathcal{E}}_0$ and ${\mathcal{E}}_1$ are unstable. Since we have already proved the local stability of ${\mathcal{E}}^{\ast }$, by Poincare Bendixon theorem, we conclude that the co-existence equilibrium ${\mathcal{E}}^{\ast }$ is globally asymptotically stable.

From this analysis, we observe that the ODE model has two regimes: if ${\mathcal{R}}_0<1$ the virus persist but the immune system does not get activated. If ${\mathcal{R}}_0>1$, both the virus and the immune system persist. The system (2(2) $\begin{array}{rl}{V}^{\prime }\left(\tau \right)& =qV\left(\tau \right)\left(1-\frac{V\left(\tau \right)}{K_v}\right)-aV\left(\tau \right)z\left(\tau \right),\\ z^{\prime }\left(\tau \right)& =\frac{b_{0}V\left(\tau \right)z\left(\tau \right)}{N_2\left(t\right)+D}-dz\left(\tau \right).\end{array}$(2) ) does not predict that the virus can be cleared. Thus, the system models chronic infections. This is reasonable since animal populations do not receive treatment and many infections for them are not cleared.

4. Equilibria of the PDE model

In this section we will explore the equilibria of the PDE model. We incorporate the within-host model with its dynamical properties; that is the within-host model is not at equilibrium.

We begin by introducing some notation that we will use throughout the article. The probability of survival of species one as infected until time τ post infection is given by $\mathrm{\Pi }\left(\tau \right)={\mathrm{e}}^{-{\int }_0^{\tau }(\mu +\alpha (r\left)\right)\mathrm{d}r}.$ We note that the probability of survival depends on the within-host viral load $V\left(\tau \right)$ through Equation (3(3) ${\beta }_1\left(\tau \right)=c_{1}V\left(\tau \right)\mathrm{and}\alpha \left(\tau \right)=c_{2}V\left(\tau \right),$(3) ). Furthermore, we denote by $G\left(N_2\right)={\int }_0^{\mathrm{\infty }}\mathrm{\Pi }\left(\theta \right)\mathrm{d}\theta .$ We note that G depends on $N_2$ through the within-host system. G in the text will also stand for $G\left(0\right)$. We note that $\mu G\left(N_2\right)<1$ gives the probability of dying from natural causes while infectious. We further adopt the following two notations: $f_s\left(N_2^{\ast }\right)=\frac{1}{{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\mathrm{\Pi }\left(\theta \right)\mathrm{d}\theta }$ and $Q=K_2\left(1-\frac{{\mu }_2}{r_2}\right).$ The existence and stability of the equilibria of the immuno-eco-epidemiological model (1(1) $\begin{array}{rl}{S}^{\prime }& =r_1\left(1-\frac{N_1+{\alpha }_{12}{N}_2}{K_1}\right)N_1-S\left(t\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i(\tau ,t)\mathrm{d}\tau -\mu S\left(t\right),\\ i_{\tau }(\tau ,t)+i_t(\tau ,t)& =-\mu i(\tau ,t)-\alpha \left(\tau \right)i(\tau ,t),\\ i(0,t)& =S\left(t\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i(\tau ,t)\mathrm{d}\tau ,\\ N_2^{\prime }\left(t\right)& =r_2\left(1-\frac{N_2+{\alpha }_{21}{N}_1}{K_2}\right)N_2-{\mu }_2{N}_2,\end{array}$(1) ) depend on a number of reproduction numbers and invasion reproductions numbers. We list these key quantities in the following definition.

Definition 4.1.

We define the reproduction numbers used in this model as follows.

1. The basic reproduction number for species one $N_1$ is defined as

   1. Species one only disease-free reproduction number $R_0^{N_{10}}=r_1/\mu$

   2. Species one only susceptible-disease reproduction number $R_0^{N_1}=r_1/(\mu +(r_1/K_1\left)S^{\ast }\right)$ where $S^{\ast }=\frac{1}{{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\mathrm{\Pi }\left(\theta \right)\mathrm{d}\theta }$

2. The basic reproduction number for species two $N_2$ is defined as $R_0^{N_2}=\frac{r_2}{{\mu }_2}.$

3. The endemic reproduction number for species one is given by $R_{N1}=\frac{r_2}{{\mu }_2+\frac{r_2{\alpha }_{21}{N}_1^{\ast }}{K_2}}$ where $N_1^{\ast }$ is the solution of equation for the species one endemic equilibrium with $N_2^{\ast }=0$. This equation will be defined in the following section.

4. The endemic reproduction number for species two is defined as $R_{N2}=\frac{r_1}{\frac{r_1{\alpha }_{12}{K}_2}{K_1}\left(1-\frac{{\mu }_2}{r_2}\right)+\mu }.$

5. The co-existence species one-species two, no disease reproduction number is given by $R_{N1N2}=\frac{r_2}{\frac{r_2{\alpha }_{21}{K}_1}{K_2}\left(1-\frac{\mu }{r_1}\right)+{\mu }_2}.$

The equilibria are obtained by equating the time derivatives to 0. Since each of the variables is independent of time, we let $S\left(t\right)=S^{\ast },i(\tau ,t)=i^{\ast }\left(\tau \right),N_1\left(t\right)=N_1^{\ast },N_2\left(t\right)=N_2^{\ast }.$ This leads to the following set of equations for the equilibria: (7) $\begin{array}{rl}0& =r_1\left(1-\frac{N_1^{\ast }+{\alpha }_{12}{N}_2^{\ast }}{K_1}\right)N_1^{\ast }-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau -\mu S^{\ast },\\ i_{\tau }^{\ast }\left(\tau \right)& =-\mu i^{\ast }\left(\tau \right)-\alpha \left(\tau \right)i^{\ast }\left(\tau \right),\\ i^{\ast }\left(0\right)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau ,\\ 0& =r_2\left(1-\frac{N_2^{\ast }+{\alpha }_{21}{N}_1^{\ast }}{K_2}\right)N_2^{\ast }-{\mu }_2{N}_2^{\ast },\\ N_1^{\ast }& =S^{\ast }+{\int }_0^{\mathrm{\infty }}{i}^{\ast }\left(\tau \right)\mathrm{d}\tau .\end{array}$(7)

We will obtain the equilibrium as an ordered set $(S^{\ast },i^{\ast }(\tau ),N_2^{\ast })$, where the elements in the ordered set satisfy the relationship, $N_1^{\ast }=S^{\ast }+{\int }_0^{\mathrm{\infty }}{i}^{\ast }\left(\tau \right)\mathrm{d}\tau$.

4.1. The eco-epidemiological extinction equilibrium

The eco-epidemiological extinction equilibrium $E_0=(0,0,0)$ always exits, since it always satisfies Equations (7(7) $\begin{array}{rl}0& =r_1\left(1-\frac{N_1^{\ast }+{\alpha }_{12}{N}_2^{\ast }}{K_1}\right)N_1^{\ast }-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau -\mu S^{\ast },\\ i_{\tau }^{\ast }\left(\tau \right)& =-\mu i^{\ast }\left(\tau \right)-\alpha \left(\tau \right)i^{\ast }\left(\tau \right),\\ i^{\ast }\left(0\right)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau ,\\ 0& =r_2\left(1-\frac{N_2^{\ast }+{\alpha }_{21}{N}_1^{\ast }}{K_2}\right)N_2^{\ast }-{\mu }_2{N}_2^{\ast },\\ N_1^{\ast }& =S^{\ast }+{\int }_0^{\mathrm{\infty }}{i}^{\ast }\left(\tau \right)\mathrm{d}\tau .\end{array}$(7) ).

4.2. Semi-trivial equilibria

Case 1: Disease-free Equilibria

(a) We observe that there exists a semi-trivial species one equilibrium in the form $E_1=(K_1(r_1-\mu )/r_1,0,0)$. Note that this equilibrium exists, if and only if $R_0^{N_{10}}>1$.

(b) Now, we will explore the semi-trivial species two equilibrium where species one is not present. Let us assume $N_1^{\ast }=0$. It is clear from this assumption that $N_1^{\ast }=0\Rightarrow S^{\ast }=0,i^{\ast }\left(\tau \right)=0$. Hence, the species two equilibrium is defined as $E_2=(0,0,N_2^{\ast })$. To find $N_2^{\ast }$, we have to solve the following equation: (8) $0=r_2\left(1-\frac{N_2^{\ast }}{K_2}\right)-{\mu }_2.$(8) Solving for $N_2^{\ast }$, we obtain $N_2^{\ast }=K_2\left(1-\frac{{\mu }_2}{r_2}\right).$ Note that this equilibrium exists if and only if the reproduction number of species two is larger than one: $R_0^{N_2}>1$.

Case 2: Semi-trivial endemic equilibria. Species one endemic equilibrium exists if and only if species one-disease reproduction number is larger than one. This result is given by the following theorem:

Theorem 4.1.

There exists an endemic equilibrium $E_{N1}$ in the form $E_{N1}=(S^{\ast },i^{\ast }(\tau ),0)$ if and only if $R_{N1}>1$.

Proof.

We examine the boundary equilibria when $N_2^{\ast }=0$ and $N_1^{\ast }\ne 0$. In this case the system for the equilibria (7(7) $\begin{array}{rl}0& =r_1\left(1-\frac{N_1^{\ast }+{\alpha }_{12}{N}_2^{\ast }}{K_1}\right)N_1^{\ast }-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau -\mu S^{\ast },\\ i_{\tau }^{\ast }\left(\tau \right)& =-\mu i^{\ast }\left(\tau \right)-\alpha \left(\tau \right)i^{\ast }\left(\tau \right),\\ i^{\ast }\left(0\right)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau ,\\ 0& =r_2\left(1-\frac{N_2^{\ast }+{\alpha }_{21}{N}_1^{\ast }}{K_2}\right)N_2^{\ast }-{\mu }_2{N}_2^{\ast },\\ N_1^{\ast }& =S^{\ast }+{\int }_0^{\mathrm{\infty }}{i}^{\ast }\left(\tau \right)\mathrm{d}\tau .\end{array}$(7) ) reduces to solving the following set of equations: (9) $\begin{array}{rl}0& =r_1\left(1-\frac{N_1^{\ast }}{K_1}\right)N_1^{\ast }-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau -\mu S^{\ast },\\ i_{\tau }^{\ast }\left(\tau \right)& =-\mu i^{\ast }\left(\tau \right)-\alpha \left(\tau \right)i^{\ast }\left(\tau \right),\\ i^{\ast }\left(0\right)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau ,\\ N_1^{\ast }& =S^{\ast }+{\int }_0^{\mathrm{\infty }}{i}^{\ast }\left(\tau \right)\mathrm{d}\tau .\end{array}$(9)

Solving for $i^{\ast }\left(\tau \right)$, we obtain $i^{\ast }\left(\tau \right)=i^{\ast }\left(0\right){\mathrm{e}}^{-{\int }_0^{\tau }(\mu +\alpha (r\left)\right)\mathrm{d}r}=i^{\ast }\left(0\right)\mathrm{\Pi }\left(\tau \right).$

Recall that $\mathrm{\Pi }\left(\tau \right)={\mathrm{e}}^{-{\int }_0^{\tau }(\mu +\alpha (r\left)\right)\mathrm{d}r}$. Using this expression in the boundary condition, we obtain the value of $S^{\ast }$: $S^{\ast }=\frac{1}{{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\mathrm{\Pi }\left(\theta \right)\mathrm{d}\theta }.$ Using the equation $N_1^{\ast }=S^{\ast }+{\int }_0^{\mathrm{\infty }}{i}^{\ast }\left(\tau \right)\mathrm{d}\tau$, we obtain the value of $i^{\ast }\left(0\right)$: $i^{\ast }\left(0\right)=\frac{N_1^{\ast }-S^{\ast }}{G}.$ Recall that $G={\int }_0^{\mathrm{\infty }}\mathrm{\Pi }\left(\theta \right)\mathrm{d}\theta$. Thus, we can express every term in the equilibrium as a function of $N_1^{\ast }$. Substituting all these in the first equation of (9(9) $\begin{array}{rl}0& =r_1\left(1-\frac{N_1^{\ast }}{K_1}\right)N_1^{\ast }-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau -\mu S^{\ast },\\ i_{\tau }^{\ast }\left(\tau \right)& =-\mu i^{\ast }\left(\tau \right)-\alpha \left(\tau \right)i^{\ast }\left(\tau \right),\\ i^{\ast }\left(0\right)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau ,\\ N_1^{\ast }& =S^{\ast }+{\int }_0^{\mathrm{\infty }}{i}^{\ast }\left(\tau \right)\mathrm{d}\tau .\end{array}$(9) ) we obtain the following quadratic equation for $N_1^{\ast }$: $\frac{r_1{N}_1^{\ast 2}}{K_1}+\left(\frac{1}{G}-r_1\right)N_1^{\ast }+\left(\mu -\frac{1}{G}\right)S^{\ast }=0.$ Note that $G={\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{-{\int }_0^{\theta }(\mu +\alpha (r\left)\right)\mathrm{d}r}\mathrm{d}\theta \le {\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{-\mu \theta }\mathrm{d}\theta =\frac{1}{\mu }\Rightarrow \mu \le \frac{1}{G}.$ Hence, the constant term in the quadratic polynomial is negative. The leading co-efficient is positive. By intermediate value theorem, it can be proved that there always exists a unique positive solution for $N_1^{\ast }$. It remains to prove that $N_1^{\ast }>S^{\ast }$.

To see that $N_1^{\ast }>S^{\ast }$, we rewrite the polynomial in a different form. We can define $N=N_1^{\ast }$ as the solution of the polynomial $H\left(N\right)=0$ where $H\left(N\right)=r_1\left(1-\frac{N}{K_1}\right)N-\mu S^{\ast }-\frac{N-S^{\ast }}{G}.$

We observe at $N=S^{\ast }$, we have $H\left(S^{\ast }\right)=S^{\ast }\left(r_1-\mu -\frac{r_1}{K_1}{S}^{\ast }\right).$ From the structure of the endemic reproduction number $R_{N1}$, we have $R_{N1}>1\Rightarrow H\left(S^{\ast }\right)>0$. Note that as $N\to \mathrm{\infty }$, $H\left(N\right)$ approaches a negative quantity. Hence, there exists a positive solution $N=N^{\ast }$ such that $N^{\ast }>S^{\ast }$. Note that in this case $i^{\ast }\left(0\right)=(N^{\ast }-S^{\ast })/G>0$. This proves the existence of an endemic species one equilibrium $E_{N1}=(S^{\ast },i^{\ast }(\tau ),0)$ where $i^{\ast }\left(\tau \right)=\left((N^{\ast }-S^{\ast })/G\right)\mathrm{\Pi }\left(\tau \right)$.

4.3. Species co-existence equilibria

There are several equilibria where the two species coexist. These are species one-species two no disease equilibrium and species one-disease-species two equilibrium.

Case 1: Species one- species two, no disease equilibrium.

This equilibrium is obtained from solving the following sets of linear equation in two variables. (10) $\begin{array}{rl}{N}_1+{\alpha }_{12}{N}_2& =K_1\left(1-\frac{\mu }{r_1}\right),\\ {\alpha }_{21}{N}_1+N_2& =K_2\left(1-\frac{{\mu }_2}{r_2}\right).\end{array}$(10)

Lemma 4.1.

System of Equations (10(10) $\begin{array}{rl}{N}_1+{\alpha }_{12}{N}_2& =K_1\left(1-\frac{\mu }{r_1}\right),\\ {\alpha }_{21}{N}_1+N_2& =K_2\left(1-\frac{{\mu }_2}{r_2}\right).\end{array}$(10) ) has a solution when $R_{N2}>1$ and $R_{N1N2}>1$.

Proof of Lemma.

Proof of Lemma

Note $R_{N1N2}>1\Rightarrow K_2(1-{\mu }_2/r_2)-{\alpha }_{21}{K}_1(1-\mu /r_1)>0$ and $R_{N2}>1\Rightarrow K_1(1-\mu /r_1)-{\alpha }_{12}{K}_2(1-{\mu }_2/r_2)>0$.

From these two inequalities, it also follows that ${\alpha }_{12}{\alpha }_{21}<1$. Hence from Cramer's rule, we have positive solution for the set of linear equation in (10(10) $\begin{array}{rl}{N}_1+{\alpha }_{12}{N}_2& =K_1\left(1-\frac{\mu }{r_1}\right),\\ {\alpha }_{21}{N}_1+N_2& =K_2\left(1-\frac{{\mu }_2}{r_2}\right).\end{array}$(10) ). This completes the proof.

Case 2: Disease-present coexistence equilibrium $E^{\ast }=(S^{\ast },i^{\ast }(\tau ),N_2^{\ast })$.

Now, we consider the disease-coexistence equilibrium. The existence is given by the following theorem.

Theorem 4.2.

Assume $R_0^{N_2}>1,$ $R_{N1N2}>1$ and $R_{N1}>1$. Then, there exists $E^{\ast }$.

Proof.

The system that has to be solved to obtain this equilibrium is given by (11) $\begin{array}{rl}0& =r_1\left(1-\frac{N_1^{\ast }+{\alpha }_{12}{N}_2^{\ast }}{K_1}\right)N_1^{\ast }-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau -\mu S^{\ast },\\ i_{\tau }^{\ast }\left(\tau \right)& =-\mu i^{\ast }\left(\tau \right)-\alpha \left(\tau \right)i^{\ast }\left(\tau \right),\\ i^{\ast }\left(0\right)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau ,\\ 0& =r_2\left(1-\frac{N_2^{\ast }+{\alpha }_{21}{N}_1^{\ast }}{K_2}\right)N_2^{\ast }-{\mu }_2{N}_2^{\ast },\\ N_1^{\ast }& =S^{\ast }+{\int }_0^{\mathrm{\infty }}{i}^{\ast }\left(\tau \right)\mathrm{d}\tau .\end{array}$(11)

Solving for $i^{\ast }\left(\tau \right)$ we obtain $i^{\ast }\left(\tau \right)=i^{\ast }\left(0\right)\mathrm{\Pi }\left(\tau \right)$. Using the boundary conditions, we obtain $S^{\ast }=\frac{1}{{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\mathrm{\Pi }\left(\theta \right)\mathrm{d}\theta }=f_s\left(N_2^{\ast }\right).$ The other equations can be solved to express $i^{\ast }\left(0\right)$ is terms of $N_1^{\ast }$ and $N_2^{\ast }$ as follows: $i^{\ast }\left(0\right)=r_1\left(1-\frac{N_1^{\ast }+{\alpha }_{12}{N}_2^{\ast }}{K_1}\right)N_1^{\ast }-\mu S^{\ast }.$ Using the last equation in (11(11) $\begin{array}{rl}0& =r_1\left(1-\frac{N_1^{\ast }+{\alpha }_{12}{N}_2^{\ast }}{K_1}\right)N_1^{\ast }-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau -\mu S^{\ast },\\ i_{\tau }^{\ast }\left(\tau \right)& =-\mu i^{\ast }\left(\tau \right)-\alpha \left(\tau \right)i^{\ast }\left(\tau \right),\\ i^{\ast }\left(0\right)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right)i^{\ast }\left(\tau \right)\mathrm{d}\tau ,\\ 0& =r_2\left(1-\frac{N_2^{\ast }+{\alpha }_{21}{N}_1^{\ast }}{K_2}\right)N_2^{\ast }-{\mu }_2{N}_2^{\ast },\\ N_1^{\ast }& =S^{\ast }+{\int }_0^{\mathrm{\infty }}{i}^{\ast }\left(\tau \right)\mathrm{d}\tau .\end{array}$(11) ) we can express $N_1^{\ast }$ as a function of $N_2^{\ast }$ as $N_1^{\ast }=\frac{K_2\left(1-\frac{{\mu }_2}{r_2}\right)-N_2}{{\alpha }_{21}}.$ This expresses every term as a function of $N_2^{\ast }$. Thus, the solution $N_2^{\ast }$ is obtained by solving the equation $\mathcal{G}\left(N\right)=0$ where $\mathcal{G}\left(N\right)=\mu f_s\left(N\right)+\frac{\frac{Q-N}{{\alpha }_{21}}-f_s\left(N\right)}{G}-r_1\left(\frac{Q-N}{{\alpha }_{21}}\right)\left(1-\frac{\frac{Q-N}{{\alpha }_{21}}+{\alpha }_{12}N}{K_1}\right)$ and $Q=K_2\left(1-\frac{{\mu }_2}{r_2}\right)=K_2\left(1-\frac{1}{R_0^{N_2}}\right).$ Note that for $R_0^{N_2}>1$ we have Q>0.

Since $\alpha \left(\tau \right)=c_{2}V\left(\tau \right)$ and $V\left(\tau \right)$, as solution of Equation (5(5) $\begin{array}{rl}{V}^{\prime }\left(\tau \right)& =qV\left(\tau \right)\left(1-\frac{V\left(\tau \right)}{K_v}\right)-aV\left(\tau \right)z\left(\tau \right),\\ z^{\prime }\left(\tau \right)& =\frac{b_{0}V\left(\tau \right)z\left(\tau \right)}{N_2^{\ast }+D}-dz\left(\tau \right).\end{array}$(5) ), depends on $N_2^{\ast }$, we can consider $G=G\left(N_2^{\ast }\right)$, a function of $N_2^{\ast }$.

At N=Q in the characteristic equation, we have (12) $\mathcal{G}\left(Q\right)=f_s\left(Q\right)\left(\mu -\frac{1}{G\left(Q\right)}\right).$(12)

From the argument presented in the Proof of Theorem 3.1, we can claim that $\mu -1/G\left(Q\right)<0$. As $f_s\left(Q\right)>0$, it clearly follows that $\mathcal{G}\left(Q\right)<0$.

Observe that $\mathcal{G}\left(0\right)=\frac{r_1{Q}^2}{{\alpha }_{21}^2{K}_1}+\left[\frac{1}{G\left(0\right)}-r_1\right]\frac{Q}{{\alpha }_{21}}+\left(\mu -\frac{1}{G\left(0\right)}\right)S.$ Hence, $\mathcal{G}\left(0\right)=\left(\mu -\frac{1}{G\left(0\right)}\right)S+\frac{Q}{{\alpha }_{21}}\left[\frac{1}{G\left(0\right)}-r_1+\frac{r_1}{K_1}\frac{Q}{{\alpha }_{21}}\right].$ From the fact that $R_{N1N2}>1$, we have $\frac{Q{r}_1}{K_1{\alpha }_{21}}-r_1>-\mu .$ Hence, $\mathcal{G}\left(0\right)>\left(\mu -\frac{1}{G\left(0\right)}\right)\left(S-\frac{Q}{{\alpha }_{21}}\right).$ The first multiple in the product is negative. From $R_{N1}>1$, we have that $\frac{Q}{{\alpha }_{21}}>N_1>S.$ Therefore, the second factor is also negative. We have that $\mathcal{G}\left(0\right)>0$. By the intermediate value theorem, we can now claim that there exist a solution $N=N_2^{\ast }$ to the equation $\mathcal{G}\left(N\right)=0$ such that $0<N_2^{\ast }<Q$. From this condition and from the definition of $N_1^{\ast }=(Q-N_2^{\ast })/{\alpha }_{21}$, it is clear that there exist a positive solution $N_1^{\ast }$ which clearly defines this equilibrium.

5. Local stability of the equilibria

In this section, we will explore the stability of the different equilibria. In general, if the equilibrium is defined as the ordered n-tuple $(S^{\ast },i^{\ast }(\tau ),N_2^{\ast })$, where $N_1^{\ast }=S^{\ast }+{\int }_0^{\mathrm{\infty }}{i}^{\ast }\left(\tau \right)\mathrm{d}\tau ,$ we will observe the change in the dynamical system when the system is perturbed from its equilibrium state. The following defines the perturbations applied to the system. (13) $\begin{array}{rl}S\left(t\right)& =S^{\ast }+s\left(t\right),\\ i^{\ast }(\tau ,t)& =i^{\ast }\left(\tau \right)+\eta (\tau ,t),\\ N_1\left(t\right)& =N_1^{\ast }+n_1\left(t\right),\\ N_2\left(t\right)& =N_2^{\ast }+n_2\left(t\right).\end{array}$(13)

With this perturbations, the system changes to the following form: (14) $\begin{array}{rl}{s}^{\prime }\left(t\right)& =r_1\left(1-\frac{N_1^{\ast }+n_1\left(t\right)+{\alpha }_{12}(N_2^{\ast }+n_2(t\left)\right)}{K_1}\left)\right(N_1^{\ast }+n_1\left(t\right)\right)\\ & -\mu (S^{\ast }+s(t\left)\right)-(S^{\ast }+s(t\left)\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\left[i^{\ast }\right(\theta )+\eta (\theta ,t\left)\right]\mathrm{d}\theta ,\\ i_{\tau }^{\ast }\left(\tau \right)+{\eta }_{\theta }+{\eta }_t& =-\alpha \left(\tau \right)(i^{\ast }+\eta (\tau ,t\left)\right)-\mu (i^{\ast }+\eta (\tau ,t\left)\right),\\ i^{\ast }\left(0\right)+\eta (0,t)& =(S^{\ast }+s(t\left)\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)(i^{\ast }+\eta (\theta ,t\left)\right)\mathrm{d}\theta ,\\ n_2^{\prime }\left(t\right)& =r_2\left(1-\frac{N_2^{\ast }+n_2+{\alpha }_{21}(N_1^{\ast }+n_1(t\left)\right)}{K_2}\right)\\ & \times (N_2^{\ast }+n_2(t\left)\right)-{\mu }_2(N_2^{\ast }+n_2(t\left)\right).\end{array}$(14)

Linearizing these equations and only retaining the terms which are linear we obtain the following set of equations for the perturbations. (15) $\begin{array}{rl}{s}^{\prime }\left(t\right)& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1^{\ast }+r_1{n}_1\left(1-\frac{N_1^{\ast }+{\alpha }_{12}{N}_2^{\ast }}{K_1}\right)-\mu s\\ & -S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta -s{\int }_0^{\mathrm{\infty }}{\beta }_1{i}^{\ast }\mathrm{d}\theta ,\\ {\eta }_{\theta }+{\eta }_t& =-\alpha \left(\tau \right)\eta (\tau ,t)-\mu \eta (\tau ,t),\\ \eta (0,t)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta (\theta ,t)\mathrm{d}\theta +s\left(t\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ n_2^{\prime }\left(t\right)& =r_2\left(-\frac{n_2+{\alpha }_{21}{n}_1}{K_2}\right)N_2^{\ast }.\end{array}$(15)

5.1. Stability of the trivial equilibrium

Theorem 5.1.

When $R_0<1,$ the trivial equilibrium is locally asymptotically stable. It is unstable when $R_0>1,$ where $R_0=max(R_0^{N_{10}},R_0^{N_2})$, i.e. the equilibrium is locally asymptotically stable when both $R_0^{N_{10}},R_0^{N_2}$ are less than one and unstable when at least one of them is greater than one.

Proof.

This equilibrium is defined by $E_0=(0,0,0)$. This leads to the following set of equations: (16) $\begin{array}{rl}{s}^{\prime }\left(t\right)& =r_1{n}_1\left(t\right)-\mu s\left(t\right),\\ {\eta }_{\tau }+{\eta }_t& =-(\alpha +\mu )\eta ,\\ \eta \left(0\right)& =0,\\ n_2^{\prime }& =(r_2-{\mu }_2)n_2.\end{array}$(16)

The solution of $n_2$ is given by $n_2=n_2\left(0\right){\mathrm{e}}^{(r_2-{\mu }_2)t}$. From the last equation in (16(16) $\begin{array}{rl}{s}^{\prime }\left(t\right)& =r_1{n}_1\left(t\right)-\mu s\left(t\right),\\ {\eta }_{\tau }+{\eta }_t& =-(\alpha +\mu )\eta ,\\ \eta \left(0\right)& =0,\\ n_2^{\prime }& =(r_2-{\mu }_2)n_2.\end{array}$(16) ) it is clear that when $R_0^{N_2}>1\Rightarrow (r_2-{\mu }_2)>0$ and hence the system is unstable. When $R_0^{N_2}<1$, we have to investigate other eigenvalues which come from the other equations in the system given by Equation (16(16) $\begin{array}{rl}{s}^{\prime }\left(t\right)& =r_1{n}_1\left(t\right)-\mu s\left(t\right),\\ {\eta }_{\tau }+{\eta }_t& =-(\alpha +\mu )\eta ,\\ \eta \left(0\right)& =0,\\ n_2^{\prime }& =(r_2-{\mu }_2)n_2.\end{array}$(16) ). Note the boundary condition on η forces the general solution to be of the form $\eta \left(t\right)=0$ for all t. This leads to the fact $s\left(t\right)=n_1\left(t\right)$.

Hence we have the solution of s as $s\left(t\right)=s\left(0\right){\mathrm{e}}^{(r_1-\mu )t}$. Note that $R_0^{N_{10}}>1\Rightarrow (r_1-\mu )>0$ and hence the system is unstable. When $R_0^{N_{10}}<1$ it follows that $(r_1-\mu )<0$ and hence the system is locally asymptotically stable.

Hence, all possible solutions to this system of equations are either negative or have negative real parts if and only if both $R_0^{N_{10}}$ and $R_0^{N_2}$ are less than one. If any one of these is greater than one, the system is unstable.

5.2. Stability of the semi-trivial equilibria

In this section, we explore the local stability of semi-trivial equilibria. There are three semi-trivial equilibria. We begin by deriving the local stability of the boundary equilibria $E_1=(S^{\ast },0,0)$ and $E_2=(0,0,N_2^{\ast })$, and then explore the stability of the other non-trivial boundary equilibrium $E_{N1}$.

5.2.1. Stability of $E_1$

Theorem 5.2.

This equilibrium $E_1$ is stable when $R_0^{N_{10}}>1$ but $R_0^{N_1}<1$ and $R_{N1N2}<1$.

Proof.

Note in this case we have $S^{\ast }=N_1^{\ast }$. Substituting $(S^{\ast },0,0)$ for this equilibrium in the equation of stability in Equation (15(15) $\begin{array}{rl}{s}^{\prime }\left(t\right)& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1^{\ast }+r_1{n}_1\left(1-\frac{N_1^{\ast }+{\alpha }_{12}{N}_2^{\ast }}{K_1}\right)-\mu s\\ & -S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta -s{\int }_0^{\mathrm{\infty }}{\beta }_1{i}^{\ast }\mathrm{d}\theta ,\\ {\eta }_{\theta }+{\eta }_t& =-\alpha \left(\tau \right)\eta (\tau ,t)-\mu \eta (\tau ,t),\\ \eta (0,t)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta (\theta ,t)\mathrm{d}\theta +s\left(t\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ n_2^{\prime }\left(t\right)& =r_2\left(-\frac{n_2+{\alpha }_{21}{n}_1}{K_2}\right)N_2^{\ast }.\end{array}$(15) ), where $S^{\ast }=K_1(1-\mu /r_1)$, we have the following set of equations: (17) $\begin{array}{rl}{s}^{\prime }\left(t\right)& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1^{\ast }+r_1{n}_1\left(1-\frac{N_1^{\ast }}{K_1}\right)-\mu s-N_1^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta ,\\ {\eta }_{\tau }+{\eta }_t& =-\alpha \left(\tau \right)\eta (\tau ,t)-\mu \eta (\tau ,t),\\ \eta (0,t)& =N_1^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta (\theta ,t)\mathrm{d}\theta ,\\ n_2^{\prime }\left(t\right)& =r_2{n}_2\left(1-\frac{{\alpha }_{21}{N}_1^{\ast }}{K_2}\right)-{\mu }_2{n}_2.\end{array}$(17)

The last equation can be solved to obtain a closed-form solution $n_2\left(t\right)=n_2\left(0\right){\mathrm{e}}^{\left(r_2\right(1-{\alpha }_{21}{N}_1^{\ast }/K_2)-{\mu }_2)t}.$ Substituting $N_1^{\ast }=S^{\ast }=K_1(1-\mu /r_1)$, we obtain the solution as, $n_2\left(t\right)=n_2\left(0\right){\mathrm{e}}^{\left(r_2/K_2\right)\left[\right(1-{\mu }_2/r_2)K_2-{\alpha }_{21}{K}_1(1-\mu /r_1\left)\right]t}.$

From the definition of $R_{N1N2}$ it follows that if $R_{N1N2}>1$ this equilibrium is unstable. When $R_{N1N2}<1$, we shall explore other eigenvalues in this system. We linearize the system using the following functions, $s\left(t\right)=s_0{\mathrm{e}}^{\lambda t},\eta (\tau ,t)=\eta \left(\tau \right){\mathrm{e}}^{\lambda t},n_1\left(t\right)=n_1{\mathrm{e}}^{\lambda t}$. This gives the following set of equations: (18) $\begin{array}{rl}\lambda s_0& =-\frac{r_1}{K_1}\left(n_1\right)N_1^{\ast }+r_1{n}_1\left(1-\frac{N_1^{\ast }}{K_1}\right)-\mu s-N_1^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta ,\\ {\eta }_{\tau }& =-(\lambda +\alpha (\tau )+\mu )\eta ,\\ \eta \left(0\right)& =N_1^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta \left(\theta \right)\mathrm{d}\theta .\end{array}$(18)

Solving the equation for η we have $\eta =\eta \left(0\right){\mathrm{e}}^{-{\int }_0^{\tau }(\lambda +\alpha (r)+\mu )\mathrm{d}r}$. Substituting this in the boundary condition, the equation leads to the following characteristic equation $\mathcal{H}\left(\lambda \right)=1,$ where $\mathcal{H}\left(\lambda \right)=N_1^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\alpha (r)+\mu )\mathrm{d}r}\mathrm{d}\theta .$

Note that as $\lambda \to \mathrm{\infty }\Rightarrow \mathcal{H}\left(\lambda \right)\to 0$. We observe that $\mathcal{H}\left(0\right)=K_1\left(1-\frac{\mu }{r_1}\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{-{\int }_0^{\theta }\left(\alpha \right(r)+\mu )\mathrm{d}r}\mathrm{d}\theta .$

When $R_0^{N_1}>1$ it is clear that $\mathcal{H}\left(0\right)>1$. Thus, there exists a positive solution for λ when $\mathcal{H}\left(0\right)>1$ and the system is unstable. For the case, when $R_0^{N_1}<1$, let us assume there exists a $\lambda =a+\mathrm{i}b$ with nonnegative real part $(a\ge 0)$. This shows that $\begin{array}{rl}|\mathcal{H}\left(\lambda \right)|& =|N_1^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\alpha (r)+\mu )\mathrm{d}r}\mathrm{d}\theta |\\ & \le K_1\left(1-\frac{\mu }{r_1}\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{-{\int }_0^{\theta }\left(\alpha \right(r)+\mu )\mathrm{d}r}\mathrm{d}\theta <1,\end{array}$ whenever $R_0^{N_1}<1$, which is a contradiction. Hence, there does not exist any root λ with nonnegative real part when $R_0^{N_1}<1$. All roots are either negative or have negative real part. We will explore another possible eigenvalue in this case. From the first equation in (18(18) $\begin{array}{rl}\lambda s_0& =-\frac{r_1}{K_1}\left(n_1\right)N_1^{\ast }+r_1{n}_1\left(1-\frac{N_1^{\ast }}{K_1}\right)-\mu s-N_1^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta ,\\ {\eta }_{\tau }& =-(\lambda +\alpha (\tau )+\mu )\eta ,\\ \eta \left(0\right)& =N_1^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta \left(\theta \right)\mathrm{d}\theta .\end{array}$(18) ) we have, $\lambda =-2\left(r_1{N}_1^{\ast }/K_1\right)+r_1-\mu$. Substituting the value of $N_1^{\ast }$, we have the following eigenvalue, $\lambda =-r_1(1-\mu /r_1)$. Note this equilibrium is unstable when $R_0^{N10}<1$ and stable when $R_0^{N10}>1$.

5.2.2. Stability of $E_2$

The equilibrium in this part is given by $E_2=(0,0,N_2^{\ast })$ where $N_2^{\ast }=K_2(1-{\mu }_2/r_2)$. This equilibrium exists only if $R_0^{N_2}>1$.

Theorem 5.3.

Assume $R_0^{N_2}>1$. The equilibrium $E_2$ is locally asymptotically stable if and only if $R_{N2}<1$.

Proof.

Substituting $S^{\ast }=0,i^{\ast }=0$ and $N_2^{\ast }$ in the linearized equations for local stability (14(14) $\begin{array}{rl}{s}^{\prime }\left(t\right)& =r_1\left(1-\frac{N_1^{\ast }+n_1\left(t\right)+{\alpha }_{12}(N_2^{\ast }+n_2(t\left)\right)}{K_1}\left)\right(N_1^{\ast }+n_1\left(t\right)\right)\\ & -\mu (S^{\ast }+s(t\left)\right)-(S^{\ast }+s(t\left)\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\left[i^{\ast }\right(\theta )+\eta (\theta ,t\left)\right]\mathrm{d}\theta ,\\ i_{\tau }^{\ast }\left(\tau \right)+{\eta }_{\theta }+{\eta }_t& =-\alpha \left(\tau \right)(i^{\ast }+\eta (\tau ,t\left)\right)-\mu (i^{\ast }+\eta (\tau ,t\left)\right),\\ i^{\ast }\left(0\right)+\eta (0,t)& =(S^{\ast }+s(t\left)\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)(i^{\ast }+\eta (\theta ,t\left)\right)\mathrm{d}\theta ,\\ n_2^{\prime }\left(t\right)& =r_2\left(1-\frac{N_2^{\ast }+n_2+{\alpha }_{21}(N_1^{\ast }+n_1(t\left)\right)}{K_2}\right)\\ & \times (N_2^{\ast }+n_2(t\left)\right)-{\mu }_2(N_2^{\ast }+n_2(t\left)\right).\end{array}$(14) ), we obtain the following set of equations: (19) $\begin{array}{rl}{s}^{\prime }& =r_1{n}_1\left(t\right)\left[1-\frac{{\alpha }_{12}{N}_2^{\ast }}{K_1}\right]-\mu s\left(t\right),\\ {\eta }_t+{\eta }_{\tau }& =-(\alpha +\mu )\eta (\theta ,t),\\ \eta (0,t)& =0,\\ n_2^{\prime }& =\frac{r_2}{K_2}{N}_2^{\ast }(n_2+{\alpha }_{21}{n}_1)+r_2{n}_2\left(1-\frac{N_2^{\ast }}{K_2}\right)-{\mu }_2{n}_2.\end{array}$(19)

We look for separable solutions of Equations (19(19) $\begin{array}{rl}{s}^{\prime }& =r_1{n}_1\left(t\right)\left[1-\frac{{\alpha }_{12}{N}_2^{\ast }}{K_1}\right]-\mu s\left(t\right),\\ {\eta }_t+{\eta }_{\tau }& =-(\alpha +\mu )\eta (\theta ,t),\\ \eta (0,t)& =0,\\ n_2^{\prime }& =\frac{r_2}{K_2}{N}_2^{\ast }(n_2+{\alpha }_{21}{n}_1)+r_2{n}_2\left(1-\frac{N_2^{\ast }}{K_2}\right)-{\mu }_2{n}_2.\end{array}$(19) ) in the form $s=s_0{\mathrm{e}}^{\lambda t},\eta (\theta ,t)=\eta \left(\theta \right){\mathrm{e}}^{\lambda t},n_2\left(t\right)=n_2{\mathrm{e}}^{\lambda t}$. Substituting these forms in the system of Equations (19(19) $\begin{array}{rl}{s}^{\prime }& =r_1{n}_1\left(t\right)\left[1-\frac{{\alpha }_{12}{N}_2^{\ast }}{K_1}\right]-\mu s\left(t\right),\\ {\eta }_t+{\eta }_{\tau }& =-(\alpha +\mu )\eta (\theta ,t),\\ \eta (0,t)& =0,\\ n_2^{\prime }& =\frac{r_2}{K_2}{N}_2^{\ast }(n_2+{\alpha }_{21}{n}_1)+r_2{n}_2\left(1-\frac{N_2^{\ast }}{K_2}\right)-{\mu }_2{n}_2.\end{array}$(19) ) and simplifying the system we obtain the following linear eigenvalue problem: (20) $\begin{array}{rl}{\eta }_{\tau }& =-(\lambda +\alpha (\tau )+\mu )\eta ,\\ \eta \left(0\right)& =0,\\ n_1& =s_0+{\int }_0^{\mathrm{\infty }}\eta \left(\tau \right)\mathrm{d}\tau ,\\ \lambda n_2& =-\frac{r_2}{K_2}{N}_2^{\ast }(n_2+n_1)+r_2{n}_2-{\mu }_2{n}_2-\frac{r_2{N}_2^{\ast }}{K_2}{n}_2.\end{array}$(20)

It is easy to observe from the boundary condition of $\eta \left(\tau \right)$ that $\eta \left(\tau \right)=0$. The eigenvalues will follow from solving the other equations. From $\eta \left(\tau \right)=0$ it follows that $n_1=s_0$. The first equation in (20(20) $\begin{array}{rl}{\eta }_{\tau }& =-(\lambda +\alpha (\tau )+\mu )\eta ,\\ \eta \left(0\right)& =0,\\ n_1& =s_0+{\int }_0^{\mathrm{\infty }}\eta \left(\tau \right)\mathrm{d}\tau ,\\ \lambda n_2& =-\frac{r_2}{K_2}{N}_2^{\ast }(n_2+n_1)+r_2{n}_2-{\mu }_2{n}_2-\frac{r_2{N}_2^{\ast }}{K_2}{n}_2.\end{array}$(20) ) changes to the following form: (21) $\lambda n_1=r_1{n}_1\left[1-\frac{{\alpha }_{12}{N}_2^{\ast }}{K_1}\right]-\mu n_1.$(21)

We will determine the eigenvalues for the following two possibilities of $n_1$.

Case 1 Let $n_1\ne 0$. This will reduce Equation (21(21) $\lambda n_1=r_1{n}_1\left[1-\frac{{\alpha }_{12}{N}_2^{\ast }}{K_1}\right]-\mu n_1.$(21) ) to $\lambda =r_1-\left[\frac{r_1{\alpha }_{12}{K}_2}{K_1}\left(1-\frac{{\mu }_2}{r_2}\right)+\mu \right].$

When the reproduction number for species two $R_{N2}<1$, we have $\lambda <0$. In this case we continue with Case 2. On the other hand, when $R_{N2}>1$ the eigenvalue λ will be positive and as a consequence the equilibrium will be unstable. Note the reproduction number as defined before is given as $R_{N2}=\frac{r_1}{\frac{r_1{\alpha }_{12}{K}_2}{K_1}\left(1-\frac{{\mu }_2}{r_2}\right)+\mu }.$

Case 2 Let $n_1=0$. In this case, the remaining eigenvalue will be given by the equation involving $n_2$ as (22) $\lambda =-\frac{r_2}{K_2}{N}_2^{\ast }+r_2-{\mu }_2-\frac{r_2{N}_2^{\ast }}{K_2}.$(22)

Substituting $N_2^{\ast }=K_2(r_2-{\mu }_2)/r_2$ this reduces to $\lambda =-r_2{N}_2^{\ast }/K_2$ which is a negative quantity. Hence, if $R_{N2}<1$, this equilibrium is locally asymptotically stable. This completes the proof.

5.2.3. Stability of $E_{N1}$

The endemic equilibrium $E_{N1}$ is given by $E_{N1}=(S^{\ast },i^{\ast }(\tau ),0)$, where $S^{\ast }=\frac{1}{{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\mathrm{\Pi }\left(\theta \right)\mathrm{d}\theta }\mathrm{and}{i}^{\ast }\left(\tau \right)=i^{\ast }\left(0\right)\mathrm{\Pi }\left(\tau \right).$

The next theorem gives the stability of $E_{N1}$. To state the theorem, define $\sigma =r_1\left(1-\frac{2N_1^{\ast }}{K_1}\right).$

Theorem 5.4.

Assume $E_{N1}$ exists. If $R_{N1}>1,$ then $E_{N1}$ is unstable. If $R_{N1}<1,$ then

• if $\sigma >0$ and $\mu -\sigma >0,$ then $E_{N1}$ is locally asymptotically stable;

• if $\sigma >0$ and $\mu -\sigma <0,$ then $E_{N1}$ is unstable.

Remark 5.1.

In the case $\sigma <0$, we could not prove that the system is stable. In addition, it is not hard to see that the linearized system does have a positive eigenvalue in this case. Thus, a possibility exists that a pair of complex conjugate eigenvalues cross the imaginary axis causing Hopf bifurcation to occur and the presence of sustained oscillations.

Proof.

Using the linearization defined in Equation (14(14) $\begin{array}{rl}{s}^{\prime }\left(t\right)& =r_1\left(1-\frac{N_1^{\ast }+n_1\left(t\right)+{\alpha }_{12}(N_2^{\ast }+n_2(t\left)\right)}{K_1}\left)\right(N_1^{\ast }+n_1\left(t\right)\right)\\ & -\mu (S^{\ast }+s(t\left)\right)-(S^{\ast }+s(t\left)\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\left[i^{\ast }\right(\theta )+\eta (\theta ,t\left)\right]\mathrm{d}\theta ,\\ i_{\tau }^{\ast }\left(\tau \right)+{\eta }_{\theta }+{\eta }_t& =-\alpha \left(\tau \right)(i^{\ast }+\eta (\tau ,t\left)\right)-\mu (i^{\ast }+\eta (\tau ,t\left)\right),\\ i^{\ast }\left(0\right)+\eta (0,t)& =(S^{\ast }+s(t\left)\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)(i^{\ast }+\eta (\theta ,t\left)\right)\mathrm{d}\theta ,\\ n_2^{\prime }\left(t\right)& =r_2\left(1-\frac{N_2^{\ast }+n_2+{\alpha }_{21}(N_1^{\ast }+n_1(t\left)\right)}{K_2}\right)\\ & \times (N_2^{\ast }+n_2(t\left)\right)-{\mu }_2(N_2^{\ast }+n_2(t\left)\right).\end{array}$(14) ) and using the values of this equilibrium, we arrive at the following set of equations: (23) $\begin{array}{rl}{s}^{\prime }& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1^{\ast }+r_1{n}_1\left[1-\frac{N_1^{\ast }}{K_1}\right]-\mu s\\ & -S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta (\theta ,t)\mathrm{d}\theta -s{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ {\eta }_{\theta }+{\eta }_t& =-\left(\alpha \right(\theta )+\mu )\eta (\theta ,t),\\ \eta (0,t)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta (\theta ,t)\mathrm{d}\theta +s\left(t\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ n_1\left(t\right)& =s\left(t\right)+{\int }_0^{\mathrm{\infty }}\eta (\theta ,t)\mathrm{d}\theta ,\\ n_2^{\prime }& =\left[r_2\left(1-\frac{{\alpha }_{21}{N}_1^{\ast }}{K_2}\right)-{\mu }_2\right]n_2.\end{array}$(23)

Note that the equation for $n_2$ separates from the remaining equations and we have $n_2=n_2\left(0\right){\mathrm{e}}^{\left(r_2\right(1-{\alpha }_{21}{N}_1^{\ast }/K_2)-{\mu }_2)t}.$ When $R_{N1}>1$, the exponent is positive and hence the system is unstable. To determine the stability of this system when $R_{N1}<1$, we have to investigate the remaining eigenvalues.

We are considering roots of the characteristic equation which are different from the previous one. We look for separable solutions in the form, $s\left(t\right)=s_0{\mathrm{e}}^{\lambda t},\eta ={\mathrm{e}}^{\lambda t}\eta \left(\tau \right),n_1\left(t\right)=n_1{\mathrm{e}}^{\lambda t}$. We solve for η to obtain the solution in the form $\eta \left(\tau \right)=\eta \left(0\right){\mathrm{e}}^{-{\int }_0^{\tau }(\lambda +\mu +\alpha (s\left)\right)\mathrm{d}s}$. Using this set of equations and substituting this in Equation (23(23) $\begin{array}{rl}{s}^{\prime }& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1^{\ast }+r_1{n}_1\left[1-\frac{N_1^{\ast }}{K_1}\right]-\mu s\\ & -S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta (\theta ,t)\mathrm{d}\theta -s{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ {\eta }_{\theta }+{\eta }_t& =-\left(\alpha \right(\theta )+\mu )\eta (\theta ,t),\\ \eta (0,t)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta (\theta ,t)\mathrm{d}\theta +s\left(t\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ n_1\left(t\right)& =s\left(t\right)+{\int }_0^{\mathrm{\infty }}\eta (\theta ,t)\mathrm{d}\theta ,\\ n_2^{\prime }& =\left[r_2\left(1-\frac{{\alpha }_{21}{N}_1^{\ast }}{K_2}\right)-{\mu }_2\right]n_2.\end{array}$(23) ), we have (24) $\begin{array}{rl}\lambda s_0& =-\frac{r_1}{K_1}{n}_1{N}_1^{\ast }+r_1{n}_1\left[1-\frac{N_1^{\ast }}{K_1}\right]-\mu s_0-\eta \left(0\right),\\ \eta \left(0\right)& =S^{\ast }\eta \left(0\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\alpha (s)+\mu )\mathrm{d}s}\mathrm{d}\theta +s_0{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ n_1& =s_0+{\int }_0^{\mathrm{\infty }}\eta \left(\theta \right)\mathrm{d}\theta .\end{array}$(24)

Solving this system for $s_0$ and $n_1$ from the first equation in Equation (24(24) $\begin{array}{rl}\lambda s_0& =-\frac{r_1}{K_1}{n}_1{N}_1^{\ast }+r_1{n}_1\left[1-\frac{N_1^{\ast }}{K_1}\right]-\mu s_0-\eta \left(0\right),\\ \eta \left(0\right)& =S^{\ast }\eta \left(0\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\alpha (s)+\mu )\mathrm{d}s}\mathrm{d}\theta +s_0{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ n_1& =s_0+{\int }_0^{\mathrm{\infty }}\eta \left(\theta \right)\mathrm{d}\theta .\end{array}$(24) ), we obtain (25) $\begin{array}{rl}{s}_0(\lambda +\mu )& =\sigma n_1-\eta \left(0\right),\\ n_1& =s_0+\eta \left(0\right){\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\mu +\alpha (r\left)\right)\mathrm{d}r}\mathrm{d}\theta .\end{array}$(25)

Combining the two equations from Equation (25(25) $\begin{array}{rl}{s}_0(\lambda +\mu )& =\sigma n_1-\eta \left(0\right),\\ n_1& =s_0+\eta \left(0\right){\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\mu +\alpha (r\left)\right)\mathrm{d}r}\mathrm{d}\theta .\end{array}$(25) ), we obtain the following solution for $s_0$: $s_0=\frac{\eta \left(0\right)}{\lambda +\mu -\sigma }\left[(\sigma {\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\mu +\alpha (r\left)\right)\mathrm{d}r}\mathrm{d}\theta -1\right].$

Substituting this equation in the second equation in Equation (24(24) $\begin{array}{rl}\lambda s_0& =-\frac{r_1}{K_1}{n}_1{N}_1^{\ast }+r_1{n}_1\left[1-\frac{N_1^{\ast }}{K_1}\right]-\mu s_0-\eta \left(0\right),\\ \eta \left(0\right)& =S^{\ast }\eta \left(0\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\alpha (s)+\mu )\mathrm{d}s}\mathrm{d}\theta +s_0{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ n_1& =s_0+{\int }_0^{\mathrm{\infty }}\eta \left(\theta \right)\mathrm{d}\theta .\end{array}$(24) ) we obtain the characteristic equation $\mathcal{G}\left(\lambda \right)=1$ where $\mathcal{G}\left(\lambda \right)=S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\mu +\alpha (s\left)\right)\mathrm{d}s}\mathrm{d}\theta +B\frac{-1+\sigma {\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\mu +\alpha )\mathrm{d}s}\mathrm{d}\theta }{\lambda +\mu -\sigma }.$ and B denotes the constant $B={\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta$. The eigenvalues of the problem (24(24) $\begin{array}{rl}\lambda s_0& =-\frac{r_1}{K_1}{n}_1{N}_1^{\ast }+r_1{n}_1\left[1-\frac{N_1^{\ast }}{K_1}\right]-\mu s_0-\eta \left(0\right),\\ \eta \left(0\right)& =S^{\ast }\eta \left(0\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\alpha (s)+\mu )\mathrm{d}s}\mathrm{d}\theta +s_0{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ n_1& =s_0+{\int }_0^{\mathrm{\infty }}\eta \left(\theta \right)\mathrm{d}\theta .\end{array}$(24) ) are the values of λ that solve the characteristic equation. Let λ be the solution for the characteristic equation $\mathcal{G}\left(\lambda \right)=1$.

Using the following notations, we define the respective functions. $\begin{array}{rl}{R}_1\left(\lambda \right)& ={\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\mu +\alpha (s\left)\right)\mathrm{d}s}\mathrm{d}\theta ,\\ \rho \left(\lambda \right)& ={\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\mu +\alpha )\mathrm{d}s}\mathrm{d}\theta ,\sigma =r_1\left(1-\frac{2N_1}{K_1}\right).\end{array}$

The notations reduce the characteristic equation to the form (26) $1+B\frac{1-\sigma \rho \left(\lambda \right)}{\lambda +\mu -\sigma }=S^{\ast }{R}_1\left(\lambda \right).$(26)

Let $\lambda =a+bi$ with $a\ge 0$. For these λ, the right-hand side of this equation satisfies $|S^{\ast }{R}_1\left(\lambda \right)|\le S^{\ast }{R}_1\left(0\right)\le 1.$ The left-hand side gives (27) $\frac{|\lambda +\mu -\sigma +B-B\sigma \rho \left(\lambda \right)|}{|\lambda +\mu -\sigma |}=\frac{\sqrt{(a+\mu -\sigma +B-B\sigma \mathrm{\Re }\rho (\lambda ){)}^2+(b+B\sigma \mathrm{\Im }\rho \left(\lambda \right){)}^2}}{\sqrt{(a+\mu -\sigma {)}^2+b^2}}>1.$(27) The last inequality holds, because $\mathrm{\Im }\rho \left(\lambda \right)>0$. In addition, $B-B\sigma \rho \left(0\right)>B-B\frac{\sigma }{\mu }\ge 0.$

In the second case when $\mu -\sigma <0$, it is not hard to see that the characteristic equation has a positive root.

5.3. Stability of the co-existence equilibria

In this section, we shall explore the stability of the co-existence equilibria of the two species both when the disease is present in species one and when it is absent.

5.3.1. Stability of the equilibrium $E_0^{\ast }=(N_1^{\ast },0,N_2^{\ast })$

In this case the equilibrium values $N_1^{\ast }$ and $N_2^{\ast }$ satisfy the following set of equations: $\begin{array}{rl}{N}_1^{\ast }+{\alpha }_{12}{N}_2^{\ast }& =K_1\left(1-\frac{\mu }{r_1}\right),\\ {\alpha }_{21}{N}_1^{\ast }+N_2^{\ast }& =K_2\left(1-\frac{{\mu }_2}{r_2}\right).\end{array}$

Theorem 5.5.

Assume $R_{N2}>1$ and $R_{N1N2}>1$. Assume also $R_0^{N_1}<1$. Then, the equilibrium $E_0^{\ast }$ is locally asymptotically stable.

Proof.

The set of equations for the equilibrium values reduce the stability equation to the following form: (28) $\begin{array}{rl}{s}^{\prime }\left(t\right)& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1^{\ast }-N_1^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta +\mu n_1-\mu s,\\ {\eta }_{\tau }+{\eta }_t& =-\left(\alpha \right(\tau )+\mu )\eta ,\\ \eta (0,t)& =N_1^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta (\theta ,t)\mathrm{d}\theta ,\\ n_2^{\prime }& =-r_2\left(\frac{n_2+{\alpha }_{21}{n}_1}{K_2}\right)N_2^{\ast },\\ n_1& =s+{\int }_0^{\mathrm{\infty }}\eta (\theta ,t)\mathrm{d}\theta .\end{array}$(28)

Using the substitution $\eta (\theta ,t)=\overline{\eta }\left(\theta \right){\mathrm{e}}^{\lambda t}$, it reduces the equation for η to ${\overline{\eta }}^{\prime }=-(\lambda +\mu +\alpha (\theta \left)\right)\overline{\eta },$ which can be integrated to obtain (29) $\overline{\eta }\left(\theta \right)=\overline{\eta }\left(0\right){\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\mu +\alpha (s\left)\right)\mathrm{d}s}.$(29)

Using the boundary condition, we obtain the characteristic equation as $G\left(\lambda \right)=1$, where $G\left(\lambda \right)=N_1^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{-{\int }_0^{\theta }(\lambda +\mu +\alpha (s\left)\right)\mathrm{d}s}\mathrm{d}\theta .$

If there exists a solution λ of the characteristic equation, such that $\lambda =a+ib,a\ge 0$, we have $|\mathrm{LHS}|=1$ whereas $\begin{array}{rl}|\mathrm{RHS}|& =|N_1^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{{\int }_0^{\theta }(\lambda +\mu +\alpha (s\left)\right)\mathrm{d}s}\mathrm{d}\theta |\\ & <K_1\left(1-\frac{\mu }{r_1}\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right){\mathrm{e}}^{-{\int }_0^{\theta }(\mu +\alpha (s\left)\right)\mathrm{d}s}\mathrm{d}\theta <1\end{array}$ which is a contradiction. This last inequality follows from the fact that $R_0^{N_1}<1$. Hence, all roots are negative or have negative real parts. Now, we explore the remaining roots. We assume $\eta \equiv 0$. Note that this results in $n_1=s$. The roots are obtained from the following set of equations: (30) $\begin{array}{rl}{n}_1^{\prime }& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1^{\ast },\\ n_2^{\prime }& =-\frac{r_2}{K_2}(n_2+{\alpha }_{21}{n}_1)N_2^{\ast }.\end{array}$(30)

Using the solution in the form, $n_1\left(t\right)={\mathrm{e}}^{\lambda t}{n}_1,n_2\left(t\right)={\mathrm{e}}^{\lambda t}{n}_2$, we have (31) $\begin{array}{rl}\left(\lambda +\frac{r_1{N}_1^{\ast }}{K_1}\right)n_1& =-\frac{r_1}{K_1}{\alpha }_{12}{N}_1^{\ast }{n}_2,\\ \left(\lambda +\frac{r_2{N}_2^{\ast }}{K_2}\right)n_2& =-\frac{r_2}{K_2}{\alpha }_{21}{N}_2^{\ast }{n}_1.\end{array}$(31)

This can be combined to obtain the characteristic equation $H\left(\lambda \right)=0$, where $H\left(\lambda \right)$ is the polynomial function given as $H\left(\lambda \right)={\lambda }^2+\lambda \left(\frac{r_1{N}_1^{\ast }}{K_1}+\frac{r_2{N}_2^{\ast }}{K_2}\right)+\frac{r_1{r}_2}{K_1{K}_2}{N}_1^{\ast }{N}_2^{\ast }(1-{\alpha }_{12}{\alpha }_{21}).$

Since we have $1-{\alpha }_{12}{\alpha }_{21}>0$ in this case, all roots λ are either negative or have negative real parts. Hence, this equilibrium is always stable whenever $R_0^{N_1}<1$.

5.3.2. Stability of the equilibria $E^{\ast }=(N_1^{\ast },i^{\ast }(\tau ),N_2^{\ast })$

Regarding the stability of the disease-present coexistence equilibrium, we have the following result:

Theorem 5.6.

Let $\alpha \left(\tau \right)=0$. If $1-{\alpha }_{12}{\alpha }_{21}>0,$ then the equilibrium $E^{\ast }$ is locally asymptotically stable whenever it exists.

Proof.

From the linearization presented in Equation (15(15) $\begin{array}{rl}{s}^{\prime }\left(t\right)& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1^{\ast }+r_1{n}_1\left(1-\frac{N_1^{\ast }+{\alpha }_{12}{N}_2^{\ast }}{K_1}\right)-\mu s\\ & -S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta -s{\int }_0^{\mathrm{\infty }}{\beta }_1{i}^{\ast }\mathrm{d}\theta ,\\ {\eta }_{\theta }+{\eta }_t& =-\alpha \left(\tau \right)\eta (\tau ,t)-\mu \eta (\tau ,t),\\ \eta (0,t)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta (\theta ,t)\mathrm{d}\theta +s\left(t\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ n_2^{\prime }\left(t\right)& =r_2\left(-\frac{n_2+{\alpha }_{21}{n}_1}{K_2}\right)N_2^{\ast }.\end{array}$(15) ), we have the following sets of equations: (32) $\begin{array}{rl}{s}^{\prime }\left(t\right)& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1^{\ast }+r_1{n}_1\left(1-\frac{N_1^{\ast }+{\alpha }_{12}{N}_2^{\ast }}{K_1}\right)\\ & -\mu s-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta -s{\int }_0^{\mathrm{\infty }}{\beta }_1{i}^{\ast }\mathrm{d}\theta ,\\ {\eta }_{\theta }+{\eta }_t& =-\alpha \left(\tau \right)\eta (\tau ,t)-\mu \eta (\tau ,t),\\ \eta (0,t)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta (\theta ,t)\mathrm{d}\theta +s\left(t\right){\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ n_2^{\prime }\left(t\right)& =r_2\left(-\frac{n_2+{\alpha }_{21}{n}_1}{K_2}\right)N_2^{\ast }+r_2\left(1-\frac{N_2^{\ast }+{\alpha }_{21}{N}_1^{\ast }}{K_2}\right)n_2-{\mu }_2{n}_2.\end{array}$(32)

We look for solutions in the form $s\left(t\right)=s{\mathrm{e}}^{\lambda t},\eta (\tau ,t)={\mathrm{e}}^{\lambda t}\eta \left(\theta \right),n_1\left(t\right)=n_1{\mathrm{e}}^{\lambda t},n_2\left(t\right)=n_2{\mathrm{e}}^{\lambda t}$. This transforms the set of equations to the following form: (33) $\begin{array}{rl}\lambda s& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1+r_1{n}_1\left(1-\frac{N_1+{\alpha }_{12}{N}_2}{K_1}\right)\\ & -\mu s-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta \left(\theta \right)\mathrm{d}\theta -s{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ \eta \left(\tau \right)& =\eta \left(0\right){\mathrm{e}}^{-{\int }_0^{\tau }(\lambda +\mu +\alpha (s\left)\right)\mathrm{d}s},\\ \eta \left(0\right)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta +s{\int }_0^{\mathrm{\infty }}{\beta }_1{i}^{\ast }\mathrm{d}\theta ,\\ \lambda n_2& =-\frac{r_2}{K_2}(n_2+{\alpha }_{21}{n}_1)N_2+A_{21}{n}_2-{\mu }_2{n}_2,\end{array}$(33) where $A_{21}=r_2(1-(N_2^{\ast }+{\alpha }_{21}{N}_1^{\ast })/K_2)$. We further note that from the equilibrium equations, we have $A_{21}={\mu }_2$. Solving the last equation in Equation (33(33) $\begin{array}{rl}\lambda s& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1+r_1{n}_1\left(1-\frac{N_1+{\alpha }_{12}{N}_2}{K_1}\right)\\ & -\mu s-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta \left(\theta \right)\mathrm{d}\theta -s{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ \eta \left(\tau \right)& =\eta \left(0\right){\mathrm{e}}^{-{\int }_0^{\tau }(\lambda +\mu +\alpha (s\left)\right)\mathrm{d}s},\\ \eta \left(0\right)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta +s{\int }_0^{\mathrm{\infty }}{\beta }_1{i}^{\ast }\mathrm{d}\theta ,\\ \lambda n_2& =-\frac{r_2}{K_2}(n_2+{\alpha }_{21}{n}_1)N_2+A_{21}{n}_2-{\mu }_2{n}_2,\end{array}$(33) ), we obtain $n_2$ in terms of $n_1$ as follows: $n_2=\frac{-\left(\frac{r_2}{K_2}{N}_2{\alpha }_{21}\right)n_1}{\lambda +\frac{r_2}{K_2}{N}_2+{\mu }_2-A_{21}},$

Substituting this $n_2$ in the first equation of (33(33) $\begin{array}{rl}\lambda s& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1+r_1{n}_1\left(1-\frac{N_1+{\alpha }_{12}{N}_2}{K_1}\right)\\ & -\mu s-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta \left(\theta \right)\mathrm{d}\theta -s{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ \eta \left(\tau \right)& =\eta \left(0\right){\mathrm{e}}^{-{\int }_0^{\tau }(\lambda +\mu +\alpha (s\left)\right)\mathrm{d}s},\\ \eta \left(0\right)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta +s{\int }_0^{\mathrm{\infty }}{\beta }_1{i}^{\ast }\mathrm{d}\theta ,\\ \lambda n_2& =-\frac{r_2}{K_2}(n_2+{\alpha }_{21}{n}_1)N_2+A_{21}{n}_2-{\mu }_2{n}_2,\end{array}$(33) ), we obtain an explicit form of s in terms of $n_1$ and η as follows: (34) $s=\frac{1}{\lambda +\mu +B}\left[n_1\left\{r_1-\frac{2r_1{N}_1}{K_1}+\frac{r_1{N}_2}{K_1}\left(\frac{{\alpha }_{21}\frac{r_2}{K_2}{N}_1}{\lambda +\frac{r_2}{K_2}{N}_2}-1\right){\alpha }_{12}\right\}-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta \right].$(34)

We replace s using the following form $s=n_1-{\int }_0^{\mathrm{\infty }}\eta \mathrm{d}\theta$. The boundary condition of the PDE in Equation (33(33) $\begin{array}{rl}\lambda s& =-\frac{r_1}{K_1}(n_1+{\alpha }_{12}{n}_2)N_1+r_1{n}_1\left(1-\frac{N_1+{\alpha }_{12}{N}_2}{K_1}\right)\\ & -\mu s-S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)\eta \left(\theta \right)\mathrm{d}\theta -s{\int }_0^{\mathrm{\infty }}{\beta }_1\left(\theta \right)i^{\ast }\left(\theta \right)\mathrm{d}\theta ,\\ \eta \left(\tau \right)& =\eta \left(0\right){\mathrm{e}}^{-{\int }_0^{\tau }(\lambda +\mu +\alpha (s\left)\right)\mathrm{d}s},\\ \eta \left(0\right)& =S^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_1\eta \mathrm{d}\theta +s{\int }_0^{\mathrm{\infty }}{\beta }_1{i}^{\ast }\mathrm{d}\theta ,\\ \lambda n_2& =-\frac{r_2}{K_2}(n_2+{\alpha }_{21}{n}_1)N_2+A_{21}{n}_2-{\mu }_2{n}_2,\end{array}$(33) ) reduces to $\eta \left(0\right)[1+B\rho (\lambda )-S^{\ast }R(\lambda \left)\right]=B{n}_1$, where $\rho \left(\lambda \right)={\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{-{\int }_0^{\tau }(\lambda +\mu +\alpha (s\left)\right)\mathrm{d}s}\mathrm{d}\tau$ and $R\left(\lambda \right)={\int }_0^{\mathrm{\infty }}{\beta }_1\left(\tau \right){\mathrm{e}}^{-{\int }_0^{\tau }(\lambda +\mu +\alpha (s\left)\right)\mathrm{d}s}\mathrm{d}\tau$. Using these substitutions in Equation (35(35) $\begin{array}{rl}{n}_1-\eta \left(0\right)\rho \left(\lambda \right)& =\frac{n_1}{\lambda +\mu +B}\left[\left\{\left(r_1-\frac{2r_1{N}_1}{K_1}\right)+\frac{r_1{N}_2{\alpha }_{12}}{K_1}\left(\frac{{\alpha }_{21}\frac{r_2}{K_2}{N}_1}{\lambda +\frac{r_2}{K_2}{N}_2}-1\right)\right\}\right]\\ & -\frac{\eta \left(0\right)}{\lambda +\mu +B}{S}^{\ast }R\left(\lambda \right).\end{array}$(35) ), we obtain the following form: (35) $\begin{array}{rl}{n}_1-\eta \left(0\right)\rho \left(\lambda \right)& =\frac{n_1}{\lambda +\mu +B}\left[\left\{\left(r_1-\frac{2r_1{N}_1}{K_1}\right)+\frac{r_1{N}_2{\alpha }_{12}}{K_1}\left(\frac{{\alpha }_{21}\frac{r_2}{K_2}{N}_1}{\lambda +\frac{r_2}{K_2}{N}_2}-1\right)\right\}\right]\\ & -\frac{\eta \left(0\right)}{\lambda +\mu +B}{S}^{\ast }R\left(\lambda \right).\end{array}$(35)

Combining all these equations, we obtain the characteristic equation $G\left(\lambda \right)=0$, where $G\left(\lambda \right)$ is given by (36) $G\left(\lambda \right)=[1+B\rho (\lambda )-S^{\ast }R(\lambda \left)\right]\left(1-\frac{J\left(\lambda \right)}{\lambda +\mu +B}\right)-B\rho \left(\lambda \right)+\frac{B{S}^{\ast }R\left(\lambda \right)}{\lambda +\mu +B}.$(36) where $J\left(\lambda \right)=\left(r_1-\frac{2r_1{N}_1}{K_1}\right)+\frac{r_1{N}_2{\alpha }_{12}}{K_1}\left(\frac{{\alpha }_{21}\frac{r_2}{K_2}{N}_1}{\lambda +\frac{r_2}{K_2}{N}_2}-1\right).$

Note that we can rearrange $G\left(\lambda \right)$ to write it in the form, $G\left(\lambda \right)=(1-S^{\ast }R(\lambda \left)\right)\left(1-\frac{J\left(\lambda \right)}{\lambda +\mu +B}\right)-\frac{B\rho \left(\lambda \right)J\left(\lambda \right)}{\lambda +\mu +B}+\frac{B{S}^{\ast }R\left(\lambda \right)}{\lambda +\mu +B}.$

It is clear that $\underset{\lambda \to \mathrm{\infty }}{lim}G\left(\lambda \right)=1$ and $G\left(0\right)=\left(B/(\mu +B)\right)(1-\rho (0\left)J\right(0\left)\right)$. Note that $J\left(0\right)=\frac{r_1{N}_1}{K_1}({\alpha }_{12}{\alpha }_{21}-1)+r_1\left(1-\frac{N_1+{\alpha }_{12}{N}_2}{K_1}\right).$

Hence, we have (37) $1-\rho \left(0\right)J\left(0\right)=1-\rho \left(0\right)r_1\left(1-\frac{N_1+{\alpha }_{12}{N}_2}{K_1}\right)+\frac{\rho \left(0\right)r_1{N}_1}{K_1}(1-{\alpha }_{12}{\alpha }_{21}),$(37) where we recall that $1>{\alpha }_{12}{\alpha }_{21}$, $\rho \left(0\right)<1/\mu$, and that $r_1(1-(N_1+{\alpha }_{12}{N}_2)/K_1)>\mu$. Consequently, the equilibrium is unstable if $G\left(0\right)<0$, or alternatively if $\rho \left(0\right)J\left(0\right)>1$.

Now, we consider the case $\alpha \left(\tau \right)=0$. Then $\rho \left(\lambda \right)=1/(\lambda +\mu )$, $\rho \left(0\right)=1/\mu$ and $r_1(1-(N_1+{\alpha }_{12}{N}_2)/K_1)=\mu$. Hence, $1-\rho \left(0\right)J\left(0\right)=\frac{r_1{N}_1}{\mu K_1}(1-{\alpha }_{12}{\alpha }_{21})>0.$ To show that in this case the characteristic equation does not have roots with nonnegative real parts, we write $J\left(\lambda \right)=\tilde{\sigma }+\frac{E}{\lambda +\frac{r_2{N}_2}{K_2}}$ where $\tilde{\sigma }=r_1\left(1-2N_1/K_1\right)-r_1{\alpha }_{12}{N}_2/K_1$ and $E=r_1{r}_2{\alpha }_{12}{\alpha }_{21}{N}_1{N}_2/K_1{K}_2$. Since the characteristic equation is given by $G\left(\lambda \right)=0$, one can take a common denominator of $\lambda +\mu +B$ to obtain $(1-S^{\ast }R(\lambda \left)\right)(\lambda +\mu +B-J(\lambda \left)\right)-B\rho \left(\lambda \right)J\left(\lambda \right)+B{S}^{\ast }R\left(\lambda \right)=0.$ Further simplification leads to $(1-S^{\ast }R(\lambda \left)\right)(\lambda +\mu -J(\lambda \left)\right)+B-B\rho \left(\lambda \right)J\left(\lambda \right)=0.$ Using the expression for $J\left(\lambda \right)$ above, we have $\begin{array}{rl}& (1-S^{\ast }R(\lambda \left)\right)\left((\lambda +\mu -\tilde{\sigma })\left(\lambda +\frac{r_2{N}_2}{K_2}\right)-E\right)\\ & +\frac{B}{\lambda +\mu }\left((\lambda +\mu -\tilde{\sigma })\left(\lambda +\frac{r_2{N}_2}{K_2}\right)-E\right)=0.\end{array}$ Hence, the characteristic equation simplifies to $\left(1-S^{\ast }R\left(\lambda \right)+\frac{B}{\lambda +\mu }\right)\left((\lambda +\mu -\tilde{\sigma })\left(\lambda +\frac{r_2{N}_2}{K_2}\right)-E\right)=0.$ Thus, in the case $\alpha \left(\tau \right)=0$ the characteristic equation splits into two equations: $1-S^{\ast }R\left(\lambda \right)+\frac{B}{\lambda +\mu }=0(\lambda +\mu -\tilde{\sigma })\left(\lambda +\frac{r_2{N}_2}{K_2}\right)-E=0.$ If we rewrite the first equation as $\frac{\lambda +\mu +B}{\lambda +\mu }=S^{\ast }R\left(\lambda \right)$ then, for λ with $\mathrm{\Re }\lambda \ge 0$ we have $|S^{\ast }R\left(\lambda \right)|\le S^{\ast }R\left(0\right)=1$. On the other hand $\frac{|\lambda +\mu +B|}{|\lambda +\mu |}>1$ and hence this equation does not have roots with nonnegative real parts. The second equation is a quadratic equation: ${\lambda }^2+\left(\mu -\tilde{\sigma }+\frac{r_2{N}_2}{K_2}\right)\lambda +(\mu -\tilde{\sigma })\frac{r_2{N}_2}{K_2}-E=0.$ In order for this equation to have only roots with negative real parts, the coefficients must be positive. We consider the constant term. The expression $(\mu -\tilde{\sigma })\left(r_2{N}_2/K_2\right)-E$ has the same sign as $\mu -\tilde{\sigma }-\frac{E}{\frac{r_2{N}_2}{K_2}}=\mu -J\left(0\right)=\mu (1-\rho (0\left)J\right(0\left)\right)>0.$ This also implies that $\mu -\tilde{\sigma }>0$. This completes the proof.

6. Discussion

In this paper we define a novel immuno-eco-epidemiological model of species competition. The model is based on a Lotka–Volterra competition model, where species one is infected by a pathogen, that does not affect species two. Competition of species in nature where only one of the species is affected by a pathogen is rare, but examples can be found. For instance, cutaneous fibromas caused by papillomaviruses affect white-tailed deer. The disease is restricted to deers and does not affect sheep, which are a natural competitor for food [9,13]. A broader application of this model is as a threshold case modelling well the circumstances when one of the species is seriously affected by the disease, while the other is affected little. This, for instance occurs often in competition between wild life and domestic animals, since domestic animals, even if subjected to a pathogen, are often treated and there is less implication of the pathogen to their well-being. Recognizing the importance pathogen affecting one species in competition with another, Anderson and May first considered this scenario [4]. The novelty here is that the second species interferes with the first not only in the competition for resources but also though lowering its ability to effectively combat the pathogen on within-host level.

Infected individuals are structured by age-since-infection, and the within-host dynamics between the pathogen and the immune response is taken into account. A novel feature of the model is the dependence of the within-host dynamics of species one on the numbers of species two. Thus, we incorporate competition on two levels: ecological and within host, where the competition from species two obstructs species one from mounting an effective immune response against the pathogen.

The within-host model, at an equilibrium of species two, has three equilibria: extinction equilibrium, which is always unstable, pathogen-only equilibrium, which is stable if the immune response reproduction number ${\mathcal{R}}_0<1$ and pathogen and immune response equilibrium, which is stable if ${\mathcal{R}}_0>1$. The immune response reproduction number depends on the equilibrial level of species two; hence, the more individuals in species two, the weaker immune response will species one mount.

The age-since-infection structured eco-epidemiological model has six equilibria, whose existence is controlled by six reproduction numbers. Eco-epidemiological extinction equilibrium $E_0$ always exists. There are also three disease-free equilibria: one corresponding to healthy species one $E_1$, one corresponding to species two, $E_2$ and one corresponding to the coexistence of healthy species one and species two $E_0^{\ast }$ which exist under appropriate conditions of the reproduction numbers. Finally, there is a unique species one-disease equilibrium $E_{N1}$ and (potentially multiple) disease-coexistence equilibria $E^{\ast }$, which also exist under appropriate conditions on the reproduction numbers.

To investigate the stability of the eco-epidemiological equilibria, we again assume that the within-host system is computed at the equilibrial level of species two. We find that the extinction equilibrium is locally asymptotically stable if the max of the species one and species two basic reproduction numbers is below one, $R_0<1$, and unstable if $R_0>1$. Species one equilibrium $E_1$ is stable if the basic reproduction number of species one is above one, but $R_0^{N_1}<1$ and $R_{N1N2}<1$. Species two equilibrium $E_2$ is stable if the basic reproduction number of species two is above one but $R_{N2}<1$. The disease-free coexistence equilibrium $E_0^{\ast }$ is stable whenever $R_0^{N_1}<1$. Species one-disease equilibrium $E_{N1}$ is locally asymptotically stable if $R_{N1}<1$ and $\sigma =r_1(1-2N_1/K_1)>0$ with $\mu -\sigma >0$. If $\sigma >0$ but $\mu -\sigma <0$ this equilibrium is unstable (saddle). In the case of a non-fatal disease, that is $\alpha \left(\tau \right)=0$, the disease-present coexistence equilibrium is locally asymptotically stable whenever it exists.

In summary, this is a complex model with multiple equilibria and conditions for existence and stability of these equilibria. The ecological competition of the species is affected by the within-host dynamics of the pathogen through the effect of the within-host dynamics exercised on the threshold quantities governing the outcome of the competition of the species and the impact of the disease.

Acknowledgments

The authors also acknowledge the extremely helpful remarks of one referee which lead to improvement of the paper. This article has benefitted from discussion with Michael Barfield and Robert D. Holt.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

Maia Martcheva acknowledges partial support from NSF grant [DMS-1515661].