Competitive exclusion in a multi-strain immuno-epidemiological influenza model with environmental transmission

ABSTRACT In this paper, a multi-strain model that links immunological and epidemiological dynamics across scales is formulated. On the within-host scale, the n strains eliminate each other with the strain having the largest immunological reproduction number persisting. However, on the population scale, we extend the competitive exclusion principle to a multi-strain model of SI-type for the dynamics of highly pathogenic flu in poultry that incorporates both the infection age of infectious individuals and biological age of pathogen in the environment. The two models are linked through the age-since-infection structure of the epidemiological variables. In addition the between-host transmission rate, the shedding rate of individuals infected by strain j and the disease-induced death rate depend on the within-host viral load. The immunological reproduction numbers and the epidemiological reproduction numbers are computed. By constructing a suitable Lyapunov function, the global stability of the infection-free equilibrium in the system is obtained if all reproduction numbers are smaller or equal to one. If , the reproduction number of strain j is larger than one, then a single-strain equilibrium, corresponding to strain j exists. This single-strain equilibrium is globally stable whenever and is the unique maximal reproduction number and all of the reproduction numbers are distinct. That is, the strain with the maximal basic reproduction number competitively excludes all other strains.


Introduction
Human infections with new avian influenza A(H7N9) virus were first reported in China in March 2013. Most of the H7N9 infections are believed to result from exposure to infected poultry or contaminated environments. These infections have raised our awareness of potential deadly pandemic that highly pathogenic flu viruses may create. Avian influenza (AI) is an infectious viral disease of birds (especially wild fowl such as ducks and geese) often causing no apparent signs of illness. AI viruses can sometimes spread to domestic poultry and cause large-scale outbreaks of serious disease. Some of these AI viruses, such as A(H5N1) and A(H7N9) viruses, have also been reported to cross the species barrier and cause disease or sub-clinical infections in humans [13] and other mammals. To understand the transmission of highly pathogenic flu among poultry and its potential ability to transform into a deadly human-to-human transmissible strain, we investigate a model of highly pathogenic flu in poultry.
We assume that the death rate of infected poultry depends on the infected within-host individual's viral load. This assumption allows us to study the impact of the viral load on the reproduction number and the disease prevalence. These are important questions in terms of providing insight for disease control and understanding characteristics of the disease. For example, this assumption was applicable to immunological and epidemiological dynamics of HIV. In addition, the persistence and the pandemic threat of avian influenza as well as the very publicized cholera outbreak in Haiti have increased the awareness of diseases that transmit both directly and indirectly through environment. Therefore, in this paper, we formulate a multi-scale immuno-epidemiological of flu viruses with direct and environmental transmission.
The importance of linking mathematical immunology and mathematical epidemiology cannot be underestimated particularly in the light of pathogen evolution and competition [2,16]. Since the seminal paper of Gilchrist and Sasaki [12] that nests a simple within-host model within an SIR epidemic model, a number of authors have considered the connection between these two scales, accounting for the epidemiological scale mostly through the epidemiological reproduction number [1,17,23]. Another approach that commonly has been applied is to link the within-host dynamics with the between-host dynamics when both are described by ordinary differential equation models. Such approach has been used in a number of articles but it typically does not separate the different time scales on which immunological and epidemiological processes progress [8].
In this paper, we use the dynamical approach in [12,21] to formulate a multi-scale immuno-epidemiological model of flu with different time scales. Here we assume that the high pathogenic avian influenza viruses come from the wild birds. HPAI viruses often kill birds within two days of onset of symptoms with mortality rate 90-100%. We then suppose the whole course of the infections is very short, so we do not include superinfection in the immune model. We further assume that all the HPAI viruses are distinct, therefore no co-infection occurs in the within-host model. In the present article, we investigate a scenario where the n strains are subjected to a complete competitive exclusion on the within-host scale. The strain with the maximal immunological reproduction number dominates. Because in our framework the strain with the higher within-host reproduction number takes over almost instantaneously, individuals in the population are infected with only one of the strains.
We aim to show the dependence of the epidemic reproduction number and disease prevalence on the infected within-host individual's viral load just when the poultry begins to be infected with the HPAI viruses and the infection begins to prevail. Thus, we do not consider super-infection in the epidemiological model. The epidemiological model we introduce here, describes the distribution of highly pathogenic flu among poultry. Since H5N1 is deadly to poultry, and chickens do not recover from it, we use an SI epidemic model with environmental transmission. The within-host model includes recruitment and clearance of target cells and allows for an infectious equilibrium.
There have been various approaches developed for continuous age-structured models that are described by first-order PDEs. The general idea is to study the nonlinear semigroup generated by the family of solutions. One approach is to use the theory of integrated semigroups [11,28]. Another method is to integrate solutions along the characteristics to obtain an equivalent integro-differential equation [6,30], and references therein. The goal of this paper is to extend the complete competitive exclusion result established by Bremermann and Thieme [5] to a multi-strain immuno-epidemiological flu model with environmental transmission. Mathematically speaking, the proof of a competitive exclusion principle is based on a proof of global stability of single-strain equilibrium. The stability analysis of nonlinear dynamical systems has always been a topic of both theoretical and practical importance since global stability is one of the most important issues related to their dynamic behaviours.
McCluskey and others, who have incorporated an integration term into a Lyapunov functional form often utilized for Lotka-Volterra type ODE models [20], have paved the way to establish global stability of the endemic equilibrium. The application of the Lyapunov function in age-structured models [6,11], and references therein, requires more delicate analysis than the case of ODE models. This often entails proving asymptotic smoothness of the semigroup generated by the family of solutions and proving existence of an interior global attractor [14], and then defining a Lyapunov function on this attractor. To show the Lyapunov function valid, uniform persistence of the single strain with the maximal reproduction number must be established. In [6,11], and references therein, authors have applied the persistence theory for infinite-dimensional systems by Hale and Waltmamm [15].
We use the immuno-epidemiological model that we formulate to address two questions: First, do periodic solutions exist in the immuno-epidemiological model of high pathogenic avian influenza (HPAI)? Our model takes the dependence of the epidemiological rates on the viral load explicitly. For the simplest dependence of the transmission and diseaseinduced death rate on the dynamic within-host viral load, the single-strain equilibrium is globally asymptotically stable and oscillations do not occur. Thus we extend the competitive exclusion principle to a multi-strain model of SI-type for the dynamics of highly pathogenic flu in poultry that incorporates both the infection age of infectious individuals and biological age of pathogen in the environment.
The second question that we address is: How the infected within-host individual's viral load affect the epidemiological reproduction number and the disease prevalence? We find that the epidemiological reproduction number increases with the viral load, and decreases with the clearance rate δ j (θ ) of the virus strain j of age θ from the environment. In addition, we give another threshold R * j linked to the epidemiological reproduction number R j . We show that the prevalence seems to increase with the increase in the viral load if R * j > R j > 1. But if R j > R * j > 1, we draw a conclusion that increasing the viral load will decrease the disease prevalence. The prevalence seems to depend on the viral load in an exactly opposite manner. It is an interesting trade-off scenario which can be observed in HIV.
This paper is organized as follows. In Section 2, we formulate a standard immunological and a standard SI epidemiological models with direct and environmental transmission. The two models are linked through age-since-infection and through the epidemiological parameters which depend on the within-host viral load. In Section 3, we discuss the trivial, semi-trivial equilibria and the global stability of the disease-free equilibria. In Section 4, we show the existence of C 0 semigroup generated by solutions of the model. We prove some important propositions of the semigroup. In Section 5, we construct a Lyapunov function to show that under the assumption of S(t) > ε, i 1 (τ , t) > ε, W 1 (θ , t) > ε, the complete orbit z(t) is the equilibrium E 1 . In Section 6, we prove that competitive exclusion occurs. How the within-host individuals' HPAI virus load impacts the disease prevalence and the reproduction number, is revealed in Section 7. Finally, a brief discussion is given in Section 8.

