## Abstract

Regular, long-distance migrations of thousands of animal species have consequences for the ecosystems that they visit, modifying trophic interactions and transporting many non-pathogenic and pathogenic organisms. The spatial structure and dynamic properties of animal migrations and population flyways largely determine those trophic and transport effects, but are yet poorly studied. As a basis, we propose a periodic Markov model on the spatial migration network of breeding, stopover and wintering sites to formally describe the process of animal migration on the population level. From seasonally changing transition rates we derived stable, seasonal densities of animals at the network nodes. We parametrized the model with high-quality GPS and satellite telemetry tracks of white storks (*Ciconia ciconia*) and greater white-fronted geese (*Anser a. albifrons*). Topological and network flow properties of the two derived networks conform to migration properties like seasonally changing connectivity and shared, directed movement. Thus, the model realistically describes the migration movement of complete populations and can become an important tool to study the effects of climate and habitat change and pathogen spread on migratory animals. Furthermore, the property of periodically changing transition rates makes it a new type of complex model and we need to understand its dynamic properties.

### 1. Introduction

Animal migration is a large-scale biological phenomenon that has fascinated humankind since ancient times. Thousands of species of birds, reptiles and mammals perform seasonal migrations, amounting to millions of individuals travelling over hundreds and often thousands of kilometres following highly predictable time schedules [1]. These mass displacements of animals have many consequences for the animals themselves, the ecosystems they visit, and also for humans [2]. While migrating, animals modify local food-web structures and predator–prey interactions (trophic effects) and transport a wide variety of non-pathogenic and pathogenic organisms, such as plant seeds and parasites (transport effects). In those respects, ‘stepping stone' migrants are especially important, because they visit and modify conditions at several stopover sites along their route where they spend days or weeks to rest and replenish their fuel reserves [3–5].

Bird migration is probably the most studied type of animal movement and much work has been devoted to understanding the physiological and ecological mechanisms of the migration process, metabolism, navigation and flight [5,6]. However, many questions remain unanswered, especially regarding the spatial structure and dynamic properties of migration and stopover usage. These properties, including route choice, number of stopovers, stopover duration and migration speed are particularly relevant in the context of the trophic and transport effects of migration [2]. They largely determine the contribution to local ecosystems processes [7], population dynamics of migrants [8] and the probability of disease transmission [9,10].

The field of animal tracking with GPS and satellite telemetry is presently growing tremendously and travels of more and more migratory species are recorded [11], improving our ability to study the temporal and spatial properties of migration. However, it is always individual animals that are tracked and usually there is large individual variation in the allocation of time and energy during their travels [12]. Thus, directly studying migration on the population or species level with tracking data is difficult. It is not straightforward how the overall migration process can be formalized with quantitative models that are able to capture its original properties with a non-infinite number of parameters.

Recently, the framework of network theory has been proposed to describe the spatio-temporal process of animal migration [13–15]. Network-based tools are already successfully used in other branches of animal ecology [16], mainly animal group behaviour [17] and population connectivity [18,19]. Furthermore, the effects of, for example, climate change on the population dynamics of migrants have been quantified with network theory methods [20–22]. However, so far, full, population-level migratory flow on a network of breeding, stopover and wintering sites (nodes) that are linked by seasonal long-distance movements (links) has rarely been modelled [23,24]. Such generalized migration models could be used to study (i) migration on the level of populations or species, (ii) the effects of climate and habitat change on migrants and (iii) disease spread phenomena.

The spread characteristics of spatial networks [25] have been extensively studied for human transportation networks like disease spread on the aviation network [26–28] or bioinvasion on the global cargo ship network [29–31]. Those networks had some very highly connected nodes and groups of tightly interconnected nodes indicating robustness to random topological changes (i.e. deletion of nodes or links) and large potential for fast spread [25]. Most studies have focused on the topological characteristics of transportation networks only, but for describing animal migration as a network, transition links need to be directed, weighted and, most importantly, include timing [14,16,22]. The latter accounts for the observation that the structure of most networks is constantly changing, with link connections evolving in time or existing only intermittently [32,33]. For animal migration networks, the explicit inclusion of timing is especially important, as migration is an inherently seasonal process.

Based on our previously published conceptual ideas [14], we present a periodic Markov model on a network that phenomenologically describes the process of animal migration in time and space. It is based on seasonally changing transition rates, reflecting migratory displacements. We provide a solution for seasonally stable densities of animals in all nodes that represent breeding, stopover and wintering sites. The model was parametrized with two example datasets: GPS and satellite telemetry tracks of white storks (*Ciconia ciconia*) and greater white-fronted geese (*Anser a. albifrons*), and seasonally changing topological and network flow properties were extracted. Node degrees, shortest path distances and motif distributions revealed that the example migration networks were highly connected and efficient for fast transport during migration, but rather stationary during breeding. Thus, our migration network model conserves important properties of the migration process, rendering it highly useful to formalize the spatio-temporal movement of complete animal populations and study the effects of habitat and climate change and disease spread.

### 2. Material and methods

#### 2.1. Migration model

We propose modelling animal migration as a non-homogeneous, periodic Markov process on a geographically embedded network. The nodes of the network are discrete areas of breeding, resting and wintering, and links represent seasonal transition events—movement between the different areas [14]. Thus, we assume there is a network of *n* nodes and a cyclic process that evolves on it (see schematic presentation of a simple version of the migration model in figure 1). The model describes the time-dependent densities *N _{i}*(

*t*) of animals in area

*i*= 1⋯

*n*at time

*t*, determined by the difference of outgoing and ingoing flows

