The Influence of Temporal Migration in the Synchronization of Populations^†

Abstracts

A discrete metapopulation model with temporal dependent migration is proposed in order to study the stability of synchronized dynamics. During each time step, we assume that there are two processes involved in the population dynamics: local patch dynamics and migration process between the patches that compose the metapopulation. We obtain an analytical criterion that depends on the local patch dynamics (Lyapunov number) and on the whole migration process. The stability of synchronized dynamics depends on how individuals disperse among the patches.

metapopulation; temporal migration; synchronization

Um modelo metapopulacional com migração dependente do tempo é proposto a fim de estudarmos a estabilidade de trajetórias sincronizadas. Durante cada geração, consi deramos que existem dois processos na dinâmica populacional: a dinâmica local e a mi gra ção entre os sítios que compõem a metapopulação. Obtemos um critério para a ocorrência de sincronização que depende da dinâmica local (número de Lyapunov) e de todo o processomigratório. A estabilidade de trajetórias sincronizadas depende de como os indivíduosmigram entre os sítios.

metapopulação; migração dependente do tempo; sincronização

1 INTRODUCTION

The forms of dispersion in a metapopulation system (populations of single-species that live in fragments called patches) can induce the whole system to multiple behaviors ^{(1) , (3) , (4) , (7) , (13)}. An interesting behavior related to the dispersal process is the synchronized dynamics where the populations in all patches evolve with the same density ⁽¹¹⁾. Its importance lies in the fact that if the whole metapopulation is not synchronized and a local population is extinct, it can be recolonized by individuals that migrate from neighboring patches ("rescue effect"), favoring the population persistence ⁽²⁾. A considerable number of populations that live in distinct regions tend to cycle in synchrony. A well-documented example is the Canadian lynxes that presents synchronized dynamics in its densities fluctuations due to weather conditions ^{(2) , (11)}. Another example is the vole populations in Norway that synchronize due to dispersal processes and birds predation ⁽⁸⁾.

Systems of discrete equations are often used to model metapopulations ^{(1) , (4) , (7) , (13)}. A metapopulation model with patches linked by migration and subjected to external perturbations was considered in ⁽¹⁾. The model is a discrete-time system composed by single species where a constant fraction disperses per generation. Through numerical simulations, it was shown that chaos can prevent global extinctions when coupling is weak. In ⁽⁴⁾ was obtained an analytical result for the stability of synchronized trajectories by considering a model with an arbitrary number of patches linked by dispersal. An analytical result examining a special case of density-dependent dispersal was obtained in ⁽¹³⁾, concluding that this mechanism reduces the stability regions of the synchronous dynamics. Nevertheless, density independent dispersal is observed in the dispersal of insects, while density-dependent dispersal is observed in such widely different invertebrates as locusts, snails and copepods ⁽⁶⁾. In this paper, we present a metapopulation model similar to the ones described in ^{(1) , (4) , (13)}. The main difference is the assumption of temporal dependent migration. This assumption can be used in order to describe the movements of species that move to other areas in different periods due to weather conditions or dependence of foraging resources.

In Section 2 we present the metapopulation model with temporal dependent migration. In Section 3 we analyze the asymptotic local stability of synchronized trajectories and obtain a criterion to its stability based on the calculation of the transversal Lyapunov numbers. In Section 4 we present numerical simulations. Final comments and discussion are done in Section 5.

2 METAPOPULATION MODEL

The metapopulation model consists of d equal patches labeled as 1, 2, ..., d. We assume that the processes of survival and reproduction which compose the local dynamics is described by a map f on [0, ∞) of class C ¹. In the absence of dispersal between patches, the time evolution of the population is given by

where x _t represents the number of individuals at time t.

We assume that a fraction of individuals leaves patch i and disperses to the neighboring patches. We assume that the migration fraction is temporal dependent, that is, it is given by a map on [0, 1] such that m _t + 1 = g(m _t), where m _t is the migration fraction at time t. Thus, the density of individuals that leaves patch i is given by m _t f( ), where denotes the population density in patch i at time t, for all I = 1, ..., d, t = 0, 1, .... Moreover, from individuals that disperse from neighboring patches k, a fraction γ_ik reaches patch i. The dispersal fractions of individuals that migrate among the patches is described by a nonnegative matrix Γ with entries γ_ik, i, k = 1, 2, ..., d. Each γ_ik represents the fraction of individuals coming from patch k that will settle in patch i (Fig. 1). We assume that there is no loss of individuals and the individuals do not return to the original patch, therefore and γ_kk = 0 for all k = 1, ..., d. Taking these into consideration, we can write a system of equations describing the dynamics of the metapopulation by