The model formulation
In this section we formulate an n-strain immuno-epidemiological model of flu with direct and environmental transmission. Individuals infected with strain j experience within-host dynamics. Within a host, the time variable is denoted by τ and signifies time-sinceinfection. Upon infection, flu virus attacks the lymphocytes, which we call the target cells and we denote them by x when healthy and by y n j when infected by the jth strain of the virus. The number of free virions of strain j is denoted by V j .
As H5N1 for poultry ultimately leads to death rather than recovery, in our study, unlike most of the models of flu [3], we consider a model incorporating target cell recruitment r and clearance m. Uninfected target cells are infected by interactions with free virions of strain j according to a mass action law with rate constant β n j . Productively infected cells, y n j , produce new virions at rate ν n j of strain j and die at rate d n j for strain j. Viral clearance rate is given by δ n j for strain j. The parameter s n j is the shedding rate of strain j. The within-host model of strain j is given by the following system of ordinary differential equations: Variations of Equation (1) have been studied extensively and have been commonly used in the investigation of HIV [25,26] and Hepatitis C dynamics [24]. Mathematically, the HIV/HCV version of Equation (1) has been completely analysed [10]. The within-host reproduction number of strain j in Equation (1) is computed as follows: The reproduction number of strain j gives the number of secondary infections that virus j-infected cell will produce in an entirely healthy cell during its lifespan as infected.
Equation (1) always has an infection-free equilibrium E 0 = (r/m, 0, . . . , 0), and if R j > 1, then the infected equilibrium E j = (x * n j , 0, . . . , 0, y * n j , 0, . . . , 0, 0, . . . , 0, V * j , 0, . . . , 0) exists with The long-term outcome of the competition of the n strains in Equation (1) is competitive exclusion. The strain with the larger reproduction number wins and eliminates the strain with the smaller reproduction number [9]. We now introduce an n-strain epidemiological model with direct and environmental transmission. We assume that the pathogen causing the disease is represented by multiple strains. To introduce the model, we denote by S(t) the number of susceptible poultry individuals, where t is the chronological time. We structure the infected individuals by agesince-infection τ . Let i j (τ , t) be the density of individuals infected by strain j. Individuals in the class i j (τ , t) experience the same within-host dynamics given by Equation (1) for strain j. We assume that for individuals in class j, the immunological reproduction number R j is maximal so strain j is the main strain that infects them. Let W j (θ , t) denote the number of virions of age θ ≥ 0 of strain j at time t in the environment.
With the above notations, we study the following infection-age-structured model with environmental transmission.
Let be the birth/recruitment rate, μ be the natural death rate of susceptible hosts at zero viral load, μ j V j (τ ) gives the additional host mortality due to the virus, and here we assume the simplest form. Gilchrist and Sasaki [12] use a different form but our within-host model does not involve the immune response. We also assume that all susceptible individuals have approximately the same equilibrium level of healthy lymphocytes. The transmission coefficient of β j (τ ) of strain j and the shedding rate of η j (τ ) of an infected individual of infection age τ infected by strain j are also dependent on the within-host viral load. In influenza, infectivity and shedding of the virus are associated with higher within-host viral load. For simplicity, we may assume that β j (τ ) and η j (τ ) are proportional to the viral load at a given age-since-infection τ , although other functional forms are also possible. Hence, where V j (τ ) is the viral load of the strain j. l j and k j represent the probability of successful transmission. s j is the shedding rate of strain j. ξ j (θ ) is the transmission rate of strain j from the environmental contamination. δ j (θ ) is the clearance rate of the virus strain j of age θ from the environment.
Here the susceptible individuals are recruited at a rate . Susceptible individuals can become infected with strain j either through a direct contact with an infected individual with strain j or through coming into contact with viral particles of strain j that are in the environment. Infection through direct contact with infected individuals can happen through contact with individuals of any age-since-infection at a specific age-specific transmission rate. As a consequence, the force of infection of susceptible individuals through direct contact is given by the integral over all ages-since-infection. So is the force of infection of susceptible individuals through the contaminated environment. Upon infection through direct or indirect transmission, the newly infected individuals move to the infected with strain j class with age-since-infection equal to zero, that is, they move to the boundary condition. Infected individuals with strain j shed the virus into the environment at a rate η j (τ ). All viral particles shed by individuals infected with strain j of all age classes are given by the integral. This is the number of virions with strain j class with age-since-infection equal to zero at time t in the environment.
Model (3) is equipped with the following initial conditions: Initial conditions are determined in advance and do not depend on the within-host model. All parameters are nonnegative, > 0, μ > 0, and μ j > 0. We make the following assumptions on the parameter functions.
Assumption 2.1. The parameter functions satisfy the following: (1) The functions β j (τ ) and ξ j (θ ) are bounded and uniformly continuous for every j. When β j (τ ) and ξ j (θ ) are of compact support, the support has non-zero Lebesgue measure.
Define the space of functions Note that X is a closed subset of a Banach Space, and hence is a complete metric space. The norm on X is taken to be: It can be verified that solutions of Equation (3) with nonnegative initial conditions belong to the positive cone for t ≥ 0. Furthermore, adding the first and all equations for i j , we have The free virus in the environment can be bounded as follows: Therefore, the following set is positively invariant for system: Finally, since the exit rate from the infectious compartment is given by μ + μ j V j (τ ), the probability of still being infectious after τ time units is given by where π 2 j (θ ) denotes the probability of the virus with strain j of age θ still being in the environment.
The reproduction number of strain j in system (3) is given by the following expression: The reproduction number of strain j gives the number of secondary infections that one strain j-infected individual will produce in an entirely susceptible population. R j gives the strength of strain j to invade when rare and alone. The reproduction number R j consists of two terms. The first term accounts for the number of secondary infections that one strain j-infected individual will produce through direct transmission. The second term accounts for the number of secondary infections that one strain-j-infected individual will produce through environmental transmission. One may define a reproduction number of the whole system: In the next section we compute explicit expressions for the equilibria and establish the global stability of the disease-free equilibrium.