*J*(

_{ji}*t*) and

*J*(

_{ij}*t*) between any two connected regions

*j*and

*i*. One can generally formulate this space discretized process as a time continuous Markov process

*J*(

_{ij}*t*) =

*r*(

_{ij}*τ*)

*N*(

_{j}*t*), with

*r*(

_{ij}*τ*) being the transition rate from node

*j*to

*i*at season

*τ*with $\tau =tmodT$ (figure 1

*a*). That is, the transition rates

*r*(

_{ij}*t*) =

*r*(

_{ij}*t*+

*T*) are assumed to be positive, non-homogeneous (i.e. time-dependent) and cyclic in time (period

*T*= 2

*π*). Time units are adapted to cyclic statistics with 2

*π*corresponding to 1 year. Note that by summation over all nodes (equation (2.1)), one can show that the total animal population size remains conserved, $\sum _{i}{N}_{i}(t)=1.$ If the process is stationary, the densities converge into seasonal stationary states ${N}_{i}^{\ast}(\tau )$.

As a first approximation, we consider the timing of movements between sites as normally distributed, thus the transition rates *r _{ij}*(

*τ*) become von Mises functions (figure 1

*b*)

*φ*is the main transition phase (season),

_{ij}*σ*characterises the length of the major transition period and

_{ij}*ω*is the maximum transition intensity (amplitude) that is typically realized during

_{ij}*τ*≈

*φ*. This functional form of

_{ij}*r*(

_{ij}*τ*) also describes the extreme cases of the

*δ*-function (

*σ*→ 0), when all animal movements on a given connection occur at the same time instance, and the uniform distribution (

_{ij}*σ*→ ∞), when the transition rates are time-independent.

_{ij}#### 2.2. Numerical implementation as a time-discrete process

To implement this process numerically, we formulate a time-discrete version of the continuous process (equation (2.1)). For this we split each year into *m* smaller time intervals of equal length Δ*τ* (*T* = *m* · Δ*τ*). The density *N _{i}*(

*t*+ Δ

*τ*) of animals at site

*i*in the next time interval is then given as the sum of the animals that migrate into the site plus the animals that stay in the site (and do not move in this time interval). This can be derived from equation (2.1) into a Markov chain equation [34]

*a*

*b*

*p*(

_{ij}*τ*;Δ

*τ*) denotes the transition probability to make a movement from node

*j*to

*i*at season

*τ*within the finite time interval of length Δ

*τ*and $\sum _{j\ne i}{p}_{ij}(\tau ;\mathrm{\Delta}\tau )+{p}_{ii}(\tau ;\mathrm{\Delta}\tau )=1$. The transition probabilities

*p*(

_{ij}*τ*;Δ

*τ*) of the time-discrete process (equation (2.3

*a*)) are related to the transition rates

*r*(

_{ij}*τ*) of the time continuous process (equation (2.1)) by an exponential transformation [34]

*P*(

*τ*;Δ

*τ*) = (

*p*)

_{ij}*represents the*

_{ij}*n*×

*n*matrix of transition probabilities and

*R*(

*τ*) = (

*r*)

_{ij}*the corresponding*

_{ij}*n*×

*n*matrix of transition rates, with diagonal elements defined as ${r}_{ii}=-\sum _{j\ne i}{r}_{ij}(\tau )$. Because of the monotone transformation by an exponential (equation (2.4)), the transition probabilities

*p*(

_{ij}*τ*;Δ

*τ*) are naturally related to the von Mises function (equation (2.2)).

To solve for the stationary distribution of densities in all sites, we plug the (empirically determined) transition rates (equation (2.2)) into equation (2.4) and determine the instantaneous transition probability for each discrete time interval, giving rise to a set of *m* probability matrices *P*(*τ*;Δ*τ*). Next, we calculate the Poincaré map *S*(*τ*) from one year to the next at each time step *τ* by matrix multiplication of the *m* subsequent transition matrices

*m*matrices

*S*(

*τ*) is a regular stochastic matrix and by the Perron–Frobenius theorem has a dominant eigenvalue

*λ*(

*τ*) = 1 [35]. Then, the seasonal stationary distribution of densities ${N}_{i}^{\ast}(\tau )$ in sites

*i*at season

*τ*is obtained as the corresponding right eigenvector of

*S*(

*τ*),

#### 2.3. Simple model example

How our model formalizes the population migration process and how seasonal stable densities are driven by network topology and the transition parameters ${\phi}_{ij}$, *σ _{ij}* and

*ω*can be followed in a simple example (figure 1

_{ij}*c*,

*d*). There are four sites (nodes) that are connected by von Mises seasonal transition rates (see figure 1 for parameters) with Δ

*τ*being about one week (

*m*= 48). In the summer months, all animals are in node 1. During autumn migration they move via the stopover site, node 2, to the wintering sites, nodes 3 and 4. Because the transition from node 2 to 3 starts somewhat earlier than from node 2 to 4 (

*φ*

_{32 }<

*φ*

_{42}) winter densities are higher in node 3. During spring migration, animals from node 4 have to migrate via the stopover site, node 2, to the summer region, node 1, whereas animals from node 3 can move via node 2 or directly to node 1. Therefore, node 3 empties earlier in spring than node 4 and the stopover site, node 2, is less used than in autumn. This example shows how we were able to quantify the complex spatio-temporal flow of migratory animals on a seasonal network with only seven model parameters.

#### 2.4. Example migration tracks

