Spatiotemporal Distributions of Migratory Birds: Patchy Models with Delay

We derive and analyze a mathematical model for the spatiotemporal distribution of a migratory bird species. The birds have specific sites for breeding and winter feeding, and usually several stopover sites along the migration route, and therefore a patch model is the natural choice. However, we also model the journeys of the birds along the flyways, and this is achieved using a continuous space model of reaction-advection type. In this way proper account is taken of flight times and in-flight mortalities which may vary from sector to sector, and this information is featured in the ordinary differential equations for the populations on the patches through the values of the time delays and the model coefficients. The seasonality of the phenomenon is accommodated by having periodic migration and birth rates. The central result of the paper is a very general theorem on the threshold dynamics, obtained using recent results on discrete monotone dynamical systems, for birth functions which are subhomogeneous. For such functions, depending on the spectral radius of a certain operator, either there is a globally attracting periodic solution, or the bird population becomes extinct. Evaluation of the spectral radius is difficult, so we also present, for the particular case of just one stopover site on the migration route, a verifiable sufficient condition for extinction or survival in the form of an attractive periodic solution. This threshold is illustrated numerically using data from the U.S. Geological Survey on the bar-headed goose and its migration to India from its main breeding sites around Lake Qinghai and Mongolia.

(2.1) . . . where Similarly, for the return migration, from the winter feeding site x n to the summer breeding site x 1 , the model is The number of birds in the air is taken care of in this model by the time delay terms which model the time of flight between one location and the next. For example, in the outward migration (model (2.1)), birds migrate from location x 1 (the breeding site) at a rate d 1,2 (t)S 1 (t). The (curvilinear) flight distance to x 2 is l 1 , so this stage of the migration takes time l 1 /U 1 = τ 1 , and therefore the rate of arrival of these birds at location x 2 is d 1,2 (t − τ 1 )S 1 (t − τ 1 ) after correction for mortality during the flight. The per capita mortality for this sector of the migration is μ 1,2 (t), and this leads to a survival probability for the sector of α 1,2 (t), as we shall show. These ideas are put onto a more rigorous mathematical footing as follows, but first we need to clarify a mathematically subtle point. Although x is a continuum, each stopover location x i is to be regarded as a patch that can accommodate a finite number of individuals (it is helpful to think of each x i not only as a point on the continuum x but also as a label for a small but finite area adjacent to the continuum at that particular point, perhaps like service stations adjacent to a freeway).
We now derive the expression shown for the first term in the right-hand side of the S 2 equation in system (2.1); this derivation will also be needed later in counting the total number of birds, both within patches and in transits. The first term in the right-hand side of the S 2 equation in system (2.1) is the flux into patch x 2 . Now, along the route between x 1 and x 2 (we write x ∈ (x 1 , x 2 ), but we remind the reader that the route could be curved), the population of birds is described by a density function s(t, x) obeying the advection equation and the flux at any point x ∈ (x 1 , x 2 ) is U 1 s(t, x). Now, the flux at x 1 , which is the rate of leaving patch x 1 , is d 1,2 (t)S 1 (t). The flux at x 2 (the rate of arriving into patch x 2 ) is U 1 s(t, x 2 ), which needs to be computed. If we define Integrating from x 1 to x 2 , and remembering that all distances are curvilinear distances from some reference point on the flight path, so that l 1 = x 2 − x 1 , we obtain Setting ξ = t − x 2 /U 1 and remembering that l 1 /U 1 = τ 1 , we have The flux into patch x 2 is U 1 s(t, x 2 ) and is given by which justifies the first term in the right-hand side of the S 2 equation in (2.1). The modeling of the outward and return migrations (systems (2.1) and (2.2)) can be combined into the following single system, valid at all times: This combined system can be obtained by assuming the functions d i,i+1 (t) and d i+1,i (t) to be zero at the time t when there is no migration from x i to x i+1 or from x i+1 to x i . More precisely, all of these functions are zero at nonmigratory times of year, those of the form d i,i+1 (t) model the outward migration and are therefore zero at all times of the year when there is no outward migration, and those of the form d i+1,i (t) are zero except for the times when return migration takes place. Unless otherwise stated, in the remainder of this paper we shall assume the following holds: continuous functions of their arguments and are periodic with respect to the time variable t with the period T > 0. The functions μ i (t), μ i,i+1 (t), μ i+1,i (t), d i,i+1 (t), and d i+1,i (t) are not identically zero. b(0, t) = 0, b(S 1 , t) ≥ 0 for all S 1 > 0 and t ∈ R, and b(S 1 , t) ≡ 0. System (2.3) can be written in the following form: for t > 0. The initial data are (2.5) It should be noted that the variables S 2 , S 3 , . . . , S n−1 each appear as both undelayed and delayed variables, with precisely two delays involved for each of these variables. This is because i = 2, 3, . . . , n−1 refer to stopover locations and birds may arrive at such locations from either side according to which migration is in progress; for example, the delayed occurrences of S 2 are S 2 (t − τ 2 ) and S 2 (t − τ 1 ). The first of these relates to outward migration from the site x 2 on to x 3 , which has journey time τ 2 , and the second relates to return migration from x 2 to x 1 , which has journey time τ 1 . The delayed occurrences of S 1 and S n each involve just one delay because these are the populations at the breeding and feeding locations x 1 and x n . Migration to or from these sites can be only to or from x 2 and x n−1 , respectively. For the study of system (2.3), with initial conditions (2.5), a natural choice for the state space is . . .
which is a Banach space when equipped with the norm with the I i defined in (2.7). As we illustrate below, it is necessary to introduce the subspace M := φ ∈ Y : φ i (0) = 0 for i = 1, . . . , n; This subspace, determined by the migration patterns, is a closed subspace of Y . It is basically the set of all initial conditions for which the solution of (2.3) remains identically zero for all t > 0. We start with the following result on positivity. Proposition 2.1. For the system (2.3), with initial conditions (2.5), we have that S i (t) ≥ 0 for all t > 0 and for all i = 1, 2, . . . , n.
Proof. The nonnegativity of each S i (t) follows immediately from Theorem 5.2.1 on page 81 of Smith [7]. Next note that, since the delayed terms in the right-hand side of system (2.3) all have nonnegative coefficients, a comparison theorem applies (for example, Theorem 5.1.1 on page 78 of [7]). Using these facts, it is easily seen that each variable remains strictly positive if it ever becomes so. For example, in the case of S 1 , it follows from S 1 (t) ≥ −(d 1,2 (t) + μ 1 (t))S 1 (t).
To complete the proof, it is therefore necessary to establish that there exists t * > 0 such that if (S 0 1 , . . . , S 0 n ) / ∈ M , then at least one component of S = (S 1 , S 2 , . . . , S n ) becomes strictly positive, and that once this happens all other components eventually become strictly positive. Suppose that (S 0 1 , . . . , S 0 n ) / ∈ M . We treat in detail the possibility that there exists In the latter case, the i = 2 equation of (2.3) gives Thus S 2 (t) > 0 for all t > t 1 . Since d 2,3 (t) is periodic and not identically zero, there exists a subsequent time between t 1 + τ 2 and t 1 + τ 2 + T at which d 2,3 (t − τ 2 )S 2 (t − τ 2 ) > 0, and therefore S 3 (t) becomes strictly positive and remains so for t ≥ t 1 + τ 2 + T . We eventually conclude the same for S 4 , S 5 , . . . , S n for t ≥ t 1 + 2T + τ 2 + τ 3 , t ≥ t 1 + 3T + τ 2 + τ 3 + τ 4 , . . . , t ≥ t 1 + (n − 2)T + τ 2 + · · · + τ n−1 . In principle, S 1 (t) could still be zero at this time, but S 1 (t) will also eventually become strictly positive through the presence of the S 2 term in the S 1 equation when return migration is in progress (see the first equation of (2.2)) for t ≥ nT + n i=1 τ i := t * . Other cases can be dealt with in a similar fashion. As the case considered in detail in the above proof shows, one way to ensure that all patches eventually become occupied is to have d 1,2 (θ)S 0 1 (θ) ≡ 0 (but not simply S 0 1 (θ) ≡ 0). This works by ensuring that some individuals are present at location x 1 in the outward migration season. The condition of nonmembership of M is more general because it allows the possibility of initially introducing individuals at other patches. Two solutions of system (2. 3) corresponding to initial data φ, ψ ∈ Y will remain identical to each other for all future times t ≥ 0 if φ − ψ ∈ M . To see this, consider two solutions S = (S 1 , S 2 , . . . , S n ) andS = (S 1 ,S 2 , . . . ,S n ) of (2.3) starting, respectively, from initial data φ, ψ ∈ Y such that φ − ψ ∈ M . Then Z := S −S is initially in M . We claim that Z ≡ 0 for all t > 0. Since all terms of system (2.3) are linear except for the b(S 1 (t), t) term, it is easily seen that Z = (Z 1 , Z 2 , . . . , Z n ) satisfies a system exactly like (2.3) but with the b(S 1 (t), t) term replaced by (∂b(ζ(t), t)/∂S 1 )Z 1 (t), where ζ(t) is a function which arises through an application of the mean value theorem. From the structure of this system and the fact that Z starts in M , it is clear that indeed Z ≡ 0 for all t > 0 as claimed.
It is natural to consider the quotient space X = Y /M consisting of the equivalence classes The importance of considering such a quotient space will become clear in the next section, when a certain notion of strong monotonicity of solutions is required. In particular, we have in The next result establishes the boundedness of solutions. We note that this result does not need the periodicity of the model coefficients.
is the total number of birds, .
Proof. The stated conditions are sufficient for the nonnegativity of each S i (t), t > 0, which is needed for this proof. At times when the outward migration is in progress, the total number of birds S tot (t) is given by where s i (t, x), i = 1, 2, . . . , n − 1, is the airborne density between x i and x i+1 , which for the outward migration is determined by Note that using (2.1). While this analysis is valid only at times when the outward migration is in progress, similar ideas lead us to conclude that, in fact, at all times, with μ given by (2.9), and this leads to (2.8). Note that S tot (t), though always equal to the total number of birds, is given by expression (2.10) only when outward migration is in progress.