Equilibria and the global stability of the DFE
System (3) always has a unique disease-free equilibrium E 0 , which is given by where 0 = (0, . . . , 0) is n-dimensional vector of zeroes. In addition, for each j there is a corresponding single-strain equilibrium E j given by where the non-zero components i * j and V * j are in positions 2j and 2j+1, respectively. The endemic equilibrium E j exists if and only if R j > 1. The non-zero components of the equilibrium E j are given by where Next, we turn to the linearized equations of the disease-free equilibrium. To introduce the linearization at the disease-free equilibrium E 0 , let To study system (7), we notice that the system for y j and z j is decoupled from the equation for x. Hence, the equations for each j are independent from the equations for the other strains. We look for solutions of the form x(t) =x e λt , y j (τ , t) =ȳ j (τ ) e λt and z j (θ , t) = z j (θ ) e λt . For each j we obtain the following eigenvalue problem: Solving each differential equation, we obtain Similarly, substituting forȳ j (τ ) andz j (θ ) in the equation forȳ j (0) and then dividingȳ j (0) from both side of the resulting equation, we get the following characteristic equation for strain j: Now we are ready to establish the following result.

Proposition 3.1: If
then the disease-free equilibrium is locally asymptotically stable. If R 0 > 1, the disease-free equilibrium is unstable.
Consider λ with λ ≥ 0. For such λ and each j, we have Hence, the equations G j (λ) = 1 do not have a solution with nonnegative real part and the disease-free equilibrium is locally asymptotically stable. Now assume For that fixed j 0 and λ real, we have Hence, the equation G j 0 (λ) = 1 has a real positive root. Therefore, the disease-free equilibrium is unstable.
We have established that the disease-free equilibrium is locally stable, i.e. given the conditions on the parameters, if the initial conditions are close enough to the equilibrium, the solution will converge to that equilibrium. In the following paragraphs our objective is to extend these results to global stability results. That is, given the conditions on the parameters, convergence to the equilibrium occurs independent of the initial conditions.
We will use a Lyapunov function to establish the global stability of the disease-free equilibrium. We need to integrate the differential equation along the characteristic lines. Denote the initial condition by B j (t): Integrating along the characteristic lines, we obtain