To examine our model's ability to portray real stepping stone migration patterns [5], we used a set of quality filtered [36] satellite telemetry and GPS tracks of white storks and greater white-fronted geese. Both are large, long-distance migrants that use several stopover sites along their migration route [37]. We selected full-year trajectories of adult birds because our model requires data that are evenly distributed throughout a full calendar year. For the white stork this yielded 32 yearly tracks of 12 adult individuals [38] and for the white-fronted goose data encompassed 8 yearly tracks of seven adult individuals [39].

#### 2.5. Determination of network nodes

For both datasets, the breeding, wintering and resting areas, i.e. the nodes of the network, had to be determined. This could generally be done using observations or environmental data, but for simplicity, we here used an approach based solely on the migration trajectories. For each track, we first selected one best daily position (based on satellite telemetry location classes and GPS horizontal accuracy, respectively) of the bird at the end of each day and then approximated flight velocities towards the previous and next day's positions (figure 2). The flight velocities showed a distinct bimodal behaviour, of either fast migratory movement or resting, allowing the determination of resting positions based on a simple threshold value.

Best daily positions for which both velocities were less than 50 km d^{−1} for the white storks [40] and 20 km d^{−1} for the white-fronted geese [41] were selected as resting localization (see in relation to data in figure 2). Those positions were then hierarchically clustered by spatial distance and grouped at distances less than 750 km, which was the average distance between migration stopover sites. To account for outliers in our dataset, only clusters that contained more than three resting positions were selected as typical breeding, resting and wintering regions of the considered species. The selected parameters and thresholds determine the resulting network topology and should, therefore, be carefully selected based on data quality, the research question and previous knowledge of the species and their migration behaviour. Per species, we designed a time-cumulative network, only containing the links between those pairs of nodes between which at least one bird had moved, i.e. using the two sites successively.

#### 2.6. Determination of seasonal transition rates

Once the network nodes and links were established, we determined the time-continuous transition rates along these connections. For this, we divided the year into *m* discrete time intervals of equal length Δ*τ* that are sufficiently small so that the transition rates can be assumed to remain constant within each interval. Because of the limited time resolution of the datasets, we selected Δ*τ* to be about one week (Δ*τ* = 0.13 if 1 year = 2*π*), giving rise to *m* = 48 seasonal time intervals. For all seasonal time intervals, we counted the numbers of animals in each of the above defined resting regions, yielding the empirical, seasonal bird densities *N _{j}*(

*τ*) in the network nodes. In a similar way, we counted the seasonal transition flows

*J*(

_{ij}*τ*,

*τ*+ Δ

*τ*) between all connected pairs of nodes for all time intervals. Based on the obtained numbers, we approximated the seasonal transition rates for each discrete time interval as

*φ*, the circular standard deviation of the transition times

_{ij}*σ*and the maximum transition intensity

_{ij}*ω*(see electronic supplementary material, for equations). Thus, we obtained a smooth model transition rate

_{ij}*r*(

_{ij}*τ*) for any season

*τ*(equation (2.2)). Finally, using the migration model (equations (2.4)–(2.6)), we obtained seasonal stationary bird densities ${N}_{i}^{\ast}(\tau )$ for each time interval

*τ*and node

*i*of the two migration networks. We provide R code for determination of the migration network model in the electronic supplementary material.

#### 2.7. Characterization of migration network models

To evaluate whether our parametrized models conform to the main characteristic of migration systems, namely shared, directed seasonal movement between breeding and wintering regions, we assessed the estimated population densities and a selection of network measures for the two cumulative networks as well as the two sets of seasonal time-dependent networks. Season-specific migration network topologies were derived such that a pair of nodes was connected by a directed link from *j* to *i* whenever *p _{ij}*(

*τ*;Δ

*t*) > 0.05. Those were then merged into time-cumulative networks by overlap.

Of several standard network measures that can help to further characterize migration networks [16], we here focus on one property of the cumulative network that can easily be calculated, the shortest directed topological path distance. This is the minimum number of links (and stopover sites) that a migrating animal has to traverse to move between two nodes. Comparatively short paths between breeding and wintering regions indicate migration. Important network measures characterizing high possible movement between nodes are node degree (number of links that start or end in a node) and the size of the giant component (size of the largest connected sub-network). For our purposes, the seasonal changes of those measures are most indicative of migration, namely times of high migration activity versus times of breeding and winter, thus we determined mean degree and giant component size of the season-specific networks.

For further confirmation that our networks allow quick directed migration flow, we determined their motif distributions [43], i.e. the distribution of differently connected triplets. If such three-node sub-networks are, for example, connected by two directed links to allow flow through the middle node, this is an indication of a transportation network, as our migration networks should be.

We determined each node's degree and betweenness centrality, which is defined as the number of topologically shortest paths in the network that pass the respective node [44]. This measure characterizes how much each node contributed to network connectivity, flow and clustering. Further, we determined the cumulative density and the staying time after the highest densities, which indicate the nodes that were used often and by a large number of birds. For more detailed flow characterization, we determined the first passage times [35] between two different nodes, as well as return times to the same node. As those values varied with starting season *τ*, we selected the first passage and return times starting in the season of highest ${N}_{i}^{\ast}(\tau )$ to provide minimum values.

#### 2.8. Model validation and sensitivity

To validate the model's appropriateness for mapping animal migration, we explored how well empirical and fitted non-zero transition rates ${\stackrel{~}{r}}_{ij}(\tau )$ and *r _{ij}*(

*τ*) matched using paired Wilcoxon's tests. Comparisons of the distribution of the number of birds in the data with the modelled, stable seasonal densities ${N}_{i}^{\ast}(\tau )$ were performed using

*G*-tests [45].