Figure 1:
Schematic represetation of the migration process. After patch local dynamics, a fraction of individuals m _t leaves patch k at time step t. From these individuals a fraction γ_ik moves to patch i. Thus, the fraction of individuals that leaves patch k and reaches patck i is m _t γ_ik.

The first term in equation (2.2) represents the individuals that did not leave patch i at time t, while the second term is the sum of all contributions of individuals from neighboring patches.

3 SYNCHRONIZATION AND TRANSVERSAL STABILITY

We assume that synchronization is achieved if the population density of all patches is the same, that is, = , for all I = 1, 2, ..., d and t = 0, 1, 2, .... Synchronized solutions of the system (2.2) may not exist if the system lacks some symmetry. Substitution of = in equation (2.2) leads us to the existence of such synchronized solution provided γ_ik. = 1, I = 1, 2, ..., d. Furthermore, he dynamics of each patch in the synchronized state satisfies which is equivalent to equation (2.1), the single patch model equation. In other words the metapopulation synchronizes with the same dynamics of a single isolated patch.

We are interested in studying the local asymptotic stability of the synchronized state, that is, whether orbits that initiate close to the synchronized state will be attracted to it. In order to achieve this goal, we linearize the equation (2.2) around the synchronized trajectory, obtaining

where Δt ∈ R d is the perturbation of the synchronized trajectory, and J( ) is the d × d Jacobian matrix of system (2.2) evaluated at , where = ( , , ..., ) ∈ Rd. Notice that the Jacobian matrix J( ) has its entries given by

Thus, it can be written as

where I is the identity matrix and B = I - Γ.

It is important to notice that the connectivity matrix Γ is doubly stochastic (all rows and columns add up to one). An application of the Perron-Frobenius Theorem ⁽⁹⁾ leads to the fact that λ₀ = 1 is the dominant eigenvalue of Γ. Moreover, its eigenspace is spanned by the vector = [1 1 ... 1]T ∈ R d which correspond to the in-phase state, that is, the diagonal of the phase space .

We assume that Γ is diagonalizable which allows us to express the Jacobian matrix as a diagonal matrix and the local stability of the synchronized state can be analyzed through the diagonal terms. With this assumption, there exists an invertible matrix Q such that

Γ = Qdiag(1, λ1, ..., λd - 1)Q-1.

It allows us to write the Jacobian matrix (3.4) in the following diagonal form

Thus, the synchronized state will be stable if transversal perturbations to the synchronized state shrink to zero. To reach this goal, we define the maximum transversal Lyapunov number, K, by

where P _{t, i} = f' ( ) (1 - m _t + λ_i m _t). Consequently, the transversal perturbation tends to zero if K ( , m ₀) < 1.

Observe that

thus, we can write the maximum transversal Lyapunov number as

where

is the Lyapunov number of f starting at and

is a quantifier that depends on the initial migration rate.

Let ρ be the natural probability measure for the local map f. Let υ be the natural probability measure for map g. Assuming the integrability of ln⁺ |f' (x)| and ln⁺ |1 - λm) with respect to ρ and υ, we can apply the Ergodic Theorem of Birkhoff ⁽⁵⁾ to guarantee the existence and uniqueness ρ-almost every of the limit defining L, and υ -almost every m ₀ of the limit defining Λ and state a criterion for the local asymptotic stability of an attractor in the synchronized invariant state given by

where

and

Notice that L depends on the patch local dynamics while Λ depends on the whole migratory process. It is important to observe that the evolution of the term that corresponds to the value 1 in the Jacobian matrix (3.5) is exactly the Lyapunov number and it gives the behavior of the synchronized trajectory within the phase space diagonal, that is, a periodic trajectory (L < 1) or a chaotic trajectory (L > 1).

In the following, we calculate the quantifier given in (3.13) to different temporal migration rules. In subsection 3.1 we assume that the map g that generates the temporal migration fractions has a stable cycle and a periodic behaviour. In subsection 3.2 we assume the temporal migration rates are given by a uniform distribution.

3.1 Temporal migration given by a Dirac measure

A Dirac measure is a measure δ_y defined on a set E such that

Let δ_m0 denote the Dirac measure centered on the fixed migration rate m ₀. Thus, we have

In this case, the criterion established in (3.11) is the same established by Earn et al. ⁽⁴⁾ that considered a metapopulation with any number of patches arbitrarily connected.

If we assume that the probability measure is concentrated in two periodic points, m ₀ and m ₁, we have

In this case, the quantifier Λ is the geometric average of (1 - m ₀λ_i) and (1 - m ₀λ_i). In fact, if the migration rates are distributed in p periodic points, m ₀, m ₁, ..., m _{p - 1}, the quantifier Λ is given by the following geometric average

