## Abstract

We study a model of bird migration between the summer breeding ground and the winter refuge site. The model involves time lags for the transition time between the patches, and the model parameters are periodic in time as the biological activities related to the migration and reproduction are seasonal. It has been shown in previous studies that the model system exhibits threshold dynamics: either all solutions converge to the trivial solution, or the system has a positive and globally attractive periodic solution. Two issues remain and will be addressed in this paper: how to express the threshold condition in terms of model parameters explicitly (rather than the abstract spectral radius of a certain monodromy operator) and how to describe the temporal pattern of the positive periodic solution. We make an interesting and surprising observation that the delay differential system is completely characterized by a finite-dimensional ordinary differential system and then a finite-dimensional map in the sense that the bird population at the initial time of spring migration determines the future status of the system. As a consequence, we derive the threshold condition, explicitly in terms of the model parameters, for the extinction and persistence of the considered bird species.

## 1. The dynamic model of migrant birds

We follow Bourouiba *et al.* (2010) and Gourley *et al.* (2010) to consider a single species bird population migrating between a summer breeding patch and a winter refuge patch. We will ignore stopover patches, though mortality during the flights to these patches and at these stopover patches will be incorporated in the model parameters indirectly. We denote by *x*_{s}(*t*) and *x*_{w}(*t*) the numbers of birds in the summer and winter sites, respectively, to obtain the following system of delay differential equations (DDEs):
1.1a
and
1.1b
Here, the functions *m*_{ws}(*t*) and *m*_{sw}(*t*) describe the migration rates between the two patches (w, winter; s, summer; ws, from the winter patch to the summer patch; sw, from the summer patch to the winter patch). The time delay *τ*_{ws} represents the time (days) flying from the winter site to the summer site, while *τ*_{sw} is the time (days) flying from the summer site to the winter site. The corresponding rates of in-flight mortality are denoted by *μ*_{ws} and *μ*_{sw}, and the death rates in the summer and winter sites, respectively, are denoted by *μ*_{s} and *μ*_{w}. The birth rate in the summer site is the intrinsic birth rate *γ*_{s}(*t*) regulated by the density through a logistic term. Naturally, all the coefficients in the system are periodic functions with the period *T*=365 days.

The above model represents a simplified version of the patchy model developed in Bourouiba *et al.* (2010) and Gourley *et al.* (2010) to capture the essential dynamic features of bird migration, and to explore the roles of bird migration in the global spread of highly pathogenic H5N1 avian influenza. We refer to these two studies for some of the basic qualitative properties of the above model system: non-negativeness, boundedness and dissipativeness, to name a few. We also note that despite the long history and intensive study of bird migration, mathematical modelling and analysis of bird migration have been difficult until recently when reliable surveillance became feasible thanks to the advancement in satellite tracking of migrant birds and with global threshold dynamics of periodic systems being possible due to the development of monotone dynamical systems and persistence theory. The aforementioned studies parametrized the model using satellite tracking records of a dozen Bar-headed geese migrating between the summer breeding site in the northern part of their migration path (e.g. Mongolia) and the wintering grounds (e.g. India), and also established the model’s threshold dynamics under very general conditions on the birth rate. In particular, they proved that either all solutions converge to the trivial solution, or the system has a positive and globally attractive periodic solution. However, for the case of applications, it is very important to develop effective techniques to express the threshold condition in terms of model parameters explicitly (rather than the abstract spectral radius of a certain monodromy operator), and to describe the temporal pattern of the positive periodic solution.

In this paper, we will develop these techniques using some of the unique features of bird migration seasonality. In what follows, let *t*_{0}=*kT* for an arbitrarily fixed . We normalize the time so *t*_{0} is the starting date when the birds begin to fly to the summer breeding site in a particular year. We let *t*_{1} be the time when the birds in the winter patch stop their spring migration to the summer breeding site. Therefore, the time when the last spring migratory bird arrives at the summer site is *t*_{1}+*τ*_{ws}. We then assume that, after summer breeding, the birds start their autumn migration at time *t*_{2} and autumn migration ends at date *t*_{3} (see figure 1).