Proposition 3.2: Assume
Then the disease-free equilibrium E 0 is globally asymptotically stable.
Proof: We will use a Lyapunov function. We adopt the logistic function used in [20,22]. Define We note that f (x) ≥ 0 for all x > 0. f (x) achieves its global minimum at one, with f (1) = 0. Let We denote S * 0 = /μ. We define the following Lyapunov function: dθ .
Before differentiating, we rewrite the middle term using Equation (10).
Changing variables, we have Differentiating above, we have dθ .
Noting that we have Merging some middle terms, we have Note that Canceling all terms that cancel, we simplify the above expression as The last inequality follows from the fact that R 0 ≤ 1. Notice that V equals zero if and only if the first term equals zero, i.e. S = S * 0 , and if each of the terms in the sum is zero. We define a set LaSalle's invariance principle [18] implies that the bounded solutions of Equation (3) converge to the largest compact invariant set of 1 . We will show that this largest compact invariant set is the singleton given by the disease-free equilibrium. First, we saw that S = S * 0 . From the first equation in (3), we have that Noting that all parameters are non-negative, we have Thus, from the solutions for the equations along the characteristic lines (10), we have that i j (τ , t) = 0 for all t > τ. Hence, for t > τ, we have Similarly, we also have We conclude that the disease-free equilibrium is globally stable. This completes the proof.
From the above discussion, we saw that if all reproduction numbers are less or equal to one, all strains are eliminated and the disease dies out. In the next section, we will show the existence and properties of semigroup.

Existence and properties of semigroup
Let us define the transformation Let x 0 ∈ . For any neighbourhood B 0 ⊂ with x 0 ∈ B 0 , in accordance with the contraction mapping theorem, for every sufficiently small ε > 0, there exists a unique contin- Here existence and uniqueness can be proved by formulating the solutions to system (3) as a fixed point of an integral operator ψ on an appropriate closed subset of For t ≥ 0, define the continuous mapping From the discussion in Section 2, solutions to system (3) can be shown to remain bounded in forward time, then existence of the semigroup {T(t)} | {t≥0} can be established. The family of functions {T(t)} | {t≥0} satisfies the properties of a C 0 semigroup on [6,14], i.e. to say, Moreover, the semigroup T(t) is point dissipative. Hence, we suppose that the solutions are forward complete, i.e. the solutions exist on the time interval [0, ∞).
In order to prove that system (3) has a global compact attractor T 0 , we need to establish asymptotic smoothness of the semigroup. As a first step, we define the semiflow ψ of the The semiflow is a mapping ψ : [0, ∞) × → . A set T 0 in is called a global compact attractor for ψ, if T 0 is a maximal compact invariant set and if for all open sets U containing T 0 and all bounded sets B of there exists some t > 0 such that ψ(t, B) ⊆ U, for all t > t. We can establish the following proposition.

Proposition 4.1: The semigroup T(t) is asymptotically smooth.
Proof: To establish this result, we will apply Lemma 3.1.3 and Theorem 3.4.6 in [14]. To show the assumptions of Lemma 3.1.3 and Theorem 3.4.6 in [14], we split the solution semiflow into two components. For an initial condition The splitting is done in such a way thatψ(t, x 0 ) → 0 as t → ∞ for every x 0 ∈ , and for a fixed t and any bounded set B in , the set {ψ(t, x 0 ) : x 0 ∈ B} is precompact. The two components of the semiflow are defined as follows: where are the solutions of the following equations (the remaining equations are as in system (3)) and System (17) is decoupled from the remaining equations. Using the formula (10) to integrate along the characteristic lines, we obtain , t < τ.
Integratingî j with respect to τ , we obtain This shows the first claim, i.e. it shows thatψ(t, x 0 ) → 0 as t → ∞ uniformly for every x 0 ∈ B ⊆ , where B is a ball of a given radius.
To show the second claim, we need to show compactness. We fix t and let x 0 ∈ . Note that is bounded. We have to show that for that fixed t the family of functions defined bỹ obtained by taking different initial conditions in is a compact family of functions. The family and, therefore, it is bounded. Thus, we have established the boundedness of the set. To show compactness, we first see that the third condition in the Frechét-Kolmogorov theorem [31] is trivially satisfied sinceĩ j (τ , t) = 0,W j (θ , t) = 0 for min{θ , τ } > t. To see the second condition, we have to bound by a constant the L 1 -norm of ∂i j /∂τ , ∂W j /∂θ . To derive that bound, first notice thatĩ First, we notice that for x 0 ∈ ,B 1 j (t) is bounded. We can observe this by recalling that S is bounded. Hence, theB 1 j (t) satisfies the following inequality: Here k 1 , k 2 and k 3 are constants that depend on the bounds of the parameters as well as the bounds of the solution. Gronwall's inequality implies that In the following, we derive that for x 0 ∈ ,B 2 j (t) is also bounded.
Next, we differentiate Equation (20) with respect to (τ , θ): We have to see that |(B 1 j (t − τ )) | and |(B 2 j (t − θ)) | are bounded. Differentiating Equation (21), we obtain Similar to Equation (22), we can rewrite and simplify the above equality as the following inequality: Gronwall's inequality then implies that, It is evident that Putting all these bounds together, we have To complete the proof, we notice that Thus, the integral can be made arbitrary small uniformly in the family of functions. That establishes the second requirement of the Frechet-Kolmogorov theorem. We conclude that the family is compact. Therefore, we prove that T(t) is asymptotically smooth.
Now we have shown that the semigroup T(t) generated by the system (3) is asymptotically smooth and point dissipative on the state space . We also notice that the forward orbit of bounded sets is bounded in . Then, in accordance with the sufficient condition by Hale [14], we have all components to establish the following proposition. Next, we recall several definitions concerning semigroup dynamics in infinite dimension.
A complete orbit through x is a function z : R → with z(0) = x and, for any s ∈ R, T(t)z(s) = z(t + s) for every t ≥ 0. The omega limit set of x, ω(x), is defined as The alpha limit set corresponding to the complete orbit z(t) through x is denoted by α z (x), and defined to be the following: If there exists a backward orbit z(t) through x, the unstable manifold is defined by Now, we will prove two propositions concerning limits sets in forward and backward time, respectively. First, we will prove a simple result about the stable manifold of the singleton E i , W s ({E i }), which will be applied later in the proof of uniform persistence for our system.
Suppose by way of contradiction, ∃ > 0, t n → ∞ such that T(t n )x − E j ≥ . As shown in the proof of Proposition 4.1, the semiflow T(t) can be written as T(t) =T(t) +T(t). Since {T(t n )} is pre-compact, there exists a convergent subsequence: which is a contradiction to the definition of the stable manifold.
Second, we consider the alpha limit set corresponding to a complete orbit corresponding to solutions of system (3). The following result is utilized in the application of a Lyapunov functional to our system.
is a complete orbit through x. Integrating along the characteristics, we get the following more general formulation: Modifying the arguments in the proof of Proposition 4.1, we can prove that Clearly the convergence is uniform ∀t ∈ R, so {z(t) : t ∈ R} is pre-compact. Then α z (x) is non-empty, compact. The remaining conclusions follow from Theorem 2.48 in [27].
In order to prove the competitive exclusion principle, we need to use a Lyapunov function to establish that the complete orbit through x is the equilibrium E 1 under some stronger condition.

