The Interaction of Migratory Birds and Domestic Poultry and Its Role in Sustaining Avian Influenza

. We investigate the role of migratory birds in the spread of H5N1 avian inﬂuenza, focusing on the interaction of a migratory bird species with nonmigratory poultry. The model is of patch type and is derived with the aid of reaction-advection equations for the migratory birds in the air along the ﬂyways. Poultry may reside at some or all of the four patches of the model, which consist of the breeding patch for the migratory birds, their winter feeding patch, and two stopover patches where birds rest and refuel on their migration. Outward and return migratory routes can be diﬀerent. The equations for the migratory birds contain time delays which represent the ﬂight times for migratory birds along particular sectors. Delays also appear in the model coeﬃcients via quantities which represent ﬂight survival probabilities for the various sectors. We establish results on positivity, boundedness, global asymptotic stability of the disease-free equilibrium, and the persistence of infection. We also discuss extensions of the model to include the seasonality of the migration phenomenon. Numerical simulations support the analytical ﬁndings; here we used data on H5N1 infected ducks in the Poyang Lake region of China.


Introduction.
Of the various H5N1 introductions into different regions of the world, at least some were most likely due to migratory birds (Kilpatrick et al. [21]). Apart from migratory birds, other main causes of H5N1 introduction include the trade in poultry and poultry products and the trade in wild birds. In previous works (Gourley, Liu, and Wu [14], Bourouiba et al. [6]), we developed mathematical models of the bird migration phenomenon with the ultimate aim of understanding the role of migratory birds in spreading the highly pathogenic avian influenza (H5N1), leading to a massive outbreak in wild birds at Lake Qinghai in central China (2005). Bourouiba et al. [6] particularly focused on bar-headed (Anser indicus) geese migration in the Central Asian Flyway from Mongolia to India. Migratory birds encounter a variety of climatic and other conditions at their breeding and wintering locations and along their flyways and resting places (stopover sites), which led us initially to consider the possibility of a partial differential equation model. It also became apparent to us in the early stages of our work that even in the absence of disease dynamics a realistic model of bird migration would necessarily be a complicated one. Eventually, we settled on models of patch type, where the patches are the breeding and wintering L. BOUROUIBA, S. A. GOURLEY, R. LIU, AND J. WU locations and the various stopover sites along the migration flyway where birds stop to rest and eat (refuel). In between patches the migration flyway is a one-dimensional continuum along which bird density is modelled using reaction-advection equations. In our previous works we have allowed for the fact that migration is essentially a periodic phenomenon. In the model derivation in [14] and [6] the reaction-advection equations are eliminated and the system reduced to a system of delay differential equations for the numbers of birds on the patches, where the delays represent the flight times between patches. One of the patches is the breeding patch, and in our models births are assumed to occur only on this patch. Another specific patch is considered to be the winter feeding patch, which is at the other end of the migration flyway. In [14] the spring and fall migration routes linking the breeding and wintering feeding grounds share their stopovers, while they are essentially distinct in [6]. A significant advantage of our approaches to modelling the migration phenomenon is that it recognizes that birds will encounter very different conditions at different stages of their migration. For example, our models easily allow mean flight velocities and in-flight mortalities to vary from sector to sector, and conditions in each patch (such as vulnerability to predation during a stopover) can be different. The aim of the present paper is to develop our previous ideas to understand the role of the interplay between migratory birds and nonmigratory birds, particularly poultry, in the persistence and recurrence of H5N1 in endemic regions.
In this paper we focus on the role of migratory birds in H5N1 spread, and we concentrate on the interaction of one migratory species with poultry. The poultry do not move, and no explicit consideration is given here to the trade in poultry and its vaccination. Unlike our previous works, in which we allowed any number of stopover patches, here we consider four patches which are the breeding and wintering patches together with two other patches which are stopover patches, one for the outward and one for the return migration. Migratory birds in the air are accounted for indirectly. We allow for the possible presence of poultry at all four patches, although in reality there may be none at one or more of them.
In section 2, we present a careful model derivation along with some baseline assumptions about the initial conditions and birth rate function at the breeding patch. We then discuss the boundedness and positiveness of solutions and the dissipativeness of the model system (section 3); the global threshold dynamics: global asymptotical stability of the disease-free nontrivial equilibrium (section 4); and the persistence of infection (section 5). In section 6, we discuss how our results and model can be extended to more biologically realistic situations involving periodic coefficients, and we present numerical simulations using available literature data on the poultry around the region of Poyang Lake in China, where live H5N1 infected ducks were found and H5N1 is thought to be endemic [10]. In this region, poultry farming includes backyard chickens and domestic ducks and geese [10,33].