Let *T*_{1}:=*t*_{1}−*t*_{0}; *T*_{2}:=*t*_{2}−*t*_{1}; *T*_{3}:=*t*_{3}−*t*_{2} and *T*_{4}:=*t*_{0}+*T*−*t*_{3} represent the durations of the aforementioned biological activities; we have *T*_{1}+*T*_{2}+*T*_{3}+*T*_{4}=*T*, and assume *t*_{1}+*τ*_{ws}<*t*_{2} and *t*_{3}+*τ*_{sw}<*t*_{0}+*T*. In what follows, we assume
1.2
1.3
and
1.4
In the final section, we will discuss how these assumptions can be relaxed.

## 2. Explicit solutions and numerical simulations

A major observation is that a solution of system (1.1) is given by an enlarged system of periodic ordinary differential equations (ODEs), with state variables corresponding to the numbers of birds in different intervals of a given year. Namely, we define (see figure 2) 2.1 and 2.2

It is easily seen that (*S*_{1}(*t*_{0}),*W*_{1}(*t*_{0}))=(*x*_{s}(*t*_{0}),*x*_{w}(*t*_{0})) and
2.3a
2.3b
2.3c
2.3d
2.3e
2.3f
2.3g
2.3h
2.3i

It is important to note that the intervals defining *S*_{i}(*t*), *i*=1,…,5, and the intervals defining *W*_{j}(*t*), *j*=1,…,4, give partitions of [*t*_{0},*t*_{0}+*T*], and these partitions are different for *S*_{i}(*t*) and *W*_{j}(*t*). The choice of these intervals is made so that solutions of equation (1.1) with the initial data (*x*_{s}(*t*_{0}),*x*_{w}(*t*_{0})) can be obtained by solving the ODEs (2.3a)–(2.3i) consequently. To be more specific, we introduce a function *f*(*a*,*b*,*c*,*α*,*f*_{0},*t*), the unique solution of the equation
2.4
subject to the initial condition *f*(0)=*f*_{0}. This equation can be solved in terms of modified Bessel functions defined as (Olver 1997)
2.5
and
2.6
Set , and
2.7
We have
2.8
The nine equations in equation (2.3) can be solved step by step, and the solutions are given as follows:
2.9a
2.9b
2.9c
2.9d
2.9e
2.9f
2.9g
2.9h
2.9i

To obtain the above expressions, we used facts such as for equation (2.9a), *m*_{sw}(*t*)=0 for *t*∈[*t*_{0}+*T*−*τ*_{ws},*t*_{0}+*T*] since *t*_{3}+*τ*_{ws}<*t*_{0}+*T*; and for equation (2.9b), *γ*_{s}(*t*)=0 for *t*∈(*t*_{0},*t*_{0}+*τ*_{ws}) and *m*_{ws}(*t*)=0 for *t*∈(*t*_{0}−*τ*_{ws},*t*_{0}) (see figure 3 for an illustration).

A significant and amazing observation is that the solutions *S*_{i}(*t*) and *W*_{j}(*t*) for *i*=1,…,5 and *j*=1,…,4 are uniquely determined as long as (*S*_{1}(*t*_{0}),*W*_{1}(*t*_{0}))=(*x*_{s}(*t*_{0}),*x*_{w}(*t*_{0})) is specified. Therefore, solutions of delay differential system (1.1) with periodic coefficients are uniquely determined by the initial values of (*x*_{s}(*t*_{0}),*x*_{w}(*t*_{0})) at a *single time* rather than on an initial interval.

To illustrate this, we provide a numerical simulation based on parameters suggested in Bourouiba *et al.* (2010): *T*_{1}=46, *T*_{2}=138, *T*_{3}=61 and *T*_{4}=120, the time delays are taken as *τ*_{ws}=16 and *τ*_{sw}=16, and all units are days. The death rates are *μ*_{s}=0.00088 and *μ*_{w}=0.00088. The values of in-flight mortalities *μ*_{sw} and *μ*_{ws} are also set to 0.00088. The birth rate in the summer breeding site is *γ*=0.00499, and the carrying capacity is *K*=60 000. We assume that 99.99 per cent of the birds depart during the migratory seasons. That is,
and
Consequently, the migratory rates are and . The simulation shows that the bird populations (*x*_{s}(*t*),*x*_{w}(*t*)) eventually converge to periodic solutions (figure 4).

## 3. Reduction to a discrete system and the threshold value