Threshold dynamics.
Let X be a real Banach space, and let C ⊂ X be a cone, that is, a nonempty closed subset of X with the properties (i) λC ⊂ C for any nonnegative λ, A cone induces a partial ordering on a Banach space. We say that x ≤ y if and only if y − x ∈ C. We define x < y to mean that x ≤ y and x = y, and we say that x y if and only if y − x ∈ int C, the interior of C.
0 when u > 0). The following version of the Krein-Rutman theorem combines the relevant parts of the two versions to be found on pages 226 and 228 of Deimling [3].
Theorem 3.1 (Krein-Rutman theorem). Let X be a Banach space, C ⊂ X be a total cone, and G : X → X be a compact linear operator that is positive. Suppose that the spectral radius , and there is no other eigenvalue with a positive eigenvector. To state and prove our main result, we recall that a natural phase space for system (2.3) is the Banach space Y . However, the consideration of a solution to eventually become positive for any nontrivial initial data leads to the introduction of the closed subspace M , which yields the quotient space X = Y /M. It is natural to introduce the cone C ⊂ X as the set of the equivalence classes [φ] so that φ ∈ Y and As noted earlier, the solution starting from ψ ∈ [φ] coincides with the solution starting from φ, so we can just talk about the solution starting from a given initial equivalence class. In particular, corresponding to a given initial datum S 0 = (S 0 1 , S 0 2 , . . . , S 0 n ) satisfying (2.5), there is one and only one solution of (2.3) S = (S 1 , S 2 , . . . , S n ) that is defined and nonnegative for all future times t ≥ 0. Furthermore, this is also the unique solution of (2.3) through any initial valueS 0 ∈ [S 0 ]. Therefore, we can define the map S t : In what follows, we will work on a representation of this operator, so we will write is essentially the state of the solution of (2.3) at time t = T corresponding to the use of S 0 as the initial datum. A fixed point of F corresponds to a periodic solution of (2.3) of period T .
We can now state the main threshold dynamics result, as long as the birth function b(S 1 , t) is subhomogeneous in S 1 .
Theorem 3.2. In system (2.3), suppose (H1) holds, b sup < ∞, and λb(S 1 , t) < b(λS 1 , t) when λ ∈ (0, 1) and S 1 > 0. Then either of the following holds: (i) every solution of (2.3), subject to (2.5), tends to zero as t → ∞, or (ii) (2.3) has a T -periodic solution which is strictly positive (componentwise) at all times, and this solution attracts all solutions with initial data which satisfy (2.5) and are not in the subspace M ; (iii) conclusion (i) (resp., (ii)) holds if the spectral radius of DF (0) is strictly less (resp., larger) than 1, where F is the aforementioned operator which maps the initial datum to the state at time T . Proof. The proof will be via the use of Theorem 2.3.4 on page 48 of Zhao [11]. Recall that C is the positive cone of X and has a nonempty interior. Clearly F : C → C; this follows from the nonnegativity of solutions (Proposition 2.1). Note also that F (0) = 0 because the solution of (2.3) corresponding to zero initial data is the zero solution. Next we shall show that F is strictly subhomogeneous, i.e., that F (λS) > λF (S) for any S ∈ C with S 0 and λ ∈ (0, 1). Now, S 0 means that S ∈ int C and therefore that Verifying F (λS) > λF (S) necessitates verifying that, for initial data satisfying (3.1), the solution S(t; λ) of (3.2) and We will multiply (3.4) by λ and subtract the result from (3.2), and then we will introduce w(t) = S(t; λ) − λS(t). This gives w i (s) ≡ 0 for s ∈ I i and, for t > 0 and i = 1, For t > 0 and i ≥ 2, we find that w i (t) satisfies an equation having exactly the same structure as the ith equation of (2.3), with w i in place of S i , since the ith equation of (2.3) is linear when i ≥ 2.
Since the delay terms all have nonnegative coefficients, a comparison theorem is applicable to the system of differential equations (inequality in the case of w 1 ) for w 1 , . . . , w n and it follows from Theorem 5.2.1 on page 81 of Smith [7] that w i (t) ≥ 0 for t ≥ 0 and all i and therefore in particular that the first statement of (3.3) holds.
We will verify the second statement of (3.3) by showing that, in fact, it holds for i = 1. Suppose, for a contradiction, that S 1 (T + θ; λ) ≡ λS 1 for all t ∈ [T − τ 1 , T ]. But then b(λS 1 (t), t) ≡ λb(S 1 (t), t) on this interval, which in turn implies that S 1 (t) ≡ 0 on [T − τ 1 , T ]. But this is impossible because in the verification of strict subhomogeneity we need only consider initial data S 0. For such initial data, S 1 (0) > 0, which implies that S 1 (t) > 0 for all t > 0, giving a contradiction. We have, therefore, shown that F is strictly subhomogeneous.
Next we shall prove that F N is strongly monotone for some suitably large integer N . To do so involves showing that F N is monotone and, additionally, that