Sensitivity of model outcomes ${N}_{i}^{\ast}(\tau )$ to changes of the phases *φ _{ij}*, widths

*σ*and amplitudes

_{ij}*ω*of single transition rates was determined as follows. For single transition rates, in turn, one of the three parameters was varied and model outcomes were compared with the original ones by

_{ij}*G*-tests. Transition rate amplitudes were modified multiplicatively by values from 0.01 to 100, phases were shifted additively by values out of [ − 1.58, + 1.58] (±3 months) and transition rate widths were increased or decreased by values of the same interval (we forced a truncation at a minimum value of 1 day = 0.017). The outcomes were considered sufficiently close to the original parametrizations' results if more than 67% of the

*G*-tests for the distributions in each season

*τ*had a

*p*> 0.1.

### 3. Results

#### 3.1. Migration network models

The determined cumulative migration networks (figure 3*a*,*c*) of breeding, resting and wintering regions showed much individual variability for both species. The network of stork migration consisted of *n = *24 nodes and *l = *96 directed links and the goose network had *n *= 17 nodes and *l *= 53 links, whereas individuals of both species used only up to nine sites per year.

The seasonal flow of birds in the networks (figure 3*b*,*d*) indicated cyclic movement in time of all individuals. At some times of the year, all birds remained in a few nodes only, whereas during migration much movement happened. Furthermore, the flow of birds was not just linear along a number of nodes, but a lot of branching occurred, especially in the stork network. Some of the nodes were frequented by a large number of birds and for a long time, but others seem to have simply been short-time resting areas. Times when birds actually used a node seemed to overlap well with model density estimations (small red lines in figure 3*c*,*d*). Six season-specific, cumulated sub-networks of the stork model (figure 4) outline how the model generally maps migratory flow characteristics.

#### 3.2. Characterization of migration network models

In comparison to the sizes of the networks, the mean shortest (directed) path lengths were very small (2.47 for storks, 2.46 for geese). When looking at shortest paths between the breeding and wintering sites, first migration properties of our networks appeared: between nodes 2 (Middle European breeding region) and 8 (Sudan wintering and staging area) of the stork network only one link had to be crossed in either direction, indicating direct, possibly fast migration. Movement between node 2 and the southernmost node 24 (South Africa) could be done through three links in autumn, but by one direct link in fast spring migration. Similarly, in the goose network, the shorted path between node 1 (Middle European wintering region) and node 11 (breeding region) had length 3 in spring (many stopovers), but only length 1 in autumn (quick directed back migration). Mean degree and the size of the giant component of the time changing networks (figure 5*a*) further illustrated times of migration movement and when birds were moving slowly.

Betweenness centrality (electronic supplementary material, table S1) revealed that for the stork network, nodes 2 (Middle European breeding region) and 8 (Sudan wintering and staging area) were especially important, whereas nodes 3 (French breeding area) and 20 (staging area close to Lake Malawi) were not central at all. For the goose network, nodes 1 and 4 were of high importance and nodes 2 and 12 not. Cumulative densities (in units of density × time intervals; range = [0; 48], 48 indicating that all individuals were in the same node all year) in highly central nodes were mostly large (12.66 for stork node 2, 12.66 for goose node 1; note that the equality of the two values is coincidental), but staying times were more intermediate (2.23 = 130 days for stork node 2, 2.50 = 145 days for goose node 1). Longer staying times were mostly obtained for less connected nodes, as calculated densities were stuck there due to small transition rates obtained from only one bird in our dataset taking that route.

The motif distributions confirm that there are a large number of sub-graphs that represent transit movement (sub-networks 3, 4), which is an important property of migration. Additionally, patterns of branching character (sub-graphs 7, 9 and 10) are overrepresented in comparison to random graphs of similar structure. It can be observed that the motif distribution (figure 5*b*) falls into the second superfamily [46] which also includes the signal-transduction and neurons networks of simpler organisms.

First passage times after maximum starting density (figure 5*c*,*d*) indicated which of the nodes were reached on average more frequently and quicker than others. They were again nodes 2 and 8 for the stork network and nodes 1 and 4 for the geese. Return times (diagonals in figure 5*c*,*d*) were generally shorter than first passage times, indicating possible backwards movements due to long transition rate distributions (i.e. large *σ _{ij}*).

#### 3.3. Model verification and sensitivity analysis

The counted and the fitted transition rates ${\stackrel{~}{r}}_{ij}(\tau )$ and *r _{ij}*(

*τ*) were very similar (Wilcoxon tests:

*W*

_{s}= 9040,

*p*

_{s}= 0.23;

*W*

_{g}= 1870,

*p*

_{g}= 0.29) with average deviations of only 10% (± 14%) for the stork network and 7% (±13%) for the goose network. Also resulting estimated densities ${N}_{i}^{\ast}(\tau )$ were very similar to the numbers of counted birds in each region at the respective seasons (

*G*-tests: ${\u27e8G\u27e9}_{\mathrm{s}}=15.77$ (

*p*= 0.76) and ${\u27e8G\u27e9}_{\mathrm{g}}=11.01$ (

*p*= 0.96); figure 3

*b*,

*d*), with average differences of 0.65 (storks) and 0.67 (geese) individuals (of

*N*= 32 and

*N*= 8, respectively). Tests of 87.5% and 100% of the seasonal stork and goose models versus densities, respectively, had

*p*> 0.1.

Both migration networks were robust to single changes in the parametrizations of the phases *φ _{ij}* and amplitudes

*ω*of the transition rates, with few exceptions only. For the stork network, only a decrease of

_{ij}*φ*

_{52}by more than 0.48 (= 28 days), an increase of