Assuming that
3.1
is sufficiently small, we obtain (by setting *t*_{0}=*iT* for integer *i* varying iteratively), from equations (2.1), (2.2), (2.8) and (2.9), the following approximations for the value of the bird populations in the summer and winter patches at the end of a particular biological activity interval:
where
with *ν*:=(*γ*−*μ*_{s})/(*M*_{ws}+*μ*_{w}). It is easily seen from equation (2.5) that as *z*→0,
3.2
It is then natural to introduce a discrete system with *X*_{0}:=*x*_{w}(0)>0, and for each ,
3.3a
3.3b
3.3c
and
3.3d
Note that in the above formulation, *T*_{1}, *T*_{2}, *T*_{3} and *T*_{4} are all defined when *t*_{0} is set to zero. Let be the mapping
This is a monotone map. The survival/extinction threshold of this discrete system is characterized by *G*′(0), the derivative of *G* at zero. It follows from equations (3.2) and (3.3) that
3.4
For the aforementioned parameter values, if we decrease the value of *γ* from 0.00499 to 0.00199 with the values of other parameters fixed, we have the value of *G*′(0) decreased from 1.63641 to 1.00324. For each *γ*, we can compute the values of *G*′(0) and *x*_{s}(*kT*+*T*_{1}+*T*_{2}) for sufficiently large ; see table 1. This numerical evidence shows that as *G*′(0) decreases to 1, the bird population declines substantially.

## 4. Threshold value of the original system

We can now rewrite *G*′(0) as
with
being the reproduction ratio in the summer;
being the survival probability during the autumn migration;
being the survival probability during the winter;
being the survival probability during the spring migration;
being the proportion of birds in transition from the winter site to the summer site, and
being the proportion of birds in transition from the summer site to the winter site. Consequently, we have a clear ecological interpretation of *G*′(0) as the annual reproduction rate. We can now mathematically show that this reproduction rate determines the global dynamics. Namely, we let be the mapping defined as
4.1
This is a monotone map, and hence the extinction threshold of the system in equation (1.1) is characterized by *ρ*(*DF*(0)), the spectral radius of *DF*(0) (Zhao 2003; Gourley *et al.* 2010).

### Theorem 4.1

*As ε→0, we have
*
4.2

Before proving this theorem, we first state two formulas about the partial derivatives of the function *f* defined in equation (2.8).

### Lemma 4.2

*Let* *f* *be defined as in equation* . *Let* *a*, *b*, *α* *and* *t* *be fixed and regard* *c* *and* *f*_{0} *as two variables. We have*
4.3
*and*
4.4

### Proof.

From the definitions of modified Bessel functions (2.5) and (2.6), we have as *z*→0,
4.5
and
4.6
Applying the above two formulas to equation (2.7) yields
as *c*→0, where
4.7a
and
4.7b
Recall that , we obtain from equations (4.5) and (4.6) that as *c*→0,
where
Set . We have from equation (2.8) that as *c*→0,
Hence,
As *f*_{0}→0, it follows from equation (4.7) that *λ*_{0}→0 and
Thus, we obtain
where we have used the fact *ν*=*a*/*α*. Now we calculate the partial derivative of *f* with respect to *f*_{0}. From equation (2.7), we have as *f*_{0}→0,
where
4.8
and
4.9
Hence, from equation (2.8), we have
Consequently,
4.10
Recall . It can be shown from equations (4.5) and (4.6) that as *c*→0,
By applying the above six formulas to equation (4.10), we obtain as *c*→0,
Therefore,
Here, we have again used the fact *ν*=*a*/*α*. This ends the proof.■

Now we are ready to prove equation (4.2) in our main theorem. From the definitions (2.1), (2.2) and (4.1), we have
where
Note that *t*_{0}=0 and *DF*(0) means the value of *DF* taken at the point (*S*_{1}(*t*_{0}),*W*_{1}(*t*_{0}))=(0,0). Similar notations apply to *M*_{ij}(0) with *i*,*j*=1 or 2. It is easily seen from equations (2.9), (4.3) and (4.4) that
As , we have *M*_{11}(0)=*O*(*ε*), *M*_{12}(0)=*O*(*ε*), *M*_{21}(0)=*O*(1) and
Therefore, the spectrum radius of *DF*(0) has the asymptotic formula
as *ε*→0. This proves equation (4.2).

## 5. Conclusion and remarks