Model derivation.
Our model is basically of patch type but is derived with the aid of reaction-advection equations which describe the migration of the birds along the flyways. There are four patches on which poultry reside and migratory birds pass through. The four patches are labelled b, o, w, r for the breeding patch, the outgoing stopover patch, the wintering patch, and the returning stopover patch, respectively. The breeding patch is that on which migratory birds breed, while the wintering patch is where they spend the winter season. Patches o and r are stopover patches along the fall and spring migratory routes, respectively. In reality migratory birds usually have multiple stopovers on both their fall and spring routes; these have been lumped together so that there is just one stopover in each direction.
We let S b m and I b m denote the numbers of susceptible and infected migratory birds in the breeding patch b; S o m and I o m the numbers in the outgoing stopover patch o; and so on. Similarly, I b p , I o p , I w p , I r p denote the numbers of infected poultry in the four patches. We shall assume that the total number of poultry (susceptible plus infected) in a patch is a constant N p with superscript to denote the patch. No such assumption is made for migratory birds.
For the model derivation, proper account needs to be taken of the birds that are in flight between two patches on their migratory route. We shall allow mean flight velocities and in-flight mortalities to depend on susceptible/infected status and on the sector. Let S bo m (t) denote the number of susceptible migratory birds in flight between patch b and patch o, and let x denote the physical space. The flyway is to be thought of as a one-dimensional closed curve starting and ending at the location of patch b, denoted x b . Therefore, x can be formally defined as arc length along the flyway for a particular sector. We illustrate the model derivation by considering susceptible migratory birds along one particular sector, and we shall focus on the sector from the breeding patch b to the outward stopover patch o. The other sectors can be treated similarly. With x denoting the arc length measured from patch b located at x b , the quantity of susceptible birds in the air along points of this flyway will be modelled by a density (with respect to the arc length) s bo m (t, x), which satisfies the linear reaction-advection equation where U bo s is the mean flight velocity of susceptible migratory birds between patches b and o, and μ bo ms is the per capita in-flight mortality for these birds as they transit between those patches. The flux of birds out of patch b is U bo s s bo m (t, x b ), and we want to use (2.1) to calculate the flux into patch o, which is U bo s s bo A calculation using (2.1) shows that Integrating from x b to x o and then choosing ξ . What we have achieved here is to rigorously derive the first term in the right-hand side of the third equation of system (2.4) below. All other flux terms into the various patches involve other time delays τ and journey survival probabilities α with appropriate sub-and superscripts, which we trust are self-explanatory. These terms can be derived similarly. The remaining terms in system (2.4) involve per capita mortalities within patches (for example, μ b ms for susceptible migratory birds in patch b) and mass action terms for the rates at which susceptible birds become infected. A susceptible migratory bird on a patch may catch the virus either from an infected migratory bird, or from infected poultry on the patch. For patch b the respective contact rates are denoted β b m and β b pm . We assume that no migratory bird can catch the disease while in flight between patches.
As far as birth of migratory birds is concerned, we assume that this happens only on the breeding patch b and that only susceptible birds can produce offspring. The birth rate is therefore taken to be of the form B m (S b m ), where the function B m (·) satisfies certain assumptions to be stated later. Based on all these assumptions, we propose the following system of delay differential equations as a model for the migratory bird dynamics: The mean time spent in a particular patch is relevant to the flowout rate of the bird species from the patch and can be written, for example, in the case of the outward stopover patch, as 1/d s ow for susceptible birds. We aim to keep all parameters free where possible, but it would be reasonable to state that the mean time spent in a particular patch is probably related to other system parameters. For example, in the case of the outward stopover patch it will very probably be related to the duration of the flight to that patch, since this will determine how much refueling is needed, and possibly also to the anticipated fuel requirement for the next segment of the journey if birds are behaving optimally, since flight costs increase considerably with increasing body mass. At stopover sites there is also a tradeoff between gaining fuel and avoiding predation. These considerations are discussed in Weber, Ens, and Houston [34]. Birds may also spend more time at sites where energy expenditure would be low, and this in turn depends on the temperature at the site. Bauer et al. [2] cite evidence that a temperature decline of 10 • C increases energy expenditure by around 40%. The equations of system (2.4) are coupled with the following four equations governing the poultry population dynamics (2.5). The poultry do not migrate, and we defer the consideration of commercial trading of poultry to a future study. We assume that on each patch the total poultry population is constant, so, for example, p is the number of susceptible poultry on that patch. On a particular patch, poultry can catch the virus either from infected poultry on that patch or from infected migratory birds that happen to be on that patch. The respective contact rates are β p and β mp with superscripts to denote the patch. Per capita mortality for poultry (including the culling of the poultry as a disease control and prevention measure) is denoted by μ p with superscript. These considerations lead to the equationsİ (2.5) All model parameters are assumed to be positive. The appropriate initial conditions for system (2.4)-(2.5) are which, as the product of the Banach spaces of continuous functions in certain initial intervals (for migratory birds) and the Euclidean spaces (for poultry), is a Banach space. In The biologically realistic set of initial data is Y , consisting of the elements of X which satisfy condition (2.6). We have Y = Y 0 ∪ ∂Y 0 , where Y 0 is basically the subset of Y corresponding to the presence of infected poultry or infected migratory birds somewhere, so that is ≥ 0 and ≡ 0 on its domain; at least one of I b p , I o p , I w p , I r p > 0; and/or at least one of I b m (·), I o m (·), I w m (·), or I r m (·) is ≥ 0 and ≡ 0 on its domain}, (2.9) where the domains for the various functions are given by (2.8). Then, ∂Y 0 is a subset of Y relating to those solutions of (2.4)-(2.5) for which all the I m and I p variables remain zero for all time, i.e., the disease-free solutions.
We now discuss the assumptions on the birth function B m (S b m ). First note that, since the α quantities are bounded by 1, We shall often assume that the birth function satisfies the following hypotheses.
Hypothesis (H1) is a natural assumption to ensure the existence of a positive equilibrium of the isolated bird species population in the absence of diseases. This assumption also implies, in the case where B m is continuously differentiable, that For the linearized stability of equilibria, and also for persistence, we need to strengthen this to strict inequality in the following hypothesis.