*ω*

_{62}(by more than *25) or a decrease of

*ω*

_{26}(by more than *0.03) led to significant model alterations. For the goose network, all phases could be decreased or increased up to three months without significant change in ${N}_{i}^{\ast}(\tau )$, but decreased amplitudes

*ω*

_{16}and

*ω*

_{15,17}(both by more than *0.03) altered model outcomes significantly.

On the contrary, model outcomes were relatively sensitive to changes in transition rate widths *σ _{ij}*. The subtraction of small values from single

*σ*often led to

_{ij}*σ*≈ 0, which is one possible cause of its greater influence on model behaviour. For the stork network, addition to

_{ij}*σ*

_{52}(+0.241 = 14 days) changed ${N}_{i}^{\ast}(\tau )$ notably, as well as addition to

*σ*

_{81}(+1.102 = 64 days) and

*σ*

_{17,15}(+1.067 = 62 days) for the goose network. Subtraction of very small values (2–14 days) to most

*σ*had great impact on model outcome. However, a few

_{ij}*σ*were more robust to changes, e.g.

_{ij}*σ*

_{8,12},

*σ*

_{13,9}and

*σ*

_{10,9}for the stork network and

*σ*

_{13}for the goose network (from −0.70 to −1.03 = 40–60 days). These sensitivities indicate that the ranges of inter-node transition times strongly affect migrant seasonal abundance, in the model as well as in migration itself.

### 4. Discussion

Combining methods from non-homogeneous Markov modelling and network theory, we have developed a spatio-temporal periodic animal migration flow model. From seasonal transition rates between breeding, stopover and wintering sites, we derived stable seasonal densities of animals, thus generalizing population migration flow. The model was parametrized with high-quality, year-round tracking data of two species of long-distance migrants. Cumulative network characteristics of the two resulting networks and flow processes are very similar and show properties of a seasonal migration system. Thus, our model phenomenologically maps migration movement in a realistic way and can be used for further applications to explore the effect of migration on traversed ecosystems and vice versa.

In comparison to already existing migration network models, we have solved one of their main limitations, namely the use of suitable empirical data (i.e. tracking data) for migratory connectivity and movement quantification [20], rather than expert knowledge [21] or ringing and count data [23]. Thus, also assumptions about migration flow characteristics like optimal flow or the general flux law based on habitat attractiveness [21,22] are not necessary. Other than previous systematic studies on network patterns and effects of topological changes [15,22], our networks are very close to the available datasets and allow for predictions of realistic change effects or spread characteristics in the modelled systems.

Different from the existing explicit population network models [20–23], our model is a pure movement model that obeys the conservation law, i.e. we neglect population growth and consider total population size constant during the timescale of the analysis. For our datasets of adult storks and geese, this assumption is sensible in relatively stable environments, as those species are long-lived. However, if individuals are short-lived and population growth is strong and not similar for individuals using different migration routes or timing, resulting movement patterns might be unrealistic in the long run. If comparing model outcomes of different, shorter time periods, one might be able to detect such changes in long-term migratory dynamics like population size, stopover usage or timing.

Further assumptions of our model are the Markov property, which implies that birds use no memory of their previous whereabouts but only season and their present state for the decision regarding where to go next or if to stay, and the independence of the transition rates from the present densities in the nodes, thus we disregard any swarming or density dependence effects. Finally, the decision for migratory departure from breeding, stopover and wintering sites has been proposed to depend on environmental characteristics like the photoperiod and most importantly food availability and uptake [3,47], that might vary between years. Such habitat quality dependence is not explicitly accounted for in our model, but is implicitly included in the time dependence of the transition probabilities.

All assumptions were made to introduce our basic, phenomenological network model of stable migration movement. It can be readily extended with additional terms of past densities, immigration, mortality, density dependence or environmental factors. For example, an additional sink/source node can be added for integration of vital rates. However, already without any additions, the novel issue that network structure changes periodically with time makes processes on our interchanging network model highly complex.

By translating bird migration into a temporal network, we have gained access to the huge arsenal of technical and statistical tools developed by network theory. We showed that application of network indices indeed allows one to obtain useful insights into the spatio-temporal flow of birds during their migration and into the small- and large-scale organization, connectivity and its temporal dependence, which would not be straightforward to obtain with other methods. Namely, the network characteristics of our cumulative stork and goose migration models were very similar to other transportation networks, especially the aviation and cargo ship networks. Like those, our migration networks are inherently non-planar, show properties of transit and have small shortest paths. However, differently from human transportation networks, our described migration flow is strictly seasonal and there are characteristic pauses during breeding, winter and in stopovers regions that differ strongly by season and species. Surprisingly, for storks it is autumn migration, but for geese spring migration that is interspersed with several stopovers and makes the movement flow irregular rather than uniform, likely driven by environmental conditions [41].

Because GPS tracks of migratory animals become more and more ubiquitous and of higher quality [11], it will become possible to parametrize our model in more detail and possibly for different, interacting populations or species. This will allow more accurate network characterization and the evaluation of processes on this network, like cross-continental interspecific disease spread. Moreover, there are other methods to determine network nodes, i.e. breeding, stopover and wintering sites from movement trajectories [48], and it should be decided by species and population which one is most applicable. Our analyses of model sensitivity to changes in the transition rates on the network links revealed that the duration of possible transition between sites strongly affected network connectivity, whereas peak transition time and magnitude could be varied without great influence. This indicates that migrants with high plasticity in migration timing are more robust to changes.