Determining a solution of a delay differential system normally requires specifying initial values of the solution in a given interval. Here we show, however, that for the patchy model of bird migrations, even in the presence of positive delays due to the transition times between patches (the summer breeding ground and the winter refuge site), we can solve the delay differential system (forward in time) with the only knowledge of the initial value of the bird population at a single point of time (the time when the spring migration starts). We made this observation by linking the delay differential system with periodic coefficients to an enlarged system of ODEs defined in two disjoint unions of the time interval of length *T*, the period. Solutions of such an enlarged ODE system can be obtained consecutively due to the nature of biological activities of bird migration and reproduction in an interwinding spatio-temporal fashion. Our study shows that specific techniques can and should be developed to investigate the population dynamics in a seasonally varying environment.

We end this paper with a few remarks of how our results and techniques can be extended to more general situations.

### (a) Finite dimensional reduction

We note the possibility of finite-dimensional reduction for periodic systems of DDEs, if the delays are smaller than the period and the periodicity arises due to seasonal biological activities. We illustrated this with the two-patch bird migration system. However, this reduction is possible for a very general class of periodic systems of DDEs. We now provide a sufficient condition for reducing a periodic delay differential system into an ODE system. Let *u*(*t*)=(*u*^{1}(*t*),…,*u*^{n}(*t*)) be a vector-valued function such that for each *i*=1,2,…,*n*,
5.1
where is defined by for *θ*∈[−*r*_{ij},0], and *h*^{i} is a functional on (in general, *u*^{i} with *i*=1,2,…,*n* can also be a vector-valued function, while *h*^{i} can be a corresponding vector-valued functional). Define , and let be the maximum of delays in the expression of *h*^{i}. We want to show that the above delay differential system (5.1) is reducible to an ODE system if there exist *t*_{0}≥0 and functions (*i*=1,…,*n*) such that for any *t*∈[*t*_{0},*t*_{0}+*r*_{i}] and any (dummy variables) ,
5.2
The above assumption requires that the DDE for *u*^{i} is reducible to an ODE on the interval [*t*_{0},*t*_{0}+*r*_{i}]. Thus, the values of *u*^{i} on this interval are determined by the initial value *u*(*t*_{0})=(*u*^{1}(*t*_{0}),…,*u*^{n}(*t*_{0})). We are left to prove that the values of *u*^{i} on the interval are determined by the initial value *u*(*t*_{0}). Without loss of generality, we may assume that *r*_{1}≤*r*_{2}≤⋯≤*r*_{n}. Given any initial value at *t*_{0}, we have solved *u*^{i} on the interval [*t*_{0},*t*_{0}+*r*_{i}]. On the interval [*t*_{0}+*r*_{1},*t*_{0}+*r*_{2}], we consider the DDE for *u*^{1}:
For any *t*∈[*t*_{0}+*r*_{1},*t*_{0}+*r*_{2}] and *θ*∈[−*r*_{1},0], we have *t*+*θ*∈[*t*_{0},*t*_{0}+*r*_{2}]. Thus, the values of with *j*≥2 is known for any *t*∈[*t*_{0}+*r*_{1},*t*_{0}+*r*_{2}] and *θ*∈[−*r*_{1},0]. Furthermore, the initial values for *θ*∈[−*r*_{1},0] are also known. Thus, the above DDE can be solved on the interval [*t*_{0}+*r*_{1},*t*_{0}+*r*_{2}]. Then we consider the DDE system for (*u*^{1},*u*^{2}) on the interval [*t*_{0}+*r*_{2},*t*_{0}+*r*_{3}]:
and
Since the values of *u*^{j}(*t*+*θ*) with *j*≥3 are known for any *t*∈[*t*_{0}+*r*_{2},*t*_{0}+*r*_{3}] and *θ*∈[−*r*_{2},0], and since the initial values and for *θ*∈[−*r*_{2},0] are also known, the above DDE system is solvable on the interval [*t*_{0}+*r*_{2},*t*_{0}+*r*_{3}]. We can proceed in a similar manner to solve the whole DDE system for *u*(*t*)=(*u*^{1}(*t*),…,*u*^{n}(*t*)) on the interval [*t*_{0},*t*_{0}+*r*_{n}] when the initial value at the single point *t*_{0} is given. For *t*≥*t*_{0}+*r*_{n}, the DDE system is also solvable because the initial values for any *θ*∈[−*r*_{n},0] with *j*=1,2,…,*n* are already obtained. This verifies that the DDE system can be reduced to an ODE system (5.1) if condition (5.2) is satisfied.