Lyapunov function
We will consider the case where all the strains have different reproduction numbers which are greater than one. If the reproduction number R j , j = 1, 2, . . . , n are smaller than one, establishing the extinction of strain j can be approached with other more direct methods. Therefore, all of the following results hold for the case where some strains have reproduction number less than unity. But in order to make the notation simpler, we assume that min i R i > 1.
Without loss of generality, we assume that We consider complete orbits for our system to analyse the global dynamics by a Lyapunov function. Suppose that a complete orbit z(t) is defined as in Proposition 4.4 and it suffice system (24). In the proof of the following proposition, we find a Lyapunov function for a complete orbit z(t), which is well-defined and bounded when z(t) satisfies certain criteria, namely z(t) is bounded from above and bounded away from an appropriate boundary set. Under these criteria, a LaSalle invariance-type argument can be invoked to show that the complete orbit z(t) must be in fact the equilibrium E 1 .

Proof:
We expect to show this result using a Lyapunov function, similar to the one used in [4,20,22]. With f (x) = x − 1 − ln x, we define the following Lyapunov function: One difficulty with the Lyapunov function U above is that the component U S is not defined if S = 0, the component U I is not defined if i 1 (τ , t) = 0 on a set of non-zero measure, and component U W is not defined if W 1 (θ , t) = 0. To show that the Lyapunov function above is valid, according to the assumption, we can get From the above equality, we have Also, we have This makes the Lyapunov function defined in Equation (26) well-defined. Because of the complexity of the expressions, we differentiate each component of the Lyapunov function separately (see Equation (26)).
Next, we need to take the time derivative of U I . Before we do that, we need to transform U I in such a way that it is more convenient for differentiation. We will use the representation of i 1 , given in Equation (10) [20]: Differentiating the last expression above, we have dτ . (29) The above equality follows from Equation (26) and the fact Hence, U I can be simplified as follows: Now we turn to the derivative of the environmental component. Similarly, we first transform U W to make it more convenient for differentiation. We will use the representation of W 1 , given in Equation (10) [20]: Similar to Equation (29), we differentiate the last expression above and we have The above equality follows from Equation (26) and the fact Hence, dθ .
(33) Finally, we consider the expression that corresponds to strains two to n. First, we rewrite U θτ (t) to prepare it for differentiation.
Similar to Equation (29), differentiating the last expression with respect to t, we have Adding all four components of the Lyapunov function, we have Notice that Using Equations (37) and (38), Equation (36) can be simplified as follows: Next, we notice that Indeed, Using Equation (40) to simplify (39), we obtain Hence, U (t) ≤ 0.
is a non-increasing map, which is bounded above, we conclude that U • (z(t)) → c < ∞ as t → −∞. Therefore, U • (x) = c for allx ∈ α z (x). Combining this with the fact that α z (x) is invariant, we get that U • (g(t)) = c for all t ∈ R, where g(t) is a complete orbit througĥ x with g(0) =x. Hence, (d/dt)U(g(t)) = 0 for all t ∈ R. This implies that {g(t), t ∈ R} is an invariant set with the property that dU/dt = 0. Therefore, g(t) = E 1 for all t, in particular when t = 0. So,x = E 1 . This show that α z (x) = {E 1 }. Thus, U • (z(t)) ≤ U • (E 1 ) for all t ∈ R. Since E 1 is the unique minimizer of U, z(t) = E 1 , ∀t ∈ R, i.e. x = E 1 .