Specific future applications of our model are plentiful, particularly in the context of trophic and transport effects of migration. By quantifying the number of individuals present at stopover sites during different seasons, one can assess the potential benefits and costs of migrating populations on local animal and plant communities. For white-fronted and other goose species this is particularly helpful in quantifying the pressure on arable fields, as well as the impact of population management at multiple spatial scales [49]. For white storks, it can help to understand where and when populations are most vulnerable to mortality and thus where and when management measures should be installed [50]. These examples can be extended to a wide range of other migrants. Obvious next steps are the combination with epidemiological models to simulate the spread of zoonotic diseases like avian influenza or West Nile virus [51], or to simulate the resilience of migratory animal populations to changes in habitat suitability at important stopover sites resulting from anthropogenic or climate-driven land-use changes [52].

In conclusion, our model provides an empirically informed quantitative tool to formalize sets of individual migration tracks into spatio-temporal dynamics of migratory populations. The periodically changing network process is a new type of complex model and we invite further theoretical studies to examine and understand its dynamic properties. We propose using our cyclic description of migration patterns in combination with the large suite of mechanistic migration models already available [53] to study different aspects of animal migration. By using recently available larger datasets and possibly more migratory species in combination, our model will show its full potential to gain important new insights into a range of current and pressing topics in the field of animal migration.

### Ethics

Permission for the equipment of white storks with telemetry tags was issued by the ‘Landesamt für Naturschutz Sachsen-Anhalt' to M.K. White-fronted geese were tagged under the umbrella of Wageningen Environmental Research (Alterra) in the Netherlands, where at the time equipment of animals with telemetry tags was not yet considered an animal experiment that needed approval.

### Data accessibility

The GPS and satellite telemetry tracks used in the study have been uploaded to Movebank (www.movebank.org), available in the studies ‘White stork full year tracks 1994–2006 MPIO' and ‘White-fronted goose full year tracks 2006–2010 Alterra IWWR', and published in the Movebank Data Repository (doi:10.5441/001/1.k29d81dh) [38] and (doi:10.5441/001/1.kk38017f) [39].

### Authors' contributions

A.K. and B.B. designed the study and developed the model. H.K. and M.K. provided the tracking data and background information. A.K. and E.K. wrote the manuscript that was subsequently revised by all other authors.

### Competing interests

We have no competing interests.

### Funding

This work was funded by the German Volkswagen-Stiftung (grant Az: 82723) and the German BMBF (grant 03F0449A).

## Acknowledgements

We are grateful to the Vogelschutz-Komitee e.V. (VsK, Hamburg), Wageningen Environmental Research (Alterra) and the Dutch Society of Goose Catchers for financial and technical support for goose catching and tracking. We want to thank D. Syga, R. Tönjes and A. Ryabov for helpful discussion on the modelling approach and B.A. Nolet for comments on an earlier version of this manuscript.

### References