Let S(t) and T (t) be the solutions of (2.3) satisfying the initial data
For t > 0, the differential equations satisfied by the components Y i of Y are of course obtained from (2.3). For i = 2, 3, . . . , n, Y i satisfies the ith equation of (2.3), with Y i in place of S i . For i = 1, we have from an argument similar to that of (3.5) that As ) for the minimal integer N such that NT ≥ nT + n i=1 τ i + max 1≤i≤n τ i . Thus, G := F N is the operator that we work with at this stage. However, a fixed point of F N assures us only of the existence of a periodic solution of (2.3) of period NT rather than T , but this issue will be dealt with later. Let us confirm that F 2 is strictly subhomogeneous (that the same is true for F N can be shown similarly). It is necessary to show that F 2 (λS) > λF 2 (S) when S 0 and λ ∈ (0, 1). Now, F (λS) > λF (S). The above arguments on monotonicity confirm that F itself is monotone, though not necessarily strongly so. Therefore, F 2 (λS) ≥ F (λF (S)). It is easily shown that F (S) 0 when S 0. Therefore, by strict subhomogeneity of F , F (λF (S)) > λF (F (S)). Hence F 2 (λS) > λF 2 (S).
It is also necessary to verify a hypothesis on the Fréchet derivative DF (0) (or DF N (0) if necessary), namely, that it must be compact and strongly positive. This derivative is the linear operator which maps S 0 to S T (or S NT in the case of F N ) as determined by the linearization of (2.3) at its zero solution. To show strong positivity, it is necessary to show that S T (or S NT ) is in int C whenever S 0 > 0. This has effectively been done already because it amounts to the study of a linear system like the system of equations satisfied by the Y i , the i = 1 equation resembling (3.6) but with ∂b(0, t)/∂S 1 in place of 1 0 ∂b(S 1 (t) + ξY 1 (t), t)/∂S 1 dξ. The issues are the same, and strong positivity can be shown for DF N (0), with the same N as that required for strong monotonicity.
We may now apply Theorem 2.3.4 in [11] to the operator F N at this stage. Letting ρ(DF N (0)) denote the spectral radius of DF N (0), we conclude that if ρ(DF N (0)) ≤ 1, then every positive orbit of F N in C converges to 0, and if ρ(DF N (0)) > 1, then there exists a unique fixed point S * 0 of F N in C such that every positive orbit of F N in C \ {0} converges to S * . A fixed point of F N corresponds to a periodic solution of (2.3) of period NT . The statements of the theorem follow if we can show that the solution we have found in the case ρ(DF N (0)) > 1 has period T as well as NT . Note that ρ(DF N (0)) = ρ((DF (0)) N ) = (ρ(DF (0))) N , so ρ(DF (0)) > 1. Moreover, DF (0) is a positive operator, though not necessarily a strongly positive one. Letting ρ * = ρ(DF (0)), we find from Theorem 3.1 that there exists a function S ev > 0 such that DF (0)S ev = ρ * S ev . Next we claim that, in fact, S ev 0. This follows from the fact that S ev will also satisfy DF N (0)S ev = ρ N * S ev and therefore is a positive eigenvector of DF N (0) corresponding to ρ N * . But DF N (0) is a strongly positive operator, so the second part of Theorem 3.1 now applies, with its statement about uniqueness, and informs us that S ev 0. From Taylor's theorem in Banach spaces, In our case, and as applied to δS ev , where δ > 0 is a small number, this gives Therefore, Since ρ * > 1, the right-hand side of the above expression is in int C when δ = 0, and therefore remains in int C when δ > 0 is sufficiently small, by continuity. Therefore, for sufficiently small δ > 0, F (δS ev ) δS ev . We now define a sequence S (k) ∈ X by S (0) = δS ev and S (k+1) = F (S (k) ), k = 0, 1, 2, . . . . We have already established that S (1) > S (0) , and we claim that S (k+1) ≥ S (k) for each integer k ≥ 1. If this is true for a particular k, then it is also true for the next k, since S (k+2) = F (S (k+1) ) ≥ F (S (k) ) = S (k+1) (using that F is monotone), so our claim follows by induction. The monotonicity and boundedness of {S (k) } ∞ k=1 in X implies the monotonicity and boundedness of {S for each fixed θ ∈ I i and each fixed i = 1, 2, . . . , n, so that there exists a function S * * ∈ X with S (k) i (θ) → S * * i (θ) as k → ∞, with θ ∈ I i and i = 1, 2, . . . , n. Note that S * * is a fixed point of F and therefore has period T . Moreover, the T -periodic solution S * * is not the trivial solution, since the above monotone iteration starts with δS ev 0. It remains to show that S * * = S * , but this is obviously so because any subsequence of {S (k) } must also converge to S * * , and so in particular this is true of the subsequence consisting of every N th term. This subsequence is precisely the iterates of F N , which converge to S * . Thus S * * = S * .