Positivity and boundedness.
We first establish a result on nonnegativity and eventual strict positivity. The birth function here does not need to satisfy all of the assumptions in (H1). (not necessarily the same) patch the number of infected or susceptible migratory birds is not identically zero on the relevant initial interval holds, then there exists T * > 0, which is independent of the initial data, such that I m (t) > 0 and I p (t) > 0 for all t > T * and on every compartment.
Proof. The first statement of the theorem follows immediately from Smith [ . Indeed, if this were false, then nonnegativity of solution components would yield that the first and third terms in the right-hand side of the fourth equation (the first term in particular) are identically zero for t ∈ [0, τ i bo ], and this contradicts A simple comparison argument therefore assures us that I o m (t) > 0 for all t ≥ t * . Having established that the variable I o m (t) becomes positive and remains so, we deduce the same for the variable I w m (t) using the sixth equation of (2.4). This equation shows that I w m (t) cannot remain zero; otherwise the same would be true for I o m (t). Once I w m (t) becomes positive it remains so. Similarly, we conclude the same for I r m (t) and I b m (t). Theorem 5.2.1 of [30, p. 81] also implies that the number of infected poultry on a patch cannot exceed N p for that patch. To see that the numbers of infected poultry become and remain positive on each patch is straightforward; for example, if I b p (t) remained zero, then the first equation of (2.5) shows that the same is true for I b m (t), a contradiction. Once I b p (t) has become positive, then, sinceİ b Finally, we prove case (ii) of the second part of the theorem. Suppose it is on patch b that the initial number of infected poultry is positive, and suppose there is another patch on which the number of susceptible migratory birds is not identically zero on the relevant initial interval (if this is the case for infected migratory birds, then case (i) would apply). Arguments similar to those described for case (i) show that I b p (t) > 0 for all t > 0 and also that, whichever patch has its susceptible migratory birds not identically zero initially, eventually all patches must have a strictly positive number of susceptible migratory birds. In particular, , we see that once I b m (t) becomes positive it remains so for all time and in particular on an interval of length τ i bo . But this puts us back into case (i), since we can translate forward in time. So from this point the argument proceeds as before.
We now consider the boundedness of the solutions of system (2.4)-(2.5).
Then every solution of (2.4)-(2.5) with initial data in Y is bounded, and the solution semiflow defined by solutions of (2.4)-(2.5) in space Y is dissipative and has a global attractor.
Proof. The total number of migratory birds N tot m (t) in the system at a given time is the sum of the total number of susceptible and infected migratory birds in the four stopovers b, o, w, and r and the total number of susceptible and infected migratory birds in each of the four flyways linking the stopovers, denoted bo, ow, wr, and ro, i.e., N tot In a given flyway, the number of birds in the air is given by the integral of the bird density. Recall that the density of susceptible birds in each of the four flyways is denoted s bo m , s ow m , s wr m , s rb m , respectively. Similarly, the density of infected birds in each of the four flyways is denoted i bo m , i ow m , i wr m , i rb m . Hence, the total number of birds in the flyway linking patches o and w is and there are similar expressions for the other flyways.
Assuming that the migration flyways are fixed over time, and that the average bird velocity in the flyways is constant, leads toṄ tot Using (2.4) and advection equations analogous to (2.1) in all the flyways leads tȯ Taking the example of infected birds in the flyway bo, each integral term in (3.3) can be expressed in the following form: Recall that the outgoing flux of infected birds leaving patch b is also equal to d i bo I b m (t) and that the incoming flux of birds entering patch o was shown to be equal to in the previous section. So all terms in (3.3) that represent fluxes cancel out in pairs. Hence, (3.3) can be further simplified tȯ This gives From (3.6), it follows that which implies the boundedness of the total number of migratory birds in the system. As the total number of poultry in the system is assumed to be constant on each patch, we obtain the boundedness of all solutions of (2.4)-(2.5) subject to condition (2.6), and the dissipativeness of the associated semiflow on the space. This dissipativeness then implies the existence of the global attractor (see, for example, [15]).

4.
Local and global stability of the disease-free equilibrium. If (H1) holds, then system (2.4)-(2.5) has a trivial equilibrium in which all components of the system are zero, and a disease-free equilibrium in which there are no infected poultry or migratory birds, but there is a population of susceptible migratory birds on each patch. In this equilibrium, the components (S b m , S o m , S w m , S r m ) are determined by the following equations: Combining these equations shows that the equilibrium values of S b m are determined from the following equation: Obviously, S b m = 0 satisfies this equation. Hypothesis (H1) assures us of the existence of one more root S b m =Ŝ b m > 0 of (4.2). The corresponding equilibrium values of the susceptible migratory birds on the other patches are denotedŜ o m ,Ŝ w m , andŜ r m ; these values are found from the last three equations of (4.1).
Our main result of this section is Theorem 4.2, which provides a set of conditions which ensure that the disease-free equilibrium is globally asymptotically stable to perturbations involving the introduction of both infected poultry and infected migratory birds. Before proceeding, we present short subsections which deal with subsystems of either migratory birds or poultry in the absence of disease.

Subsystem of migratory birds.
It is necessary for what follows to understand completely the dynamics of the subsystem of four equations consisting of just the equations for the susceptible migratory birds, when the numbers of infected migratory birds and infected poultry remain zero for all time on every patch. This and if (H1) holds, then it has exactly two equilibria, (0, 0, 0, 0) and , the latter being found by solving (4.1). The initial data for system (4.3) is just the S m components of (2.6), and each of the four initial functions is assumed to be nonnegative. for which at least one of the initial components is not identically zero on its initial interval (for example, Similarly to the proof of Theorem 3.1, we see that all variables in system (4.3) remain nonnegative for t ≥ 0 and that in fact they all eventually become strictly positive and remain so. The structure of system (4.3) allows us to apply Theorem 5.4.1 of Smith [30] on convergence to equilibrium, to conclude that all solutions of (4.3) must converge to one of the equilibria. To decide which one is the attractor, we can use linear stability theory. The linearization of (4.3) about the equilibrium r has a structure involving positive coefficients for all delayed variables, so that it suffices to consider only the real roots of the characteristic equation. Solutions of the linearized system with temporal dependence exp(λt) lead to the characteristic equation The left-hand side of (4.4) is increasing in λ, while the right-hand side is decreasing in λ when λ ≥ 0. Therefore, if the left-hand side exceeds (is smaller than, respectively) the right-hand side when λ = 0, we are assured that all real roots including the dominant root are negative (or positive, respectively). Therefore, (H2) ensures that the equilibrium (0, 0, 0, 0) is linearly unstable and the equilibrium ( is the global attractor, we need to show that no solution satisfying the requirements on the initial conditions may approach the unstable equilibrium (0, 0, 0, 0) as t → ∞. This unstable equilibrium still has an infinite-dimensional stable manifold due to the presence of delay. We briefly explain why there cannot be a solution (S b m (t), S o m (t), S w m (t), S r m (t)) → (0, 0, 0, 0). Another similar calculation is explained much more thoroughly later, in the proof of our Theorem 5.2. Since each component of the solution becomes and remains strictly positive, without loss of generality we can assume each component has this property initially. A contradiction can be reached by using a comparison argument to show that the solution which supposedly approaches the zero solution is also bounded below by the exponentially growing function δ exp(λt)d, where δ > 0 is a suitably chosen small constant, λ > 0 is the dominant eigenvalue of the linearization at (0, 0, 0, 0), and exp(λt)d is the solution of the linearized system corresponding to the dominant eigenvalue. The vector d is known to have positive components [30, Theorem 5.5.1].

Stability with migratory birds absent.
We next consider the stability of the disease-free equilibrium in the situation when migratory birds are completely absent. Then we need only consider (2.5) when the I m variables are all zero. But in this case system (2.5) decouples completely, and each equation determines the number of infected poultry on a particular patch. It follows immediately that the conditions for the infected poultry to die out on every patch when no migratory birds are introduced are

Stability (of the positive disease-free equilibrium) with poultry absent.
If there are no poultry (susceptible or infected), then we can calculate necessary and sufficient conditions for the disease-free equilibrium to be linearly stable to perturbations involving the introduction of migratory birds. The subsystem that needs to be studied is the linear systeṁ (4.6) The structure of this system enables us to apply the theory due to Smith [29]. Corollary 3.2 of that paper assures us that the zero solution of system (4.6) is asymptotically stable if and only if the zero solution of the corresponding undelayed system is asymptotically stable. It needs to be carefully noted, however, that the corresponding undelayed system is system (4.6) with the delays set to zero where they appear in the arguments of the state variables. Delays also feature indirectly in system (4.6) via the α coefficients such as α i rb . Those delays are not set to zero. Corollary 3.2 of [29] also provides a useful set of necessary and sufficient conditions, involving determinants, for the stability modulus (the eigenvalue of largest real part, which is known to be real here) to be strictly negative. Applying these conditions, noting carefully what constitutes the undelayed system here, leads to the following conditions for the asymptotic stability of the nontrivial disease-free equilibrium when there are no poultry: [29] actually provides just four inequalities and, when applied to (4.6), furnishes three of the first four inequalities shown above, together with the fifth. However, it can be shown that these inequalities hold if and only if all five inequalities of (4.7) hold.

Stability with infection in both poultry and migratory birds.
The next theorem presents conditions for global asymptotic stability of the disease-free equilibrium for solutions involving both infected poultry and infected migratory birds. Not surprisingly, one needs to impose an additional inequality involving the crossinfection parameters β c pm , with c = b, o, w, r, as follows: In condition (4.8), we implicitly assume that each factor on the left-hand side is positive. This is equivalent to the requirement Since the I m and I p variables remain nonnegative, the four S m variables satisfy a system of differential inequalities. This system is precisely (4.3) with every = sign replaced by ≤. A comparison argument using Theorem 5.1.1 of Smith [30] implies that each S m component is bounded above by the corresponding S m component of the solution of (4.3), where the two components share the same initial data. In view of Theorem 4.1, it follows that, for each c = b, o, w, r, Therefore, for any > 0, it is true that, for sufficiently large t, S c m (t) ≤Ŝ c m + for each c = b, o, w, r. Therefore, for sufficiently large t, the I m and I p variables obey the inequalitieṡ (4.10) Again by Theorem 5.1.1 of Smith [30], each I m and I p variable is bounded above by a system of differential equations associated with (4.10), namely, system (4.10) with every ≤ replaced by =. We need to show that every I m (t) and I p (t) tends to 0 as t → ∞, and it suffices to show that this is the case for the corresponding system of differential equations, which is of course linear. Purely for economy of notation, we shall in fact finish the proof for the case when = 0. The proof has to work for a small positive , but the usual continuity arguments L. BOUROUIBA, S. A. GOURLEY, R. LIU, AND J. WU assure us of this. We shall briefly comment on this point again later. The fact that inequality (4.8) is a strict inequality is crucial.
If we seek trial solutions of the differential equation system associated with (4.10) of the form , we obtain a characteristic equation of the form G(λ) = det (λI − A(λ)) = 0, where A(λ) is a matrix whose off-diagonal entries are nonnegative and nonincreasing in λ. We may apply Theorem 5.5.1 of [30, p. 92] to conclude that the linear stability of the equilibrium is completely determined by the real roots of the characteristic equation. The algebra is very complicated, and it is necessary to make full use of the relationships between the components ξ b , η b , . . . to simplify the determinant as much as possible. Eventually, we can show that G(λ) is given by (4.11) We shall prove that all real roots of the equation G(λ) = 0 are strictly negative. This will be achieved by using a monotonicity argument, but, before proceeding, we note that we are in the situation in which Corollary 5.5.2 of [30, p. 93] applies, which assures us that the disease-free equilibrium is linearly asymptotically stable if and only if it is linearly asymptotically stable in the corresponding system without delays in the arguments of the variables (the model coefficients still involve delays through the α coefficients). Removal of the delays only in the arguments of the variables, leaving their presence via the α coefficients undisturbed, corresponds to formally setting the exponential term in (4.11) to 1 and making no other changes. We write the modified characteristic equation as G 0 (λ) = 0. This equation is no longer transcendental but an eighth order polynomial. After some algebra, the equation G 0 (λ) = 0 can be rewritten as and (4.13) with (4.14) By inequalities (4.5) and (4.7), all quantities in (4.14) are positive. Therefore F 2 (λ) is decreasing for λ ≥ 0. Also, by (4.9), each factor in the expression for F 1 (λ) is monotonically increasing and positive when λ ≥ 0. Therefore, F 1 (λ) is monotonically increasing for λ ≥ 0. Moreover, by inequality (4.8), F 1 (0) > F 2 (0). By continuity, this strict inequality will still hold if we repeat the analysis for a sufficiently small positive in (4.10). For such an it is therefore clear that all real roots including the dominant root of F 1 (λ) = F 2 (λ), and therefore of G 0 (λ) = 0, are strictly negative.
We have shown via a comparison argument involving an associated linear system that the I m and I p variables all approach zero as t → ∞. The asymptotic behavior of the S m variables is then determined by subsystem (4.3), the dynamics of which is described in Theorem 4.1. The proof is complete.

Disease persistence.
In this section we prove that, under certain conditions, all the variables I m (t), I p (t) representing the infected compartments are persistent in the strong, uniform sense, i.e., that there exists some constant η > 0 which is independent of the initial conditions such that, for each c = b, o, w, r, This is achieved using the persistence theory due to Hale and Waltman [15] and, in particular, Theorem 4.1 of their paper. There are other approaches to persistence; see, for example, Samanta [26] for a nonautonomous epidemic model. Let a metric space Y be the closure of an open set Y 0 , so that Y = Y 0 ∪ ∂Y 0 , where ∂Y 0 is the boundary of Y 0 . Let T (t) be a C 0 -semigroup on Y satisfying Assume that the restricted semiflow has the global attractor A ∂ , and assume that where ω(x) is the ω-limit set of x. Then we have the following result from Hale and Waltman [15].

Theorem 5.1 (Hale and Waltman). Suppose that T (t) satisfies (5.2) and that
Hale and Waltman permit the M i to be chosen as isolated invariant sets and not necessarily as equilibria, but for our purposes here each M i is a boundary equilibrium of our model. Assuming (H1) holds, these are the two equilibria, denoted M 1 and M 2 , in which the variables I m and I p are zero on every patch. The set M is the union M = M 1 ∪ M 2 . In our setting, acyclicity of the covering M 1 ∪ M 2 means that neither of the equilibria M 1 , M 2 is connected to itself via a homoclinic connection within ∂Y 0 , and that there is no cycle with is the stable manifold of M i . Note that the superscript s on any letter other than W has a different meaning (susceptible), but there is no risk of confusion.
In this paper Y and Y 0 are defined in the paragraph containing (2.9). Next we define the C 0 -semigroup T (t). Recall the initial data has the form (2.6). For any initial datum

) be the solution of system (2.4)-(2.5) and define T (t) by
where again the domains for the functions are inferred from (2.8); for example, We prove the following theorem on disease persistence.

(5.5)
Assume further that the initial data (2.6) satisfies (i) or (ii) of Theorem 3.1, in either case with the additional requirement that at least one of the S m components not be identically zero on its initial interval. Then there exists η > 0, which is independent of the initial conditions, such that the solution of (2.4)-(2.5), subject to (2.6), satisfies (5.1) for each c = b, o, w, r.
Proof. We will apply Theorem 5.1, and we first seek to prove that T (t) satisfies (5.2). Let us prove that T (t) : Y 0 → Y 0 . Choose some initial datum in Y 0 . We know from Theorem 3.1 that all solution components must remain nonnegative, but suppose that during the evolution there exists a time t * > 0 when the solution leaves Y 0 and arrives at ∂Y 0 . At the time t * we will have Since t * can serve as the initial time, it follows from uniqueness of solutions that, for all t > t * , But this is impossible because if the initial datum lies in Y 0 , then, by Theorem 3.1, every I m and every I p variable eventually becomes strictly positive and remains so. We now prove that We include the full details of only the latter calculation. Suppose, for a contradiction, that W s (M 2 ) ∩ Y 0 = Φ. Then a solution exists which is initially in Y 0 and in W s (M 2 ). The former implies that this solution initially has one of its I p variables positive or one of its I m variables not identically zero, and Theorem 3.1 then yields that, eventually, I m (t) > 0 and I p (t) > 0 for this solution, for every superscript b, o, w, r. We may translate forward in time and assume this to be so initially. Also, since this solution starts in W s (M 2 ), I m (t) → 0 and I p (t) → 0 as t → ∞, for each superscript b, o, w, r. The asymptotic behavior of the S m variables is then determined by the subsystem This is a linear differential inequality system. All delay terms have positive coefficients. By further shrinking if necessary, we can arrange it so that all off-diagonal terms also have positive coefficients. We may then assert, by Theorem 5.1.1 of Smith [30], that each component of a solution of (5.8) is greater than or equal to the corresponding component of a solution of the system of differential equations associated with (5.8) when ≥ is replaced by =, provided that the two solutions are initially ordered in this way.
Calculations similar to those used in the proof of Theorem 4.2 can be used to study the characteristic equation of the system of differential equations associated with (5.8). The characteristic equation can be put in the form F 1 (λ) = F 2 (λ), where F 1 (λ) and F 2 (λ) have the same structure as (4.12) and (4.13), respectively, but involve the small number . Inequality (5.7) implies that F 1 (0) < F 2 (0). Since F 2 (∞) = 0 and F 1 (∞) = 1, it follows that the characteristic equation has a real positive root λ * , with a corresponding solution to the system of differential equations of the form δ exp(λ * t) c for any real δ (but we take δ > 0), in which (crucially) every component of the constant vector c is strictly positive. Note that F 1 (λ) and F 2 (λ) do not need to have the monotonicity properties of F 1 (λ) and F 2 (λ).
We have a solution that supposedly lies in W s (M 2 ) and should approach M 2 as t → ∞. Yet, by the above-mentioned comparison argument, its I m and I p components satisfy, for a suitably chosen δ > 0, This is clearly a contradiction, and so W s (M 2 ) ∩ X 0 = Φ. The positive number δ is chosen so that (5.9) holds initially, which for some variables means at all times in their initial intervals. For example, one requirement is that δ should satisfy , and since we noted earlier that we may translate forward in time and assume that I b m0 (θ) > 0 for all θ ∈ [−τ i bo , 0], this and all other requirements are clearly met if δ > 0 is sufficiently small.
We briefly discuss the proof that W s (M 1 )∩Y 0 = Φ, which is similar. If a solution exists starting in Y 0 and in W s (M 1 ), then the former implies that this solution has one of its S m components starting off not identically zero. Arguments similar to those of the proof of Theorem 3.1 show that, for the full system (2.4)-(2.5), each S m variable eventually becomes and remains strictly positive. Without loss of generality, it may be assumed that each S m variable is strictly positive at all points of its initial interval, and this makes it possible to choose the positive numberδ mentioned below. Since this solution starts in W s (M 1 ), all S m , I m , and I p components tend to zero as t → ∞. For any > 0 it follows that, for t sufficiently large, with similar inequalities for S 0 m , S w m , and S r m . By another comparison argument, each S m variable is bounded below by the corresponding component of the solution of the associated differential equation system. However, for a sufficiently small > 0, the latter has a real positive dominant eigenvalue λ * * -this follows from the inequality , a consequence of (H2)-and corresponding solutionδ exp(λ * * t)c, wherec is a fourdimensional vector with strictly positive components andδ > 0 is chosen advantageously. The rest of the argument proceeds as for the verification that W s (M 2 )∩X 0 = Φ. Theorem 5.1 now yields persistence of the disease, but we still need to establish (5.1), i.e., that it persists in both the poultry and migratory bird populations. We briefly sketch a proof that this is so. If the disease were to die out in the migratory birds, then, by (4.5), it would die out in the poultry also. On the other hand, if it were to die out in the poultry, then, by (2.5), it would die out in the migratory birds. Both possibilities contradict the persistence just proved.

Discussions: Numerical simulations and extensions to including seasonality.
6.1. Threshold and generic conditions. We have established the threshold dynamics of the model system under a generic set of threshold conditions (4.5), (4.7), and (5.5). This set of generic conditions is sharp in the sense that either the disease extinction or the disease persistence occurs generically, and here "generic" means strict inequalities. As an illustration, we carried out numerical simulations using parameter L. BOUROUIBA, S. A. GOURLEY, R. LIU, AND J. WU Table 1 The parameter definitions and values estimated from sources including Bourouiba et al. [6], Gourley, Liu, and Wu [14], Javed et al. [18], and Prins and van Wieren [24].

Para
Meaning Value values available in the literature [6,14,18,24]. The parameters are defined in Table 1 Figure 1 shows that if only migratory birds are considered in all patches and no diseased migratory bird is introduced, then every nontrivial solution of system (4.3) evolves to the positive equilibrium under Hypothesis (H1). If in addition inequality (4.7) holds, then, as Figure 2 shows, the disease-free nontrivial equilibrium of the system (again, with poultry absent) is stable and the disease dies out.
When the poultry is considered, Figure 3 illustrates a situation where (4.7) is not satisfied (the last of the five inequalities does not hold), and in this case the disease-free nontrivial equilibrium loses its stability and there is a disease outbreak. See also Figures 4 and 5 to see how Hypotheses (H1), (H2) and inequalities (4.5), (4.7), and (4.8) combined decide the outcome of H5N1 infection: either the diseasefree nontrivial equilibrium of system (2.4)-(2.5) is globally asymptotically stable or the disease persists.

Seasonality and recurring outbreaks.
The Poyang Lake region in Jiangxi Province, China, is one of the largest areas of fresh water in China, making it one of the largest wintering sites for waterbirds in Asia, attracting tens of thousands of migratory birds of more than 300 species [19]. The Poyang Lake National Nature  Reserve was created in 1983 for migratory birds wintering in this region [19]. It is also one of the densest areas of poultry farming in China, surrounded by about 10,000 square kilometers of farmland and more than 10 million people [28]. Poultry farming typically includes backyard chickens and domestic ducks and geese [10,33]. Densities of poultry around Poyang Lake can be as high as or higher than 100-250 heads of poultry per square kilometer (2007 data on avian influenza from [9]). It is in the Poyang Lake region that Chen et al. [10] identified H5N1 infected ducks, which were labeled as "migratory." This finding raised the question of the role of migratory birds as silent spreaders of H5N1. The lack of proper identification of subspecies and location of capture of the ducks were pointed out on several occasions as limitations to the conclusion that the birds were indeed wild or migratory [12]. For example, the migratory birds sampled around the lake were said to involve falcated teal, spotbill, and mallard ducks; however, only the first is surely migratory, the second is a common breeder around the lake, and mallards and their descendants are the most common domesticated species released in the vicinity of the lake by local farmers [12]. The role of migratory wild birds in spreading and maintaining the endemicity of H5N1 is still to some extent an open question, particularly due to conflicting evidence of pathogenicity of H5N1 in domesticated and wild ducks [13]. In order to account for this complex interplay between poultry and migratory birds in an area as dense as Poyang Lake, we now apply the simplified model of four patches to migratory birds wintering in the Poyang Lake region, where H5N1 is assumed to be endemic among poultry. We proceed by first accounting for the seasonality of the migration, which can be modelled by modifying the basic model (2.4) to allow temporal periodicity as follows: Here, the time variable t is explicitly incorporated into the birth rate B m (S b m , t) at the breeding patch, the in-flight survival probabilities α s c (t), α i c (t) (which are no longer given by expressions of the form (2.3)), and the migration rates d s c (t) and d i c (t) between patches (c = rb, bo, ow, wr, respectively). For annual migration, we assume these functions are all T -periodic in t. For now, we assume the migration functions are all positive and continuous, as this will substantially simplify the discussions below. Under the assumption that B m (S b m , t) is a nonnegative continuous and uniformly bounded function for all t ∈ R and S b m ≥ 0 (with B m (0, t) = 0 for all t ∈ R), we can use arguments similar to those in section 3 to establish the nonnegativity and boundedness of solutions of system (6.1) coupled with (2.5) and subject to the initial conditions from Y 0 satisfying (2.6).
As for the existence and global attractivity of a positive T -periodic solution of the subsystem of migratory birds in the absence of disease infection, namely, we need to assume that the birth rate function B m (S b m , t) is continuously differentiable with respect to S b m , and we need to consider the spectral radius r(T M ) of the time-T operator T M of the linearized system at the trivial equilibrium. If the birth rate function B m (S b m , t) is sublinear in the sense that B m (λS b m , t) > λB m (S b m , t) for all λ ∈ (0, 1), S b m > 0, t ∈ R, we can conclude that when r(T M ) < 1 all solutions of system (6.2) converge to zero, and when r(T M ) > 1 system (6.2) has a unique positive T -periodic solution (Ŝ b m (t),Ŝ o m (t),Ŝ w m (t),Ŝ r m (t)) and every solution as specified above converges to this periodic solution. This threshold dynamics for subsystem (6.2) can be obtained using an argument similar to that of Theorem 3.2 in [14] based on a general threshold dynamics theorem for order-preserving maps (Theorem 2.3.4 of [35, p. 48]).
The argument for Theorem 4.2 can be adapted to establish the global attractivity of the disease-free periodic solution (Ŝ b m (t), 0,Ŝ o m (t), 0,Ŝ w m (t), 0,Ŝ r m (t), 0, 0, 0, 0, 0) by considering the linear periodic delay systeṁ of healthy birds (μ ms = 2.2 × 10 −3 day −1 ). Concerning the H5N1 infected teals, we refer to various sources of inoculation studies of wild birds (e.g., [8,20]) combined, for example, in Bourouiba, Teslya, and Wu [5] for a duck species of small size. In Bourouiba, Teslya, and Wu [5] a susceptible-asymptotic-infected-recovered (SAIR) epidemic model modelling the highly pathogenic avian influenza (HPAI) infection of wood ducks was considered with two infectious stages of increasing associated force of infection. These were associated with a mean duration of infectious period of 9.86 days, transmission parameter in the range β m = 1.5-3 × 10 −3 day −1 , and death rate of infected migratory birds of 0.133 day −1 [5,8], leading to a removal rate of infectious migratory birds of μ mi = 0.23 day −1 , which will be used in the present simulations. On the other hand, the study of Roche et al. [25] focusing on the persistence of avian influenza in wildlife used transmission parameters in the range β m = 5 × 10 −5 -0.1 year −1 . We used values of β m in the range that is discussed in the previous literature. The densities of poultry in the area where common teals were surveyed is considered to be 100 to 250 poultry per square kilometer. When considering the bird protection area of Poyang Lake, which measures 224 km 2 , these densities correspond to 22,400 to 56,000 poultry on the wintering ground. The density of poultry along the migratory route (passing by the region of Shenyang and North Korea) remains comparable to that surrounding Poyang Lake; hence we consider similar numbers of poultry on the migratory route. However, on the northernmost breeding areas (teals and wigeons) in Russia and Siberia there was clearly a drop of density of poultry [9,32]; hence, no poultry is considered on the breeding ground. Concerning the disease dynamics of HPAI, Bouma et al. [4] estimated various transmission parameters for poultry. These were determined from inoculation experiments using HPAI H5N1. Pairs of chickens, one directly inoculated with HPAI H5N1, and the other infected by H5N1 through contact, were examined. The estimated transmission parameter of HPAI in poultry is 4.78 × 10 −4 -0.4 day −1 [4,17]. The mean infectious period of contact-infected birds of μ p = 1/2.5 day −1 [4]. Note that the cross-species transmission parameters β pm and β mp are assumed to be smaller than those characterizing transmission within a species. We assume that β mp = β pm = β m /10. Figures 6 and 7 show the dynamics of the migratory bird population in the absence of poultry. In the absence of avian influenza, the solutions of the subsystem of migratory birds converge to a positive solution. In the presence of avian influenza introduced by migratory birds, the system either evolves to a stable disease-free periodic solution ( Figure 6) or the migratory bird population becomes extinct, which happens, for example, when β m = 1.5 × 10 −3 day −1 (not shown). In contrast to the situation in [14], nonperiodic disease persistent solutions are also observed for transmission coefficients moderately large; however, this particular regime is particularly sensitive to changes in parameter values (Figure 7). Figure 8 shows the convergence of the migratory bird subsystem to a periodic solution with disease persistence induced by the endemicity of the disease in poultry in the wintering patch. We find that the periodicity of the number of infected poultry as migratory birds return during the year is possible. In turn, it is the endemicity of poultry on the wintering patches that is sustaining the epidemic in the migratory bird population.

Environmental contamination and other factors.
In this paper we use only direct transmission models. Breban et al. [7] examine the additional role of environmental transmission, i.e., virions in the environment that have been shed by infectious birds. Although environmental transmission rates are hundreds of times lower than direct transmission rates, it is known that the virions can persist in the environment for a long time. Indeed, environmental transmission apparently provides a possible persistence mechanism in situations where an epidemic would not be sustained by direct transmission alone [7,22]. The model setting in our study here does not incorporate the environmental contamination and thus may underestimate the likelihood of an outbreak for a given set of parameter values. How to extend our model and analysis to address the role of environmental transmission in avian influenza spread would be an interesting challenge for future study.

Poultry trading.
Iwami et al. [17] considered the effect of the illegal trade of poultry on the efficiency of the control of avian influenza outbreaks. In their model, only poultry was considered and the trade was accounted for with movement of poultry between two patches characterized by two different control configurations showing that only the complete eradication in the vaccinated area can lead to the complete eradication in the connected other area.