Seasonal migration dynamics: periodicity, transition delay and finite-dimensional reduction

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.


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 *Author for correspondence (wujh@mathstat.yorku.ca). summer and winter sites, respectively, to obtain the following system of delay differential equations (DDEs): x s (t) = −(m sw (t) + m s )x s (t) + e −m ws t ws m ws (t − t ws )x w (t − t ws ) andẋ w (t) = −(m ws (t) + m w )x w (t) + e −m sw t sw m sw (t − t sw )x s (t − t sw ). (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 t ws represents the time (days) flying from the winter site to the summer site, while t sw is the time (days) flying from the summer site to the winter site. The corresponding rates of in-flight mortality are denoted by m ws and m sw , and the death rates in the summer and winter sites, respectively, are denoted by m s and m w . The birth rate in the summer site is the intrinsic birth rate g 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 k ∈ N. 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 t 0 t 1 t 2 t 3 s p r i n g m i g r a t i o n b e g i n s s p r i n g m i g r a t i o n e n d s a u t u m n m i g r a t i o n b e g i n s a u t u m n m i g r a t i o n e n d s site is t 1 + t 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 + t ws < t 2 and t 3 + t sw < t 0 + T . In what follows, we assume and g s (t) = g > 0, t ∈ (t 0 + t ws , t 2 ), 0, t ∈ (t 0 , t 0 + t ws ) ∪ (t 2 , t 0 + T ). (1.4) In the final section, we will discuss how these assumptions can be relaxed.

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) S 2 (t), t ∈ (t 0 + t ws , t 1 + t ws ), S 3 (t), t ∈ (t 1 + t ws , t 2 ), S 4 (t), t ∈ (t 2 , t 3 ), on July 20, 2018 http://rspa.royalsocietypublishing.org/ Downloaded from Various curves given by x s (t) and x w (t) in different subintervals, these curves are obtained by solving the system of ODEs (2.3a)-(2.3i). and It is easily seen that (S 1 (t 0 ), W 1 (t 0 )) = (x s (t 0 ), x w (t 0 )) anḋ 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, a, f 0 , t), the unique solution of the equatioṅ 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.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: . Illustration for the derivation of the system of ordinary differential equations equivalent to the bird migration model which is time-periodic and involves multiple delays.
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 t ws = 16 and t sw = 16, and all units are days. The death rates are m s = 0.00088 and m w = 0.00088. The values of in-flight mortalities m sw and m ws are also set to 0.00088. The birth rate in the summer breeding site is g = 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,

Reduction to a discrete system and the threshold value
Assuming that 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: x s (kT + T 1 + t ws )

∼
(1 − m s /g)K 1 − e −(g−m s )T 1 · g(2/(M ws + m w ) gx w (kT )M ws /K exp(m ws t ws )) , where g(z) := G(−n + 1)I −n−1 (z)(z/2) 2n G(n + 1)I n+1 (z) with n := (g − m s )/(M ws + m w ). It is easily seen from equation (2.5) that as z → 0, It is then natural to introduce a discrete system {X n } ∞ n=0 with X 0 := x w (0) > 0, and for each k ∈ N, . (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 G : R → R 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 For the aforementioned parameter values, if we decrease the value of g 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 g, we can compute the values of G (0) and x s (kT + T 1 + T 2 ) for sufficiently large k ∈ N; see table 1. This numerical evidence shows that as G (0) decreases to 1, the bird population declines substantially.

Threshold value of the original system
We can now rewrite G (0) as G (0) = P s · P sw · P w · P ws · t ws · t sw with P s := exp[(g − m s )(T 1 + T 2 − t ws )] being the reproduction ratio in the summer; on July 20, 2018 http://rspa.royalsocietypublishing.org/ Downloaded from Table 1. The bird population declines significantly as G (0) decreases to 1. This, in agreement of the threshold phenomena at G (0) = 1 for the discrete map G, indicates the close relationship between G (0) and the threshold dynamics of model system (1.1). being the survival probability during the autumn migration; being the survival probability during the winter; being the survival probability during the spring migration; t ws := M ws M ws + m w + g − m s being the proportion of birds in transition from the winter site to the summer site, and t sw := M sw M sw + m s − m w 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 F : R 2 → R 2 be the mapping defined as F (x s (0), x w (0)) := (x s (T ), x w (T )).
By applying the above six formulas to equation (4.10), we obtain as c → 0, Therefore, Here, we have again used the fact n = a/a. This ends the proof.

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, u i (t) = h i (t, u i1 t , . . . , u in t ), (5.1) where u ij t ∈ C [−r ij , 0] is defined by u ij t (q) = u i (t + q) for q ∈ [−r ij , 0], and h i is a functional on R × n j=1 C [−r ij , 0] (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 t j = max 1≤i≤n r ij , and let r i = max 1≤j≤n r ij 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 h i (i = 1, . . . , n) such that for any t ∈ [t 0 , t 0 + r i ] and any (dummy variables) (4 i1 , . . . , 4 in ) ∈ n j=1 C [−r ij , 0], h i (t, 4 i1 , . . . , 4 in ) =h i (t, 4 ii (0)). (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 [t 0 + r i , ∞) 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 that the finite-dimensional reduction holds for the general model of Bourouiba