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 xs(t) and xw(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 mws(t) and msw(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 t0=kT for an arbitrarily fixed . We normalize the time so t0 is the starting date when the birds begin to fly to the summer breeding site in a particular year. We let t1 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 t1+τws. We then assume that, after summer breeding, the birds start their autumn migration at time t2 and autumn migration ends at date t3 (see figure 1).
Let T1:=t1−t0; T2:=t2−t1; T3:=t3−t2 and T4:=t0+T−t3 represent the durations of the aforementioned biological activities; we have T1+T2+T3+T4=T, and assume t1+τws<t2 and t3+τsw<t0+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 (S1(t0),W1(t0))=(xs(t0),xw(t0)) 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 Si(t), i=1,…,5, and the intervals defining Wj(t), j=1,…,4, give partitions of [t0,t0+T], and these partitions are different for Si(t) and Wj(t). The choice of these intervals is made so that solutions of equation (1.1) with the initial data (xs(t0),xw(t0)) can be obtained by solving the ODEs (2.3a)–(2.3i) consequently. To be more specific, we introduce a function f(a,b,c,α,f0,t), the unique solution of the equation 2.4 subject to the initial condition f(0)=f0. 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), msw(t)=0 for t∈[t0+T−τws,t0+T] since t3+τws<t0+T; and for equation (2.9b), γs(t)=0 for t∈(t0,t0+τws) and mws(t)=0 for t∈(t0−τws,t0) (see figure 3 for an illustration).
A significant and amazing observation is that the solutions Si(t) and Wj(t) for i=1,…,5 and j=1,…,4 are uniquely determined as long as (S1(t0),W1(t0))=(xs(t0),xw(t0)) is specified. Therefore, solutions of delay differential system (1.1) with periodic coefficients are uniquely determined by the initial values of (xs(t0),xw(t0)) 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): T1=46, T2=138, T3=61 and T4=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 (xs(t),xw(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 t0=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)/(Mws+μ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 X0:=xw(0)>0, and for each , 3.3a 3.3b 3.3c and 3.3d Note that in the above formulation, T1, T2, T3 and T4 are all defined when t0 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 xs(kT+T1+T2) 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).
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).
Let f be defined as in equation (2.8). Let a, b, α and t be fixed and regard c and f0 as two variables. We have 4.3 and 4.4
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 f0→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 f0. From equation (2.7), we have as f0→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 t0=0 and DF(0) means the value of DF taken at the point (S1(t0),W1(t0))=(0,0). Similar notations apply to Mij(0) with i,j=1 or 2. It is easily seen from equations (2.9), (4.3) and (4.4) that As , we have M11(0)=O(ε), M12(0)=O(ε), M21(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)=(u1(t),…,un(t)) be a vector-valued function such that for each i=1,2,…,n, 5.1 where is defined by for θ∈[−rij,0], and hi is a functional on (in general, ui with i=1,2,…,n can also be a vector-valued function, while hi can be a corresponding vector-valued functional). Define , and let be the maximum of delays in the expression of hi. We want to show that the above delay differential system (5.1) is reducible to an ODE system if there exist t0≥0 and functions (i=1,…,n) such that for any t∈[t0,t0+ri] and any (dummy variables) , 5.2 The above assumption requires that the DDE for ui is reducible to an ODE on the interval [t0,t0+ri]. Thus, the values of ui on this interval are determined by the initial value u(t0)=(u1(t0),…,un(t0)). We are left to prove that the values of ui on the interval are determined by the initial value u(t0). Without loss of generality, we may assume that r1≤r2≤⋯≤rn. Given any initial value at t0, we have solved ui on the interval [t0,t0+ri]. On the interval [t0+r1,t0+r2], we consider the DDE for u1: For any t∈[t0+r1,t0+r2] and θ∈[−r1,0], we have t+θ∈[t0,t0+r2]. Thus, the values of with j≥2 is known for any t∈[t0+r1,t0+r2] and θ∈[−r1,0]. Furthermore, the initial values for θ∈[−r1,0] are also known. Thus, the above DDE can be solved on the interval [t0+r1,t0+r2]. Then we consider the DDE system for (u1,u2) on the interval [t0+r2,t0+r3]: and Since the values of uj(t+θ) with j≥3 are known for any t∈[t0+r2,t0+r3] and θ∈[−r2,0], and since the initial values and for θ∈[−r2,0] are also known, the above DDE system is solvable on the interval [t0+r2,t0+r3]. We can proceed in a similar manner to solve the whole DDE system for u(t)=(u1(t),…,un(t)) on the interval [t0,t0+rn] when the initial value at the single point t0 is given. For t≥t0+rn, the DDE system is also solvable because the initial values for any θ∈[−rn,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 t1+τws<t2 and t3+τsw<t0+T, and mws, msw are piecewise constants as given in equations (1.2) and (1.3), respectively. In this case, h=(h1,h2), where and Here, for any , and It is easy to verify that r11=0, r1=r12=τws, r2=r21=τsw and r22=0. For t∈[t0,t0+r1], we have mws(t−τws)=0 on account of equation (1.2) and t1+τws<t2<t0+T; noting that mws(t) is a periodic function with period T. Thus, on the interval [t0,t0+r1]. Similarly, by virtue of equation (1.3), inequality t3+τsw<t0+T and periodic property of msw, we can show that on the interval [t0,t0+r2]. 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 (t0+τws,t2), equation (1.1a) involves a logistic birth term γxs(t)(1−xs(t)/K) and a death term −μsxs(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 (t0+τws,t2), 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 (t0+τws,t2); recalling that t2−t0=T1+T2. 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 mws(t) in equation (1.2) (resp. msw(t) in equation (1.3)) to be constant in its support interval (t0,t1) (resp. (t2,t3)). 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 Mws(t)≡Mws is a constant, then Similarly, if Msw(t)≡Msw 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.
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.