Threshold parameter estimation.
It remains to calculate the spectral radius of DF (0). Calculating this radius in a closed form involving the model parameters seems to be a challenging task, so we turn to a qualitatively equivalent issue when this radius is less than or larger than 1. In this section, we consider the case where n = 3, so the linearization of model (2.3) at the zero solution takes the form (4.1) Next we will show that for each fixed t the first three terms in the above expression all tend to zero as k → ∞. We concentrate on the second term; the others can be treated similarly. The second term, if constants are omitted, is which, using the definitions of X λ i (t, s) and d λ i,j (t) and the periodicity (4.6), equals which, after obvious simplification, can be bounded by for each integer m by hypothesis and using the periodicity of the integrands. Also note that z 2 ∈ [−kT, t]. Therefore, Also, if we let M (depending on k and z 2 ) denote the integer such that This inequality can be rewritten in the form The implication of this is that, for the bird species to survive, the birth rate in the breeding patch x 1 plus the rate of return to patch x 1 of birds which left it the previous summer and survived a full year cycle of migration must be larger than the total of the death rate μ 1 at x 1 and the migration rate d 1,2 from patch x 1 to the stopover patch x 2 . In (4.13), the ratio with a numerator of d 2,1 is the rate of departure from x 2 back to x 1 of those birds which made it to x 2 on their outward migration and have since been at x 3 during the winter. The terms in the denominator of that fraction account for the natural death at x 2 , the migration rate from x 2 to x 1 , and what we could call the effective migration rate from x 2 to x 3 . A portion of the birds which left the stopover patch x 2 in the fall migration (with rate d 2,3 ) will not make it back to x 2 in the spring, because they will die at x 3 or in transit. The quantity 1 − α 2,3 α 3,2 d 3,2 /(μ 3 + d 3,2 ) is the fraction that do make the journey from x 2 back to x 2 .