Competitive exclusion
Proposition 5.1 states that the only complete orbit z(t) in an appropriate subset is the equi- (3). If we find a global attractor on this appropriate subset, then due to its invariance, the global attractor will reduce to the equilibrium E 1 . To follow this strategy, we apply the persistence theory by Hale and Waltmann [15] for infinite-dimensional systems and a result from Magal and Zhao [19] on existence of an interior global attractor. The methods and techniques have been recently employed by other authors. [6,11].
Persistence theory provides a mathematical formalism for determining whether a species will ultimately go extinct or persist in a dynamical model. Consider as the closure of an open set X 1 , that is, = X 1 ∪ ∂X 1 , where ∂X 1 (assumed to be non-empty) is the boundary of X 1 . Assume that the C 0 semigroup T(t) on satisfies Let T| ∂ := T(t)| ∂X 1 and A ∂ be the global attractor for T| ∂ . The particular invariant sets of interest areÃ The following result is discussed in [15, Theorem 4.2]: Theorem 6.1: Suppose that T(t) satisfies Equation (44) and the following conditions: (iv)Ã ∂ is isolated and has an acyclic coveringM, wherẽ Then T(t) is a uniform repeller with respect to ∂X 1 , i.e. there is an ε > 0 such that The following theorem relates uniform persistence to existence of a global attractor in X 1 . [19]): Assume that the semigroup T(t) satisfying Condition (44), is asymptotically smooth and uniformly persistent, and has a global attractor T 0 . Then the restriction of T(t) to X 1 , T(t)| X 1 , has a global attractor T.