3.2 Temporal migration given by a uniform distribution

Now we assume that the temporal migration rates are given by a uniform distribution. In our case, the probability density function with a uniform distribution on a set [a, b] ⊂[0, 1] is

Observe that . Besides that, we can write (3.13) as

The above integral can be solved analytically resulting

In the following, we perform numerical simulations to illustrate the stability regions of synchronized solutions and present the values of the quantifier Λ to different temporal migration rules.

4 NUMERICAL SIMULATIONS

We perform numerical simulation to illustrate the behavior of our network of patches connected by temporal dependent migration. In order to determine whether synchronization occurs wedefine the synchronization error, e _t, by

where . Synchronization is detected when et → 0.

We consider that the local dynamics is given by the following Ricker function

where x represents the patch density and r is the intrinsic growth rate (r > 0). The dynamics of a local habitat with the Ricker model is well-known and exhibits equilibrium, periodic and chaotic dynamics depending on the growth rate ⁽¹⁰⁾. For 0 < r < 2, the local dynamics become a state of equilibrium. For 2 < r < 2.526, the equilibrium point is unstable and a two-periodic trajectory takes its place. As r is increased, there appears a four-periodic trajectory and the two-periodic become unstable and we can observe period doubling bifurcations to chaos. In order to simulate chaotic within patch dynamics we assume r = 3.1, which implies that the isolated patch model has a chaotic trajectory with Lyapunov number L ≈ 1.327.

The configuration matrix Γ can be defined in different forms. Two well-known configurations are the nearest neighbor coupling and the global coupling ⁽⁴⁾. We illustrate our results considering the patches linked in a ring format with the two-nearest neighbors coupling, whose configuration matrix, Γ, is given by

In this case, the eigenvalues of Γ are given by λ₀ = 1 and , i = 1, 2, ..., d - 1. Of course, other local patch dynamics and configuration network topologies could be used, but our main concernment is to show the different behavior generated by temporal migration rates.

Figure 2 shows the bifurcation diagram of the synchronization error and the respective maximum transversal Lyapunov number versus the migration rate. In all cases, individuals migrate with the same rate at time t. We can observe that the non-synchronization region is characterized by weak interaction. Moreover, the increase in the number of patches decreases the synchronization region and increase the maximum transversal Lyapunov number. In fact, the subdominant eigenvalue of the matrix Γ tends to one as n → +∞ (λi → 1 as d → +∞, i = 1, 2, ..., d - 1), thus the quantifier Λ given in (3.15) also tends to one as d → +∞ (see ⁽¹²⁾). It means that the stability criterion will approach to the Lyapunov number if we increase the number of patches. Moreover, the synchronized attractors will be unstable for any value of the migration rate if the local dynamics of a single isolated patch is chaotic (L > 1).

Figure 2:
Synchronization error ((a), (b) and (c)) and respectably maximum transversal Lyapunov number ((d), (e) and (f)) vs m. Local dynamics is given by the Ricker function f(x) = x exp(r(1 - x)) with r = 3. 1. The patches are coupled with the two-nearest neighbors coupling. (a) 2 patches. (b) 5 patches. (c) 10 patches.

Figure 3 shows the metapopulation behavior with periodic migration. We consider 5 patches and three different scenarios. In all cases, individuals migrate according to a two periodic rule. In the first case, individuals migrate with a periodic rate given by 0. 1 and m (Fig. 3 (a)). In the second case, the migration rates are given by 0. 5 and m (Fig. 3 (b)), while in the third case it is given by 0. 9 and m (Fig. 3 (c)). We can observe that weak interactions between the patches decreases the region of synchronization, while intermediate and high migration rates have an oppositive effect.

Figure 3:
Synchronization error ((a), (b) and (c)) and respective maximum Lyapunov number ((d), (e) and (f)) vs m. Local dynamics is given by the Ricker function f(x) = x exp(r(1 - x)) with r = 3. 1. Five patches are coupled with the two-nearest neighbors. (a) m ₀ = m and m ₁ = 0. 1. (b) m ₀ = m and m ₁ = 0. 5. (c) m ₀ = m and m ₁ = 0. 9.

Table 1 shows different migration rules and the values of the quantifier Λ. We observe that, if the temporal migration rates are distributed around an average, the values of the quantifier Λ won't change its values significantly. In Table 1 the migrations rates are distributed around an average of m = 0. 3. In all cases, we observe that the values of Λ do not change its value significantly due to different migration rules when compared with the average of the migration rates.