- 1
Milner-Gulland EJ, Fryxell JM, Sinclair ARE . 2011**Anima migration—a synthesis**. Oxford, UK: Oxford University Press. Crossref, Google Scholar - 2
Bauer S, Hoye BJ . 2014 Migratory animals couple biodiversity and ecosystem functioning worldwide.**Science**, 1242552. (doi:10.1126/science.1242552) Crossref, PubMed, Web of Science, Google Scholar**344** - 3
Berthold P, Terrill SB . 1991 Recent advances in studies of bird migration.**Ann. Rev. Ecol. Syst.**, 357–378. (doi:10.1146/annurev.es.22.110191.002041) Crossref, Google Scholar**22** - 4
Drent RH, Eichhorn G, Flagstad A, van der Graaf AJ, Litvin KE, Stahl J . 2007 Migratory connectivity in Arctic geese: spring stopovers are the weak links in meeting targets for breeding.**J. Ornithol.**, 501–514. (doi:10.1007/s10336-007-0223-4) Crossref, Web of Science, Google Scholar**148** - 5
Hedenström A, Alerstam T . 1997 Optimum fuel loads in migratory birds: distinguishing between time and energy minimization.**J. Theor. Biol.**, 227–234. (doi:10.1006/jtbi.1997.0505) Crossref, PubMed, Web of Science, Google Scholar**189** - 6
Wiltschko R, Wiltschko W . 2003 Avian navigation: from historical to modern concepts.**Anim. Behav.**, 257–272. (doi:10.1006/anbe.2003.2054) Crossref, Web of Science, Google Scholar**65** - 7
Van Turnhout CAM, Foppen RPB, Leuven R, Van Strien A, Siepel H . 2010 Life-history and ecological correlates of population change in Dutch breeding birds.**Biol. Conserv.**, 173–181. (doi:10.1016/j.biocon.2009.09.023) Crossref, Web of Science, Google Scholar**143** - 8
Saino N 2011 Climate warming, ecological mismatch at arrival and population decline in migratory birds.**Proc. R. Soc. B**, 835–842. (doi:10.1098/rspb.2010.1778) Link, Web of Science, Google Scholar**278** - 9
Altizer S, Bartel R, Han BA . 2011 Animal migration and infectious disease risk.**Science**, 296–302. (doi:10.1126/science.1194694) Crossref, PubMed, Web of Science, Google Scholar**331** - 10
Verhagen JH, Herfst S, Fouchier RAM . 2015 How a virus travels the world.**Science**, 616–617. (doi:10.1126/science.aaa6724) Crossref, PubMed, Web of Science, Google Scholar**347** - 11
Bridge ES 2011 Technology on the move: recent and forthcoming innovations for tracking migratory birds.**Bioscience**, 689–698. (doi:10.1525/bio.2011.61.9.7) Crossref, Web of Science, Google Scholar**61** - 12
Vardanis Y, Klaassen RHG, Strandberg R, Alerstam T . 2011 Individuality in bird migration: routes and timing.**Biol. Lett.**, 502–505. (doi:10.1098/rsbl.2010.1180) Link, Web of Science, Google Scholar**7** - 13
Ambrosini R 2014 Modelling the progression of bird migration with conditional autoregressive models applied to ringing data.**PLoS ONE**, e102440. (doi:10.1371/journal.pone.0102440) Crossref, PubMed, Web of Science, Google Scholar**9** - 14
Kölzsch A, Blasius B . 2008 Theoretical approaches to bird migration. The white stork as a case study.**Eur. Phys. J. Spec. Top.**, 191–208. (doi:10.1140/epjst/e2008-00641-y) Crossref, Web of Science, Google Scholar**157** - 15
Taylor CM, Norris DR . 2010 Population dynamics in migratory networks.**Theor. Ecol.**, 65–73. (doi:10.1007/s12080-009-0054-4) Crossref, Web of Science, Google Scholar**3** - 16
Jacoby DMP, Freeman R . 2016 Emerging network-based tools in movement ecology.**Trends Ecol. Evol.**, 301–314. (doi:10.1016/j.tree.2016.01.011) Crossref, PubMed, Web of Science, Google Scholar**31** - 17
Farine DR, Whitehead H . 2015 Constructing, conducting and interpreting animal social network analysis.**J. Anim. Ecol.**, 1144–1163. (doi:10.1111/1365-2656.12418) Crossref, PubMed, Web of Science, Google Scholar**84** - 18
Urban D, Keitt T . 2001 Landscape connectivity: a graph-theoretic perspective.**Ecology**, 1205–1218. (doi:10.1890/0012-9658(2001)082[1205:LCAGTP]2.0.CO;2) Crossref, Web of Science, Google Scholar**82** - 19
Wells K, Dolich T, Wahl J, O'Hara RB . 2013 Spatio-temporal dynamics in waterbirds during the non-breeding season: effects of local movements, migration and weather are monthly, not yearly.**Basic Appl. Ecol.**, 523–531. (doi:10.1016/j.baae.2013.07.001) Crossref, Web of Science, Google Scholar**14** - 20
Hostetler JA, Sillett TS, Marra PP . 2015 Full-annual-cycle population models for migratory birds.**Auk**, 433–449. (doi:10.1642/AUK-14-211.1) Crossref, Web of Science, Google Scholar**132** - 21
Iwamura T, Possingham H, Chades I, Minton C, Murray N, Rogers D, Treml E, Fuller R . 2013 Migratory connectivity magnifies the consequences of habitat loss from sea-level rise for shorebird populations.**Proc. R. Soc. B**, 20130325. (doi:10.1098/rspb.2013.0325) Link, Web of Science, Google Scholar**280** - 22
Taylor CM, Laughlin AJ, Hall RJ . 2016 The response of migratory populations to phenological change: a migratory flow network modelling approach.**J. Anim. Ecol.**, 648–659. (doi:10.1111/1365-2656.12494) Crossref, PubMed, Web of Science, Google Scholar**85** - 23
Jain N, Dilkina B . 2015 Coarse models for bird migrations using clustering and non-stationary Markov chains. InConference on Artificial Intelligence (AAAI) ,Austin, Texas ,25–30 January , pp. 63–68. Palo Alto, CA: The AAAI Press. Google Scholar - 24
Sheldon D, Elmohamed MAS, Kozen D . 2007 Collective inference on Markov models for modeling bird migration. In**Neural information processing systems conference**(ed.Platt JC ), pp. 1321–1328, Cambridge, MA: MIT Press. Google Scholar - 25
Barthelemy M . 2011 Spatial networks.**Phys. Rep. Rev. Sec. Phys. Lett.**, 1–101. Web of Science, Google Scholar**499** - 26
Colizza V, Barrat A, Barthelemy M, Vespignani A . 2006 The role of the airline transportation network in the prediction and predictability of global epidemics.**Proc. Natl Acad. Sci. USA**, 2015–2020. (doi:10.1073/pnas.0510525103) Crossref, PubMed, Web of Science, Google Scholar**103** - 27
Guimera R, Mossa S, Turtschi A, Amaral LAN . 2005 The worldwide air transportation network: anomalous centrality, community structure, and cities' global roles.**Proc. Natl Acad. Sci. USA**, 7794–7799. (doi:10.1073/pnas.0407994102) Crossref, PubMed, Web of Science, Google Scholar**102** - 28
Hufnagel L, Brockmann D, Geisel T . 2004 Forecast and control of epidemics in a globalized world.**Proc. Natl Acad. Sci. USA**, 15 124–15 129. (doi:10.1073/pnas.0308344101) Crossref, Web of Science, Google Scholar**101** - 29
Kaluza P, Kölzsch A, Gastner MT, Blasius B . 2010 The complex network of global cargo ship movements.**J. R. Soc. Interface**, 1093–1103. (doi:10.1098/rsif.2009.0495) Link, Web of Science, Google Scholar**7** - 30
Kölzsch A, Blasius B . 2011 Indications of marine bioinvasion from network theory.**Eur. Phys. J. B**, 601–612. (doi:10.1140/epjb/e2011-20228-5) Crossref, Web of Science, Google Scholar**84** - 31
Seebens H, Gastner MT, Blasius B . 2013 The risk of marine bioinvasion caused by global shipping.**Ecol. Lett.**, 782–790. (doi:10.1111/ele.12111) Crossref, PubMed, Web of Science, Google Scholar**16** - 32
Gross T, Blasius B . 2008 Adaptive coevolutionary networks: a review.**J. R. Soc. Interface**, 259–271. (doi:10.1098/rsif.2007.1229) Link, Web of Science, Google Scholar**5** - 33
Li A, Cornelius SP, Liu YY, Wang L, Barabasi AL . 2017 The fundamental advantages of temporal networks.**Science**, 1042–1046. (doi:10.1126/science.aai7488) Crossref, PubMed, Web of Science, Google Scholar**358** - 34
Logofet DO, Lesnaya EV . 2000 The mathematics of Markov models: what Markov chains can really predict in forest succession.**Ecol. Modell.**, 285–298. (doi:10.1016/S0304-3800(00)00269-6) Crossref, Web of Science, Google Scholar**126** - 35
- 36
Austin D, McMillan JI . 2003 A three-stage algorithm for filtering erroneous Argos satellite locations.**Mar. Mamm. Sci.**, 371–383. (doi:10.1111/j.1748-7692.2003.tb01115.x) Crossref, Web of Science, Google Scholar**19** - 37
Beekman JH, Nolet BA, Klaassen M . 2002 Skipping swans: fuelling rates and wind conditions determine differential use of migratory stopover sites of Bewick's swans*Cygnus bewickii*.**Ardea**, 437–460. Web of Science, Google Scholar**90** - 38
Kaatz M. 2018Data from: A periodic Markov model to formalise animal migration on a network [white stork data] .**Movebank Data Repository**. (doi:10.5441/001/1.k29d81dh) Google Scholar - 39
Kruckenberg H, Müskens GJDM, Ebbinge BS. 2018Data from: A periodic Markov model to formalise animal migration on a network [white-fronted goose data] .**Movebank Data Repository**. (doi:10.5441/001/1.kk38017f) Google Scholar - 40
Gerkmann B, Riede K . 2005 Use of satellite telemetry and remote sensing data to identify important habitats of migratory birds (*Ciconia ciconia*). In**African biodiversity: molecules, organisms, ecosystems**(edsSinclair BJ, Humber BA, Lampe K-H ), pp. 261–269. Bonn, Germany: Springer. Crossref, Google Scholar - 41
Kölzsch A, Müskens GJDM, Kruckenberg H, Glazov P, Weinzierl R, Nolet BA, Wikelski M . 2016 Towards a new understanding of migration timing: slower spring than autumn migration in geese reflects different decision rules for stopover use and departure.**Oikos**, 1496–1507. (doi:10.1111/oik.03121) Crossref, Web of Science, Google Scholar**125** - 42
- 43
Milo R, Shen-Orr S, Itzkovitz S, Kashtan N, Chklovskii D, Alon U . 2002 Network motifs: simple building blocks of complex networks.**Science**, 824–827. (doi:10.1126/science.298.5594.824) Crossref, PubMed, Web of Science, Google Scholar**298** - 44
Freeman L . 1997 A set of measures of centrality based upon betweenness.**Sociometry**, 35–41. (doi:10.2307/3033543) Crossref, Google Scholar**40** - 45
Sokal RR, Rohlf FJ . 1995**Biometry: the principles and practice of statistics in biological research**. New York, NY: Freemann and Company. Google Scholar - 46
Milo R, Itzkovitz S, Kashtan N, Levitt R, Shen-Orr S, Ayzenshtat I, Sheffer M, Alon U . 2004 Superfamilies of evolved and designed networks.**Science**, 1538–1542. (doi:10.1126/science.1089167) Crossref, PubMed, Web of Science, Google Scholar**303** - 47
Gwinner E, Schwabl H, Schwabl-Bensinger I . 1988 Effects of food deprivation on migratory restlessness and diurnal activity in the garden warbler*Sylvia borin*.**Oecologia**, 321–326. (doi:10.1007/BF00378037) Crossref, PubMed, Web of Science, Google Scholar**77** - 48
Buchin M, Kruckenberg H, Kölzsch A . 2013 Segmenting trajectories by movement states. In**Advances in spatial data handling**(edsTimpf S, Laube P ), pp. 15–25. Berlin, Germany: Springer. Crossref, Google Scholar - 49
Klaassen M, Bauer S, Madsen J, Ingunn T . 2006 Modelling behavioural and fitness consequences of disturbance for geese along their spring flyway.**J. Appl. Ecol.**, 92–100. (doi:10.1111/j.1365-2664.2005.01109.x) Crossref, Web of Science, Google Scholar**43** - 50
Flack A 2016 Costs of migratory decisions: a comparison across eight white stork populations.**Sci. Adv.**, e1500931. (doi:10.1126/sciadv.1500931) Crossref, PubMed, Web of Science, Google Scholar**2** - 51
Reed KD, Meece JK, Henkel JS, Shukla SK . 2003 Birds, migration and emerging zoonoses: West Nile virus, Lyme disease, influenza A and enteropathogens.**Clin. Med. Res.**, 5–12. (doi:10.3121/cmr.1.1.5) Crossref, PubMed, Google Scholar**1** - 52
Wilcove DS, Wikelski M . 2008 Going, going, gone: is animal migration disappearing.**PLoS Biol.**, e188. (doi:10.1371/journal.pbio.0060188) Crossref, PubMed, Web of Science, Google Scholar**6** - 53
Bauer S, Klaassen M . 2013 Mechanistic models of animal migration behaviour their diversity, structure and use.**J. Anim. Ecol.**, 498–508. (doi:10.1111/1365-2656.12054) Crossref, PubMed, Web of Science, Google Scholar**82**