Theorem 6.2 (Magal and Zhao
In order to proceed, we need to be precise about considering various forward invariant subset of X. Then, we can define our uniformly persistent set and complementary boundary, and utilize mathematical induction to characterize the dynamics on the boundary set.
For j = 1, 2, . . . , n, define the following sets: Note thatτ j ,θ j are allowed to be infinity. Proposition 6.1: For j = 1, 2, . . . , n, X j and ∂X j are forward invariant under the semigroup T(t). Also, X 0 is forward invariant, and if x ∈ X 0 , then T(t)x → E 0 as t → ∞. In addition, . From the properties of the C 0 -semigroup T(t), we have T(t 2 + ε) = T(t 2 )T(ε), i.e. T(t 2 + ε)x is a solution to system (3) with initial condition T(ε)x for every t 2 ≥ 0. In accordance with the following inequality, we have .
As a result, we get In case (ii), similarly to the above proof, we can get By i j (0, t) ∈ M 0 j , we can obtain i j (τ , t) ∈ M 0 j using the same argument as in case (i). This shows forward invariance for j = 1. When j > 1, notice that ∂X j−1 is forward invariant, which implies forward invariance of X j . Further, we have T(t)X j ⊂ (X j ) + , ∀t > 0, j = 1, 2, . . . , n.
Since X j ∈ X 0 , j = 1, 2, . . . , n, we conclude that X 0 is forward invariant. Also, X 0 = ∂X n is forward invariant. In view of our system, it is clear that ∀x ∈ X 0 , we have T(t)x → E 0 as t → ∞, where E 0 is the infection-free equilibrium.
Our next goal is to prove the main result, i.e. solutions with initial conditions corresponding to positive productive infected individuals i j (τ , t) or positive concentration of virion in the environment W j (θ , t) will converge to the equilibrium E 1 (the single-strain equilibrium belonging to the strain with maximal basic reproduction number).

Proposition 6.2:
Assume that R 1 > R 2 > · · · > R n . Then, E 1 is globally asymptotically stable for system (3) with respect to initial conditions satisfying Proof: We prove the proposition by induction on the number of strains, n, in system (3). n = 1: We note that the proof for this case is contained inside following arguments. We will not repeat that proof but we will simply state that the case n = 1 was proved in [7]. Induction Step: Assume that Proposition 6.2 is true for all n < m. We will prove the proposition is true for n = m. Lemma 1: If x ∈ X l where l ≥ 2, then T(t)x → E l as t → ∞.
Proof: For l = 2, 3, . . . , m, define the projection operator Furthermore, define a semigroup of the projected system as follows: T l (t) is the semigroup on R + × m−l+1 1 (L 1 + × L 1 + ) generated by the solutions to system (3) with n = m−l+1 strains, which matches the m strains model except that the first l−1 strains are eliminated. Then, X l is projection invariant with respect to system (3), i.e. P l (T(t)x) = T l (t)P l (x), ∀x ∈ X l , t ≥ 0. Following from our induction hypothesis, we get that for

Lemma 2:
The semigroup T(t) is uniformly persistent with respect to X 1 and ∂X 1 . Moreover, there exists a compact set T ⊂ (X 1 ) + , which is a global attractor for {T(t)} t≥0 in X 1 , and ∃u > 0 such that Proof: We will apply Theorem 6.1. Let A ∂ be the strong global attractor of ∂X 1 . Noting that ∂X 1 = X 0 ∪ m l=2 X l andÃ ∂ := ∪ A ∂ w(x), from Lemma 1 and Proposition 6.1, we obtain thatÃ ∂ = {E 0 , E 2 , E 3 , . . . , E m }. We will show that each {E k } ⊂Ã ∂ , k = 0, 2, 3, . . . , m is an isolated invariant set. Let B := B r (E k ) be an open ball of sufficiently small radius r around E k . We claim that B is an isolating neighbourhood. Suppose by way of contradiction that there exists some κ ∈ {2, 3, . . . , m} or κ = 0 contenting with {E k } being not a maximal invariant set. Let M ⊂ B be an invariant set with M = {E k }. Then, there exists a complete orbit γ (x) ⊂ M for x ∈ M\{E k }. If x ∈ X κ ∪ X κ+1 ∪ · · · ∪ X m , by Proposition 5.1, x = E κ , which is a contradiction. If x ∈ Xk,k = 0, 2, 3, . . . , κ − 1, then x = E˙k by Lemma 1 and Proposition 6.1, again a contradiction since E˙k∈B.
To show thatÃ ∂ is acyclic, we need to institute that no subset ofÃ ∂ forms a cycle. From Lemma 1, we know that when x ∈ X l , T(t)x → E l as t → ∞. Hence, by Proposition 4.3 and the definition of stable manifold, ∀l = 2, 3, First, Let us consider the possibility of a cycle with length greater than or equal to 2. This cycle must include a chain with For the above b and l, if the cycle is a chain with {E l } → {E b }, then, for x ∈ W u ({E l }) ∩ W s ({E b }), x ∈ X 0 or x ∈ X l , l = 2, 3, . . . , m − 1. If x ∈ X 0 , it means i l (τ , 0) ∈ ∂M 0 l and W l (θ , 0) ∈ ∂N 0 l . Similarly to the argument in the above paragraph, we have x∈W u ({E l }). If x ∈ X l , i l (τ , 0) ∈ M 0 l or W l (θ , 0) ∈ N 0 l . Since X l is forward invariant and T(t)X l ⊂ (X l ) + , for any negative t on a backward orbit through x, i l (τ , t) ∈ M 0 l or W l (θ , t) ∈ N 0 l . In accordance with the continuity of i l (τ , t) and W l (θ , t), lim t→−∞ i l (τ , t) ∈ M 0 l or lim t→−∞ W l (θ , t) ∈ N 0 l . That is to say, x∈W u ({E l }), ∀l < l. This excludes the possibility of cycles of length greater than or equal to 2 for T(t)| ∂X 1 . Now we consider the possibility of a 1-cycle for T(t)| ∂X 1 . Then, {E k } → {E k } for some k = 0, 2, 3, . . . , m. First, we show that we can not have a 1-cycle for E 0 . It suffice to show that (X 0 \{E 0 }) ∩ W u ({E 0 }) = ∅. Let x ∈ X 0 \{E 0 }. Any backward orbit of x must stay in X 0 since X 0 is forward invariant. If i e (τ , 0) ∈ ∂M 0 e , W e (θ , 0) ∈ ∂N 0 e , ∀e = 1, 2, . . . , m, then we have a scalar ODE with a unique positive equilibrium and lim t→∞ T(t) = 0 or ∞.  0 W e (θ , t) dθ > 0 for some negative t on a backward orbit through x, which is a contradiction. Therefore, there can be no backward orbit through x if ∞ τ e i e (τ , 0) dτ > 0 or ∞ θ e W e (θ , 0) dθ > 0. Hence, E 0 can not be an α−limit point of any x ∈ X 0 \{E 0 }. Now consider the case x ∈ X l for some l ∈ {2, 3, . . . , m}. We prove by contradic- In a similar fashion, we can show S(t) ≥ ε for all t ∈ R. Since both X l and \X l are forward invariant, we can get z(t) ∈ X l , ∀t ∈ R. X l is a projection invariant with respect to P j and T j as defined earlier. In other words, in X l , we can consider an equivalent n = m−l+1 strain model with R l as the maximal reproduction number. In this case, Proposition 5.1 applies to P l (z(t)) and semigroup T l (t). We can conclude that x = E l , which is a contradiction. Hence,Ã ∂ is acyclic for T(t) | ∂X 1 .
To complete the proof of uniform persistence, we need to prove: Assume that on the contrary, there exists x ∈ X 1 with i 1 (τ , 0) ∈ M 0 1 or W 1 (θ , 0) ∈ N 0 1 such that x ∈ W s ({E k }) where k = 0, 2, 3, . . . , m. Then, T(t)x → E k as t → ∞. By proposition 6.1, X 1 is forward invariant under the semigroup T(t) and T(t)X 1 ⊂ (X 1 ) + , ∀t > 0. That is, i 1 (τ , t) ∈ M 0 1 or W 1 (θ , t) ∈ N 0 1 , ∀t > 0. Furthermore, we can get as t → ∞, It means that T(t)x E k . It is a contradiction. Thus, W s ({E k }) ∩ X 1 = ∅. By Theorem 6.1, we find that T(t) is uniformly persistent with respect to X 1 and ∂X 1 , i.e. ∂X 1 is uniform strong repeller. Then, by Theorem 6.2, we can conclude that there exists a compact set T ⊂ X 1 which is a global attractor for {T(t)}| t≥0 in X 1 . Since T(t)X 1 ⊂ (X 1 ) + , the global attractor, T, is actually contained in (X 1 ) + . Because of this, ∃u > 0 such that Because the interior global attractor T is invariant, we can find a complete orbit through any point contained in T. For any complete orbit z(t) : t ∈ R ⊂ T ⊂ (X 1 ) + , there exists ε > 0 such that S(t) > ε, i 1 (τ , t) > ε and W 1 (θ , t) > ε for all t ∈ R. Hence, by Proposition 5.1, T = E 1 . Thus, E 1 is the globally attractor for system (3) with respect to initial conditions satisfying A global attractor is also locally stable by definition, therefore, E 1 is indeed globally asymptotically stable.
Since the equilibrium value of the susceptible in the single-strain equilibrium E j is inversely proportional to the reproduction number of strain j, then the strain that can persist on the smallest number of the susceptible is the strain that persists in the population. That is, the competitive exclusion principle has been established in the context of system (3) and the strain with the maximal reproduction number eliminates all the rest.
We have established the global stability of strain one. This is in contrast with many other age-since-infection models where the age-since-infection has been found to destabilize the equilibrium and sustained oscillations were present [29]. It should be noted that the dramatic difference in the dynamical behaviour is due to the mass action incidence that we assume in model (3). In [29] and other articles where destabilization occurs, standard incidence has been used.

How does the infected individual's within-host viral load affect the disease prevalence?
In this section we show how the infected individual's within-host viral load affects the disease prevalence and the epidemiological reproduction number. First we need to transform the epidemiological reproduction number R j , to make it more convenient to visualize the link between it and the within-host viral load. We will use the representation of R j given in Equations (4) and (5).

Thus we have
It is evident that the epidemic reproduction number R j increases with the increase of the viral load and decreases with the decrease of the viral load. It is also easy to see that the epidemic reproduction number R j increases with the decrease of the clearance rate of δ j (θ ) of age θ from the environment and decreases with the increase of δ j (θ ). Next, we need to take the derivative of I j with respect to the viral load. Before doing that, we need to transform I j to make it more convenient for differentiation. I * j is given by we can obtain Taking the ρ derivative of I * j dI * Define R * j as follows: If R * j > R j > 1, we have dI * j /dρ < 0. So we obtain that increasing ρ decreases the disease prevalence I * j , while ρ decreases with the increase of the viral load V j (τ ). As a result, the disease prevalence I * j will increase with the increase of the viral load V j (τ ). If R j > R * j , we have dI * j /dρ > 0. So we obtain that increasing ρ increases the disease prevalence I * j , while ρ decreases with the increase of the viral load V j (τ ). We draw a conclusion that the disease prevalence I * j will decrease with the increase of the viral load V j (τ ). It is an interesting trade-off scenario.
The competitive exclusion result [5,22], i.e. the strain with the maximal reproduction number eliminates all the rest, has been successfully extended to a multi-strain immunoepidemiological flu model with environmental transmission. The results show that the basic reproduction number R 0 determines the transmission dynamics of HPAI diseases: the disease-free equilibrium is globally asymptotically stable if R 0 ≤ 1. If R 0 = R j > 1 is the maximal reproduction number, the dominant single-strain equilibrium E j is globally asymptotically stable, and the other strains die out. The results suggest that an effective strategy to contain HPAI disease is decreasing the reproduction number R 0 below one.
From the perspective of public health, public health efforts will work best if directed to monitoring the single strain j whose basic reproduction number is maximal. We can largely decrease the infected individuals by culling as culling is one of the main control measures used to decrease avian influenza in the poultry population. Then the reproduction number R j will decrease and the decrease in the AI disease prevalence I * j can be obtained if R * j > R j > 1. But under the condition of R * j < R j , although R j decreases with the decrease of the viral load V j (τ ), the disease prevalence I * j increases. This may explain why AI infection is still around despite implementation of integrated control measures.

Discussion
In this paper, we study a multi-scale immuno-epidemiological model of influenza viruses including direct and environmental transmission. This work extends prior results in [22] through the linking to within-host model and the integration of arbitrarily distributed clearance rates of the virus in the environment. Mathematically this extends the model in [22] to a model with two time-since-infection structural variables in contrast to most models which only have one structural variable.
We assume that the high pathogenic avian influenza viruses of the infected poultry are obtained from the wild birds. The within-host and the between-host models are linked through the age-since-infection structure of the epidemiological variables. In addition, rates of between-host transmission, shedding by individuals infected with strain j, and disease-induced death are dependent on the within-host viral load. The immunological and epidemiological reproduction number (R j , R j ), respectively, are then computed. We show that if all epidemiological reproduction numbers are lesser or equal to one, the disease-free equilibrium in the system is locally and globally asymptotically stable and unstable if at least one of them is greater than one. Furthermore, the dominance equilibrium E j of HPAI virus strain exists and is unique when the reproduction number R j is greater than one. Without loss of generality, we assume that R 1 > 1 is the maximal reproduction number and all of the reproduction numbers are distinct, and then we prove that the dominance equilibrium E 1 is globally asymptotically stable. In brief, under the assumption of the absence of super-infection and co-infection, on the within-host scale, the n strains eliminate each other with the strain with the largest immunological reproduction persisting. However, on the population scale, we extend the competitive exclusion principle to a multi-strain model of SI-type for the dynamics of highly pathogenic flu in poultry that incorporates both the infection age of infectious individuals and biological age of pathogen in the environment. We note that in the presence of super infection the complete competitive exclusion result that we obtain here may not hold.
In this study, we also address the question: How does the viral load of the infected within-host individuals affect the epidemic reproduction number and the disease prevalence? These are important questions in terms of providing insight on disease control and understanding characteristic of the disease. Calculations suggest that increasing the viral load V j (τ ) increases the transmission coefficient β j (τ ) and the reproduction number R j . If R * j > R j > 1, the disease prevalence I * j will increase with the increase of the viral load V j (τ ). But if R j > R * j , the disease prevalence I * j will decrease with the increase of the viral load V j (τ ). The increase of prevalence with decreasing viral load can be observed in the use of HIV medications in practice. In addition, the reproduction number R j decreases with the clearance rate δ j (θ ) of the virus strain j of age θ from the environment.
Under the assumptions that have been considered by us, our results suggest that the AI infection can be controlled well at the beginning of the disease if the infected poultry with HPAI virus who has already begun to show symptoms are all slaughtered in a timely fashion and also their carcasses are handled appropriately. It is important to keep the stall clean all the time and that helps in preventing the disease from spreading to some degree.
The primary risk factor for human infection appears to be a direct or an indirect exposure to the infected live or dead poultry or contaminated environments. The control of circulation of the H5N1 and H7N9 viruses in poultry is essential to reduce the overall risk to human infection.