Table 1:
Quantifier Λ for different temporal migration rates. Temporal migration rates are distributed around m = 0. 3. We can observe that the values of Λ do not change significantly.

	Fixed Point	Period 2	Period 4	Uniform	Uniform
m	0. 3	0. 2 and 0. 4	0. 15, 0. 25, 0. 35 and 0. 45	[0. 2, 0. 4]	[0. 1, 0. 5]
Λ	0. 7927	0. 78968	0. 78892	0. 7916	0. 7887

5 DISCUSSION

In this paper we develop a model of a network of equal patches linked by temporal dependent migration. The time evolution of the system involves two processes: local patch dynamics and migration between the patches. We then analyze the phenomenon of synchronization between the patches. We obtain an analytical criterion for the local asymptotic stability of synchronized trajectories based on the computation of the transversal Lyapunov numbers of attractors on the synchronous invariant manifold. The criterion is obtained via linearization process around the synchronized trajectories. The criterion is given by the product of two quantifiers: the separation rate of two nearby orbits in the single isolated patch measured by the Lyapunov number, L, and a quantifier that depends on the whole migration process, Λ (eq. 3.11). We then calculate the value of this quantifier to different migration rules. At first, we describe Λ assuming the migration rate with a periodic behavior (eq. 3.17), we then consider the migration rate uniform distributed on the interval [a, b] ⊆[0, 1] (eq. 3.20). The quantifier Λ in the case of migration rates with a periodic behavior involves the migration rates and the eigenvalues of the matrix that inform the network topology between the patches, while in the case of uniform distribution also involves the size of the interval.

Our observation based on theoretical results and on numerical simulations reveals the importance of analyzing a metapopulation model with temporal migration. We performed numerical simulation assuming each patch dynamics given by Ricker map with chaotic behavior. Then, we analyze the influence that the migration process has over synchronized dynamics. We observe that an increase in the number of patches decrease the stability regions (Fig. 2). Besides that, weak interactions between patches, decreases the size of stability regions, while intermediate and high migration rates have an opposite effect (Fig. 3). We observe that, if the temporal migration rates are distributed around an average, the values of quantifier Λ do not change significantly (Table 1).

• ¹
  [1] J. C. Allen, W. M. Schauffer & D. Rosko. Chaos reduces species extinction by amplifying local population noise. Nature, 364 (1993), 229-232.
• ²
  [2] B. Blasius, A. Huppert & L. Stone. Complex dynamics and phase synchronization in saptially extended ecological systems. Nature, 399 (1999), 353-359.
• ³
  [3] M. Doebeli. Dispersal and Dynamics. Theor. Pop. Biol. , 47 (1995), 82-106.
• ⁴
  [4] D. J., Earn, S. A. Levin & P. Rohani. Coherence and Conservation. Science, 290 (2000), 1360-1364.
• ⁵
  [5] J. P. Eckmann & D. Ruelle. Ergodic Theory of chaos and strange attractors. Rev. Modern Phys. , 57 (1985), 617-656.
• ⁶
  [6] L. Hansson. Dispersal and connectivity in metapopulations. Biol. J. Linn. Soc. , 42 (1991), 89-103.
• ⁷
  [7] M. Heino, V. Kaitala, E. Ranta & J. Lindström. Synchronous dynamics and rates of extinction in spatially structured populations. Proc. R. Soc. London B, 264 (1997), 481-486.
• ⁸
  [8] R. A. Ims & H. P. Andreassen. Spatial synchronization of vole population dynamics by predatory birds. Nature, 408 (2000), 194-196.
• ⁹
  [9] P. Lancaster & M. Tismenetsky. The Theory of Matrices, Academic Press, London (1985).
• ¹⁰
  [10] R. M. May. Biological populations with nonoverlapping generations: stable points, stable cycles, and chaos. Science, 186 (1974), 645-647.
• ¹¹
  [11] R. M. May& A. L. Lloyd. Synchronicity, chaos and population cycles: spatial coherence in an uncertain world. Trends Ecol. Evol. , 14 (1999), 417-418.
• ¹²
  [12] J. A. L. Silva, M. L. De Castro & D. A. R. Justo. Synchronism in a metapopulation model. Bull. Math. Biol. , 62 (2000), 337-349.
• ¹³
  [13] J. A. L. Silva& F. T. Giordani. Density-dependent migration and synchronism in metapopulations. Bull. Math. Biol. , 68 (2006), 451-465.

• V. Manica thanks CAPES for supporting his PhD at UFRGS.
• †
  The paper was presented at CMAC-CO 2013.

Publication Dates

• Publication in this collection
  Apr 2015

History

• Received
  28 Nov 2013
• Accepted
  04 Feb 2015