We now demonstrate that our system (1.1) satisfies the condition mentioned above if *t*_{1}+*τ*_{ws}<*t*_{2} and *t*_{3}+*τ*_{sw}<*t*_{0}+*T*, and *m*_{ws}, *m*_{sw} are piecewise constants as given in equations (1.2) and (1.3), respectively. In this case, *h*=(*h*^{1},*h*^{2}), where
and
Here, for any ,
and
It is easy to verify that *r*_{11}=0, *r*_{1}=*r*_{12}=*τ*_{ws}, *r*_{2}=*r*_{21}=*τ*_{sw} and *r*_{22}=0. For *t*∈[*t*_{0},*t*_{0}+*r*_{1}], we have *m*_{ws}(*t*−*τ*_{ws})=0 on account of equation (1.2) and *t*_{1}+*τ*_{ws}<*t*_{2}<*t*_{0}+*T*; noting that *m*_{ws}(*t*) is a periodic function with period *T*. Thus,
on the interval [*t*_{0},*t*_{0}+*r*_{1}]. Similarly, by virtue of equation (1.3), inequality *t*_{3}+*τ*_{sw}<*t*_{0}+*T* and periodic property of *m*_{sw}, we can show that
on the interval [*t*_{0},*t*_{0}+*r*_{2}]. Therefore, the condition (5.2) is satisfied and the DDE system (1.1) is reducible to a finite-dimensional system. We note, in particular, that the finite-dimensional reduction holds for the general model of Bourouiba *et al.* (2010) and Gourley *et al.* (2010) with stopovers during migrations.

### (b) Carrying capacity and general birth rate

We also remark that under our assumptions, during the breeding time (*t*_{0}+*τ*_{ws},*t*_{2}), equation (1.1a) involves a logistic birth term *γx*_{s}(*t*)(1−*x*_{s}(*t*)/*K*) and a death term −*μ*_{s}*x*_{s}(*t*), a combination of which gives , where and . In this sense, can be regarded as the carrying capacity of the summer breeding site.

We assume the birth rate *γ*_{s}(*t*) in equation (1.4) to be constant in the interval (*t*_{0}+*τ*_{ws},*t*_{2}), which helps us to express the solution of equation (2.4) in terms of Bessel function. If, in general,
5.3
then the coefficients *a* and *b* in equation (2.4) are no longer constants and we can no longer find an explicit solution to this general equation. However, we are still able to derive an asymptotic formula for the threshold. In this case, the reproduction ratio in the summer (cf. §4)
should be replaced by
where
is the average of *γ*(*t*) in the interval (*t*_{0}+*τ*_{ws},*t*_{2}); recalling that *t*_{2}−*t*_{0}=*T*_{1}+*T*_{2}. Another change should be made on the transition proportion from the winter site to the summer site (cf. §4), where
should be replaced by
To compare this formula with that of constant birth rate, we repeat integration by parts and obtain
Here, we have omitted the terms involving which is small; see equation (3.1).

### (c) General migration rate

We assumed the migration rate *m*_{ws}(*t*) in equation (1.2) (resp. *m*_{sw}(*t*) in equation (1.3)) to be constant in its support interval (*t*_{0},*t*_{1}) (resp. (*t*_{2},*t*_{3})). If, in general,
5.4
and
5.5
then we can also obtain the approximated threshold similar to that given in the beginning of §4. Actually, we only need to replace the transition proportion from the winter site to the summer site
by
and the transition proportion from the summer site to the winter site
by
Note that if *M*_{ws}(*t*)≡*M*_{ws} is a constant, then
Similarly, if *M*_{sw}(*t*)≡*M*_{sw} is a constant, then
It is also easy to extend our result to the case when the assumptions (1.2), (1.3) and (1.4) are replaced by equations (5.4), (5.5) and (5.3), respectively.

## Acknowledgements

This work is partially supported by CRC, GEOIDE, MITACS, NSERC and Mprime. We are grateful to the referees for their insightful comments and constructive criticism, which have been answered in detail in the last section.

- Received April 14, 2011.
- Accepted August 15, 2011.

- This journal is © 2011 The Royal Society

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.