Simulations and discussions.
Migratory birds typically breed in the summer in the northern part of their migratory route and initiate their southward migration route in the fall until reaching their warmer wintering grounds, where they feed during the winter. They initiate a northward migration, returning to their breeding location in the spring. The spring and fall migration routes, which in reality are not necessarily the same, each involve a number of stopovers. The case when n = 3 captures the essential feature of this migration pattern so that, in addition to the summer breeding site x 1 and the wintering location x 3 , we lump together all stopovers as a single site denoted by x 2 . If n = 3 and all periodic coefficients are constants, inequality (4.12) provides the condition for survival of the birds and encapsulates the role of all parameters, including, via α 1,2 , α 2,1 , α 2,3 , and α 3,2 , the in-flight mortalities, distances, and mean flight velocities. Note that if μ 1 > b, then (4.12) is automatically violated, since its left-hand side is less than d 2,1 + d 2,3 (recall that the α parameters are all bounded by 1), so extinction is the outcome in this case. This is not surprising, since births occur only in the breeding patch at x 1 . Survival is therefore possible only if μ 1 < b and (4.12) holds. Even in this case, note that (4.12) is violated if d 2,3 is sufficiently large, since the coefficient of d 2,3 in the left-hand side is less than 1. This requires more interpretation and discussion but can be understood as follows. The right-hand side of (4.12) is the per capita rate of departure from the stopover patch at x 2 , whether by migration or death. The second term in the left-hand side is the per capita rate of arrival back at x 2 of birds returning from the wintering location x 3 . These birds were at x 2 previously, but some will have died at x 3 or during the flight. The term is therefore d 2,3 multiplied by the probabilities of surviving the flight from x 2 to x 3 , of not dying at x 3 , and of surviving the flight from x 3 back to x 2 . These probabilities are α 2,3 , d 3,2 /(μ 3 + d 3,2 ), and α 3,2 , respectively. The prediction of extinction for sufficiently large d 2,3 , all other things being equal, would seem to suggest that it is better for a bird at the stopover location x 2 on its outward migration to simply retrace its steps and go directly back to x 1 , rather than continue to x 3 , at which location there is a possibility of death but not of producing offspring. Here we remind the reader that (4.12) relates only to the autonomous (nonperiodic) version of the model when n = 3. It is only in the autonomous case, and not for realistically chosen periodic parameters, that a bird could fly from x 1 to x 2 and then straight back to x 1 . In nature, each bird would attempt to complete the full migration circuit. These observations highlight the importance of the fact that migration is a periodic phenomenon and that therefore the migration rates should not be constants but realistically chosen periodic functions. Each should be zero at certain times of the year, and such that a bird at x 2 can fly on to x 3 only in the outward migration season, and back to x 1 only in the return season, and not given the choice, as the autonomous model implies. Condition (4.12) is, however, sharp if migration is temporally homogeneous.
We carried out some numerical simulations to illustrate the threshold dynamics by using some data from the U.S. Geological Survey [8], which tracked three of a dozen bar-headed geese using satellites. Bar-headed geese can fly up to 50,000 km from northern Mongolia to India during their migration cycle. The flocks follow various routes, but their major stopover locations seem to be common, and both the spring and fall migrations involve between five and six stopovers for the tracked bar-headed geese traveling from Mongolia to India. In our simulations, we use a single patch to denote the collection of these stopovers. The birth function was chosen as b(S 1 , t) = rS 1 (1 − S 1 /K). Solutions always remain below K. The other parameter values were taken from Bourouiba and Wu [2], Javed et al. [5], and Prins and van Wieren [6]. We divide each year into the periods (0, T 1 ), (T 1 , T 2 ), (T 2 , T 3 ), and (T 3 , 365), with time measured in days, and define May 1 to be the beginning of each year. The breeding season is from time t = 0 to t = T 1 = 138. The outward migration (fall migration) is from time T 1 to T 2 , where T 2 − T 1 = 61 days. From time T 2 to T 3 , birds stay at a winter feeding patch. From time t = T 3 to t = 365, inward migration (spring migration) is in progress. The parameters associated with the birth function were taken as r = 0.005 and K = 60,000. The mortality rates in the first (summer breeding) patch, the second (stopover) patch, and the third (winter feeding) patch were taken as μ 1 = 0.00088, μ 2 = 0.00088, and μ 3 = 0.001. We chose the migration rates for outward migration to be d i,i+1 (t) = 0.3, i = 1, 2, when t mod 365 ∈ (T 1 , T 2 ) and zero at other times. For return migration, d i+1,i (t) = 0.3, i = 1, 2, when t mod 365 ∈ (T 3 , 365) and is zero at other times. The values of the in-flight mortalities μ i−1,i and μ i+1,i were 0.002. For the flight times τ i = l i /U i , we chose τ i = 5 for each i. These values enter the model through the α parameters α 12 = α 23 = α 32 = α 21 = 0.99. Figure 1 shows a periodic solution of a period of one year, reflecting the one-year periodicity of the migration rates. In the plots, only birds occupying patches are counted, but even the total bird population, including those in the air, may oscillate, as there is no single equation governing the total number of birds. In another simulation, the result of which is shown in Figure 2, the death rate on the summer breeding patch has been increased tenfold with other parameters remaining the same. The result is that the bird population goes to extinction. Other simulations suggest that it is necessary only to increase this death parameter between fourfold and fivefold to realize extinction rather than a positive periodic solution, but extinction occurs very slowly unless the parameter is increased considerably more.
In this paper, we have concentrated solely on the issue of the mathematical modeling of Increasing the death rate tenfold on the first patch (summer breeding) and keeping the other parameters the same as in Figure 1. The bird population dies out regardless of the initial conditions. the bird migration phenomenon and in particular on how to correctly deal with the flying times between stopovers and the issue of in-flight mortality. Our longer term goal is to understand the possible role of migratory birds in the spread of H5N1 between continents. In further papers presently in preparation, we are further developing our modeling methodology to include the interaction of migratory birds with nonmigratory species and especially with poultry. Some migratory bird species follow different routes for outward and return journeys and therefore incur different local conditions, such as density of poultry farms and different local control strategies for H5N1 control, and we think it is especially important to compare the impact of the outward and return journeys on disease transmission.