Mathematical Biosciences 250 (2014) 1–9

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Population distribution and synchronized dynamics in a metapopulation model in two geographic scales Vanderlei Manica ⇑, Jacques A.L. Silva Departmento de Matemática Pura e Aplicada-IM-UFRGS, Av. Bento Gonçalves, CEP 91509-900 Porto Alegre, RS, Brazil

a r t i c l e

i n f o

Article history: Received 26 March 2013 Received in revised form 17 September 2013 Accepted 4 February 2014 Available online 12 February 2014 Keywords: Metapopulation Distribution matrix Synchronization Lyapunov number Short-and long-distance dispersal

a b s t r a c t In this paper, a metapopulation model composed of patches distributed in two spatial scales is proposed in order to study the stability of synchronous dynamics. Clusters of patches connected by short-range dispersal are assumed to be formed. Long distance dispersal is responsible to link patches that are in different clusters. During each time step, we assume that there are three processes involved in the population dynamics: (a) the local dynamics, which consists of reproduction and survival; (b) short-range dispersal of individuals between the patches of each cluster; and (c) the movement between the clusters. First we present an analytic criterion for regional synchronization, where the clusters evolve with the same dynamics. We then discuss the possibility of a full synchronism, where all patches in the network follow the same time evolution. The existence of such a state is not always ensured, even considering that all patches have the same local dynamics. It depends on how the individuals are distributed among the local patches that compose a cluster after long-range dispersal takes place in the regional scale. An analytic criterion for the stability of synchronized trajectories in this case is obtained. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction A metapopulation is composed of populations that live in fragments called patches and are often connected via dispersal processes [20]. An interesting phenomenon related to the dispersal process is the synchronized dynamics where the population densities in all patches evolve with the same amplitude and phase [24]. Its importance lies in the fact that synchronized dynamics can make the whole metapopulation vulnerable to extinction. On the other hand, if the 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 [2,3]. The synchronization phenomenon is important from the ecological and epidemiological point of view [8]. In ecology, synchronization may have a deleterious effect on population persistence, because it may lead to the impossibility of a recolonization, and can be dangerous for species that need to be preserved. In epidemiology, synchronization can be beneficial in order to control and improve efforts against a disease. In Thailand [5], time series analyzes of dengue cases between 1984 and 1996 showed a spatial synchrony in the number of dengue cases between the cities indicating that the increase in number ⇑ Corresponding author. Tel.: +55 5198115293. E-mail address: [email protected] (V. Manica). http://dx.doi.org/10.1016/j.mbs.2014.02.002 0025-5564/Ó 2014 Elsevier Inc. All rights reserved.

of cases in a city may be reflected over the entire country. Another example is the seasonal dynamics of the influenza virus which epidemics occur annually with the highest activity occurring during winter months [38]. Systems of discrete equations are often used to model metapopulations [1,7,13,32,34]. A metapopulation model with patches linked by migration and subjected to external perturbations was considered in [1]. Through numerical simulations, it was shown that chaos can prevent global extinctions when coupling is weak. In [13] the patches were linked by considering dispersal process and distance and it was concluded that this asynchrony reduces the synchronization likelihood. In [7] was established a positive correlation between the degree of coherence of the oscillations in each patch and the risk of extinction of the metapopulation. They also 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 [34], concluding that this mechanism reduces the stability of the synchronous dynamics. These results are based on the calculation of the transversal Lyapunov number that depends on the local dynamics of each patch, and a parameter that depends on the whole migration process. Symmetric and asymmetric dispersion [6,12,16,18,21], density-dependent dispersion [29,33,34,37] studies were done in order to show the different trends of

2

V. Manica, J.A.L. Silva / Mathematical Biosciences 250 (2014) 1–9

metapopulation dynamics and its implications on the persistence and conservation of species. The dispersal process often involves two modes, short-distance and long-distance dispersal [30]. Short-distance dispersal refers to movement to other sites on a local scale usually by the individual’s own way of transport such as flying or swimming. The long-distance dispersal occurs in a regional scale and it takes place usually through passive transport on wind, water flow, birds or artificial transportation. Examples of long-range dispersal have been reported on insects [22,26,27], and in aquatic invertebrates [11,17]. Long-range dispersal events may be rare but they are crucial to population spread since it allows recolonization and counteracts genetic drift and inbreeding in fragmented environments [36]. In [14] the genetic structure of the medfly Ceratitis capitata in South Africa was analyzed and molecular approaches were employed to obtain estimates of the dispersal ability of the Mediterranean fruit fly at three spatial scales. The results in [14] suggested that the structure observed in South African medflies maybe the result of complex interactions at local scales (limited dispersal ability) and at broad scales (human-mediated and other forms of long-range dispersal). In this paper we propose a model of a network of local populations linked by a dispersal process that takes into account shortand long-distance movements. The environment is assumed to be fragmented in such way that clusters of patches connected by short-range dispersal are formed. We assume that these clusters are too far away to admit connection by short-range dispersal and thus, long-range dispersal is responsible for establishing connections between some patches that are distant apart. We then describe the synchronization phenomenon. The analysis is done by linearizing the equations of the model around the synchronized trajectories and further decomposition of perturbation vectors into components in the synchronized manifold and other components that are transversal to it, obtaining conditions to its local asymptotic stability. These conditions are obtained from the block decomposition of the Jacobian Matrix that presents a fundamental role in the stability analysis of the synchronous manifold [4,28]. In Section 2 we introduce the metapopulation model with patches distributed in two geographic scales. 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 considering different distributions of individuals in the patches that compose a cluster. Final comments and discussion are done in Section 5. 2. The mathematical model

of patches linked by short-range dispersal as in [1,7,13]. Assume that the processes of survival and reproduction which compose the local dynamics is described by a map f on ½0; 1Þ of class C 1 . In the absence of dispersal between patches, the time evolution of the population is given by

xtþ1 ¼ f ðxt Þ;

t ¼ 0; 1; 2; . . . ;

ð1Þ

where xt represents the number of individuals at time t. Important examples of f used in ecology are given in [25,35], and the single patch model (1) can display rich dynamical behavior including stable cycles, periodic-doubling cascades and chaos. We assume that a fraction m leave patch i and disperse to the neighboring patches. Thus, the density of individuals that leave patch i is given by mf ðxit Þ, where xit denote the population density in patch i at time t, for all i ¼ 1; . . . ; d; t ¼ 0; 1; . . .. Moreover, from the individuals that disperse from the neighboring patches k, a P fraction cik reach patch i. Clearly ckk ¼ 0, and dk¼1 cki ¼ 1, for all i ¼ 1; . . . ; d. We now can write the equations describing the dynamics of the isolated cluster as

xitþ1 ¼ ð1  mÞf ðxit Þ þ

d X

cik mf ðxkt Þ:

ð2Þ

k¼1

The first term in Eq. (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 of the neighboring patches. 2.2. A system of linked clusters The metapopulation model in two spatial scales consist of n equal clusters, as described in the previous subsection, labelled as 1; 2; . . . ; n, each one containing d equal patches. The total number of patches in the whole network is nd. Let xlit the number of l2 ld individuals in patch i of cluster l at time t. Let xlt ¼ ðxl1 t ; xt ; . . . ; xt Þ the d-dimensional vector representing the distribution of individuals among the patches of cluster l. In the absence of long-range dispersal the dynamics of cluster l is given by xltþ1 ¼ fðxlt Þ, where f : Rd ! Rd is given by

1 d X l1 c1k mf ðxlkt Þ C B ð1  mÞf ðxt Þ þ C B k¼1 C B C B d X C B l2 lk C B ð1  mÞf ðxt Þ þ c mf ðx Þ 2k t C B l fðxt Þ ¼ B C: k¼1 C B C B . .. C B C B C B d X A @ ld lk cdk mf ðxt Þ ð1  mÞf ðxt Þ þ 0

ð3Þ

k¼1

We consider a network of patches distributed in two spatial scales. In a local scale, nearby patches are connected by short-range dispersal forming clusters or conglomerates. We assume that these clusters are too far away from each other to be linked by short-range dispersal. Thus, in a regional scale, patches of different clusters are allowed to be connected by long-range dispersal (see Fig. 1). Fig. 1 is just a schematic representation but it resembles real networks topologies [15]. In [15] the network topology showing patch conglomerates was obtained with data on the grasshopper Stethophyma grossum distributed in a fragmented agricultural landscape. 2.1. The isolated cluster We start describing the basic habit unit of our network viewed in a broad scale, a cluster of patches linked by short-range dispersal. If there is no long range dispersal, the clusters are isolated. This will be a metapopulation in a more classical sense, a collection

fðxlt Þ

l2 ld ðyl1 t ; yt ; . . . ; yt Þ.

Let ¼ Thus after the local dynamics (within patch dynamics) and short-range dispersal the number of individuals of patch i of cluster l is ylit . We now describe the long distance movement between patches in different clusters. Let li be the fraction of individuals that leave patch i in any cluster in a long-range movement to establish in another cluster. Of course 0 6 li 6 1, i ¼ 1; . . . ; d. Thus, the number of individuals that leave patch i of cluster l at time t is li ylit . From these individuals only a fraction will move to cluster j. This process P is governed by a nonnegative n  n matrix C, satisfying nj¼1 cjl ¼ 1 and cll ¼ 0 for all l ¼ 1; . . . ; n. Thus, the number of organisms that leave patch i of cluster l and reach cluster j at time t is cjl li ylit . Let patch k of cluster l be the final destination of these migrants. Only part of them will settle in patch k with proportion wki . This process of distribution of migrants among the patches of the new cluster is governed by the d  d matrix W, with entries wki , with 0 6 wki 6 1 for all i; k ¼ 1; . . . ; d. Adding the contribution of all patches in clus-

V. Manica, J.A.L. Silva / Mathematical Biosciences 250 (2014) 1–9

Pd

ter l to the migration to patch k of cluster j we get cjl i¼1 wki li ylit (see Fig. 2). The distribution of migrants coming from all patches of cluster l among the patches of cluster j is given by the vector cjl WMfðxlt Þ, j ¼ 1; 2; . . . ; n, where M is d  d diagonal matrix defined by M ¼ diagðl1 ; l2 ; . . . ; ld Þ. Now, treating cluster l as a generic cluster and adding all the contributions to cluster j from migrants P of all clusters we obtain nl¼1 cjl WMfðxlt Þ as the distribution of migrants from all clusters among the patches of cluster j. Of course, the distribution of individuals that did not engage in this long distance dispersal process is given by ðI  MÞfðxjt Þ, where I is the d  d identity matrix. Therefore, we can write the dynamical system with n vectorial equations governing the time evolution of the network of nd patches as

xjtþ1 ¼ ½I  Mfðxjt Þ þ

n X cjl WMfðxlt Þ;

Further, we linearize the system of equations (4) around the synchronized trajectories obtaining a stability criterion based on the computation of the transversal Lyapunov numbers of attractors on the synchronous invariant manifold. Definition 1. The synchronized mode is given by S ¼ fXst ¼ ðxst ; xst ; . . . ; xst Þ 2 Rnd j xst ¼ ðx1t ; x2t ; . . . ; xdt Þ 2 Rd ; t ¼ 0; 1; 2; . . .g. The Definition 1 means that synchronization is achieved if the density of all clusters is the same at time t. Substitution of xjt ¼ xst ; j ¼ 1; . . . ; n, in Eq. (4) leads us to the existence of such synP chronized solution provided nl¼1 cjl ¼ 1; j ¼ 1; . . . ; n. Furthermore, the dynamics of each cluster in the synchronized mode satisfies

xstþ1 ¼ ðI  M þ WMÞfðxst Þ: j ¼ 1; . . . ; n:

ð4Þ

l¼1

We call W the distribution matrix. It is of central importance in our model. The matrix W tell us how the migrants in the departure cluster are distributed as well as how they distribute in the destination cluster after the long-range dispersal. Notice that the sum of any column of W must be either zero or one. The case where the ith column of W is zero means that the ith patch does not contribute to the long-range dispersal. If there is long distance dispersal from patch i, the ith column of W must add up to one. For example, if the ith row of W is zero then we are in the case where patch i does not receive contribution from the long-range dispersal process. All entries in the ith row of W are equal to one represents the case that patch i receives all the migrants coming from long distances. A complete uniform distribution is obtained in the case all the entries of W are equal 1d. In this case all patches in the departure cluster have equal contribution to the long distance dispersal process. Moreover, all the migrants distribute uniformly in all patches of the destination cluster. Notice that these considerations about how migrants distribute in the departure and destination clusters hold independently on which are these departure and destination clusters. More realistic models should include such dependence thus we would have to consider a family of distribution matrices W jl ; j – l ¼ 1; 2; . . . ; n. Of course, in such models the basic habitat units (the clusters of patches connected by short-range dispersal) would not be all equal and therefore, the network would lack some essential symmetries to ensure the existence of synchronism, which is the issue discussed next. 3. Synchronization and transversal stability In the following, we define the synchronized mode for our metapopulation model with patches in two geographic scales.

Fig. 1. Schematic representation of the network topology with 5 clusters and 4 patches.

3

ð5Þ

The solution set of (5) represents our synchronized attractor and it depends on the ‘‘local’’ cluster dynamics given by f and the dispersal process in both scales. Observe that if individuals migrate preferentially to patch of the same index in the neighboring cluster, the matrix W is exactly the identity matrix, and the synchronized attractor is given by the local cluster dynamics described by f. We are interested in studying the local asymptotic stability of the synchronized mode, that is, whether orbits that initiate close to the synchronized mode will be attracted to it. In order to achieve this goal, we assume that C is diagonalizable and use the following theorem that allows us to decompose the Jacobian matrix of system (4) evaluated at the synchronized mode into n blocks of dimension d  d. Theorem 1. Assume that C is diagonalizable with eigenvalues kj , j ¼ 0; 1; . . . ; n  1. Let Dt 2 Rnd be a perturbation from the synchronized mode. Then, there is a change of variables of the type Y t ¼ P Dt , where P is an invertible nd  nd matrix, such that, the linearization of system (4) around the synchronous mode take the form n1

Y tþ1 ¼ a½I  M þ kj WMDfðxst ÞY t ;

ð6Þ

j¼0

where Dfðxst Þ is the Jacobian matrix of the isolated cluster system (3). Proof. Standard linearization around Xst yields

Dtþ1 ¼ JðXst ÞDt ;

ð7Þ

nd

where Dt 2 R is the perturbation of the synchronized mode, and JðXst Þ is the nd  nd Jacobian matrix of system (4) evaluated at Xst .

Fig. 2. Long-range dispersal process. A fraction li of individuals leave patch i in the cluster l. From this individuals a fraction cjl will move to cluster j and only part of them will settle in patch k with proportion wki ; i; k ¼ 1; 2; . . . ; d; j; l ¼ 1; 2; . . . ; n. Thus, the fraction of individuals that leave patch i in cluster l and reach patch k in cluster j is cjl wki li .

4

V. Manica, J.A.L. Silva / Mathematical Biosciences 250 (2014) 1–9

Using the Kronecker product1 allow us to write JðXst Þ as

JðXst Þ ¼ I  Dfðxst Þ  I  MDfðxst Þ þ C  WMDfðxst Þ:

ð8Þ

We now proceed with the decomposition of the system (7). Since C is diagonalizable, there exists a single invertible matrix Q such that K ¼ QCQ 1 , where K ¼ diagðk0 ; k1 ; . . . ; kn1 Þ, and kj are the eigenvalues of C; j ¼ 0; 1; . . . ; n  1. Making the change of variables Y t ¼ ðQ  IÞDt , and taking (7) and (8) into account we may write

rion for the transversal stability of the synchronous mode (K < 1). The value h is associated with the synchronized attractor, and chaotic trajectories are observed if h > 1, while periodic trajectories at h < 1. We summarize the results in the following theorem. Theorem 2. Assume that C is diagonalizable. Let q be the natural probability measure of (15). Assume that þ ln k½I  M þ kj WMDfðxÞk 2 L1 ðqÞ; j ¼ 1; . . . ; n  1. Then, the criterion for the local stability of attractors of the synchronized mode generated by (5) is K < 1.

Y tþ1 ¼ ðQ  IÞðI  Dfðxst Þ  I  MDfðxst Þ þ C  WMDfðxst ÞÞDt ¼ ðQ  Dfðxst Þ  Q  MDfðxst Þ þ QC  WMDfðxst ÞÞDt :

ð9Þ

But Dt ¼ ðQ  IÞ1 Y t ¼ ðQ 1  IÞY t , thus

Dfðxst Þ

Y tþ1 ¼ ðI 

I

MDfðxst Þ



þK 

WMDfðxst ÞÞY t :

ð10Þ



Since K ¼ diagðk0 ; k1 ; . . . ; kn1 Þ, we have the block decomposition matrix, that is,

3.1. Full synchronization A particular case of the synchronized attractor generated by (5) is the fully synchronized mode, that is, all patches in the network evolving with the same density. It is a more restricted manifold and more symmetry must be imposed. We start with a definition of this state.

n1

Y tþ1 ¼ a½I  M þ kj WMDfðxst ÞY t :

ð11Þ

Definition 2. The fully synchronized mode is given by S ¼ fXst ¼ ðxst ; xst ; . . . ; xst Þ 2 Rnd j xst ¼ ðxst ; xst ; . . . ; xst Þ 2 Rd ; t ¼ 0; 1; 2; . . .g.

The importance of this theorem lies in the fact that the local stability of the synchronized mode of system (4) can be studied by analyzing the spectrum of the individual blocks in equation (11). It is important to notice that the connectivity matrix C is doubly stochastic (all rows and columns add up to one). A simple application of the Gershgorin Theorem [19] leads to the fact that k0 ¼ 1 is the dominant eigenvalue of C. The block corresponding to such eigenvalue corresponds to the variational matrix of the system Eq. (5) that generates the synchronized attractor, while the other n  1 blocks correspond to perturbations to the transversal directions of the synchronous invariant manifold and govern its local asymptotic stability. In order to describe the synchronized attractor behavior and its local asymptotic stability, we define the maximum parallel Lyapunov number, h, by

Substitution of xjt ¼ xst ; j ¼ 1; . . . ; n, in Eq. (4) leads us to the existence of fully synchronized solution provided Pn Pd l¼1 c jl ¼ 1; j ¼ 1; . . . ; n, and k¼1 cik ¼ 1; i ¼ 1; . . . ; d. Notice that (3) and Definition 2 yield fðxst Þ ¼ f ðxst Þ~ 1, where ~ 1 ¼ ½1 1 . . . 1T 2 Rd . Thus (14) and (15) lead to



j¼0

hðxs0 Þ ¼ lim kP0;s1  . . .  P0;1 P0;0 k1=s ; s!1

ð12Þ

where P0;t ¼ ½I  M þ WMDfðxst Þ, and the maximum transversal Lyapunov number, K, by

Kðxs0 Þ ¼ max

j¼1;...;n1



 lim kPj;s1  . . .  Pj;1 Pj;0 k1=s ;

s!1

ð13Þ

where Pj;t ¼ ½I  M þ kj WMDfðxst Þ; j ¼ 1; . . . ; n  1, for all initial points xs0 . Denoting the map that generates the synchronized attractor by g, i.e.,

gðxÞ ¼ ðI  M þ WMÞfðxÞ:

ð14Þ

The synchronized attractor dynamics can be written as

xstþ1 ¼ gðxst Þ:

ð16Þ

t¼0

n

1 ¼ a~ 1 and notice that associated to it. Now suppose A~ T T~ ~ ~ ~ On the other hand, we have ðA1Þ 1 ¼ a1 1 ¼ ad. T ðA~ 1Þ ~ 1¼~ 1T AT ~ 1¼~ 1T ~ 1 ¼ d and thus, a ¼ 1. Therefore, if ~ 1 is an eigenvector of A, the fully synchronized mode exists and all patches will evolve in time according to xstþ1 ¼ f ðxst Þ, which is equivalent to Eq. (1), the single patch model equation. In other words the whole network is fully synchronized and its dynamics is identical to the dynamics of a single isolated patch. We are interested in studying the local asymptotic stability of fully synchronized attractors. In order to do it, we can use Theorem 1 to get the block decomposition and proceed in an similar way to the previous subsection to calculate the maximum transversal Lyapunov number. Observe that the variational matrix Dfðxst Þ of (3), applied in xst ¼ ðxst ; xst ; . . . ; xst Þ 2 Rd , has its entries given by



ð1  mÞf 0 ðxst Þ; if k ¼ i;

cki mf 0 ðxst Þ;

if k – i;

Dfðxst Þ ¼ ðI  mBÞf 0 ðxst Þ;

where dx is the Dirac point mass at x. Assuming the integrability of þ ln k½I  M þ kj WMDfðxst Þk with the ergodicity of q, we can apply the Ergodic Theorem of Oseledec [9] to guarantee the existence and uniqueness q-almost everywhere of the limits defining h in (12) and K in (13). The maximum transversal Lyapunov number K of an attractor on the synchronous invariant manifold gives a criteBe A ¼ ½aij m 2 Rmm and B ¼ ½bij i;j¼1 2 Rnn , the Kronecker product is defined by i;j¼1 A  B ¼ ½aij Bnij¼1 2 Rmnmn . 1

We observe that the vector xst given by the above equation

which can be written as

Let q be the natural probability measure for g,

s!1

ð17Þ

has all identical components if and only if ~ 1 is an eigenvector of the matrix A ¼ I  M þ WM. Notice that all the columns of A add up to one, thus k0 ¼ 1 is an eigenvalue of AT and ~ 1 is an eigenvector

aki ¼ ð15Þ

s1 X q ¼ lim dxt ;

xstþ1 ¼ f ðxst ÞðI  M þ WMÞ~ 1:

where B ¼ I  C, and C is the coupling matrix between the patches in the local scale. Thus, P j;t ¼ ðI  M þ kj WMÞðI  mBÞf 0 ðxst Þ; j ¼ 1; . . . ; n  1. Observing that

kPj;s1  . . .  Pj;1 Pj;0 k ¼

s1 Y

!

jf

0

ðxst Þj

kððI  M þ kj WMÞðI  mBÞÞs k;

t¼0

ð18Þ the maximum transversal Lyapunov number can be written as

Kðxs0 Þ ¼ max

j¼1;...;n1



 lim kP j;s1  . . .  Pj;1 Pj;0 k1=s ¼ Lðxs0 ÞK;

s!1

ð19Þ

5

V. Manica, J.A.L. Silva / Mathematical Biosciences 250 (2014) 1–9

Qs1

1=s

where Lðxs0 Þ ¼ lims!1 ð t¼0 jf 0 ðxst ÞjÞ is the Lyapunov number of f starting in xs0 , and K is the maximum of the spectral radius of ðI  M þ kj WMÞðI  mBÞ, j ¼ 1; . . . ; n  1. In this case, the natural probability measure can be defined for þ the local map f. Assuming the integrability of ln jf 0 ðxÞj, we can apply the Ergodic Theorem of Birkhoff to guarantee the existence and uniqueness q-almost every xs0 of the limit defining L. Observe that L depends on the local patch dynamics, while K depends on the migration process in both scales. Moreover, full chaotic synchronization is observed when L > 1 and K < 1L . We summarize the results in the following theorem. Theorem 3. Assume that C is diagonalizable. Let q be natural probability measure generated by Eq. (1). Assume that þ ln jf 0 ðxÞj 2 L1 ðqÞ and ½1 1 . . . 1T is an eigenvector of A ¼ I  M þ WM. Then, the criterion for the local stability of attractors of the fully synchronized mode generated by (1) is K ¼ LK < 1. 4. Applications and numerical simulations We perform numerical simulation to illustrate the behavior of our network of patches distributed in two geographical scales. The ensemble trajectories are calculated using (4), while the dynamics of the synchronized modes are obtained using (5). In order to determine whether synchronization occurs in the metapopulation we define the synchronization error, et , given by

et ¼

x1;i t .

f ðxÞ ¼ xerð1xÞ ;

ð21Þ

where x represents the patch density and r is the intrinsic growth rate (r > 0). The dynamics of a local patch with the Ricker model exhibits equilibrium, periodic or chaotic dynamics depending on the intrinsic growth rate [23]. In order to simulate chaotic within patch dynamics we assume r ¼ 3:1, which implies that the isolated patch model have typical chaotic trajectories with Lyapunov number L  1:327. The configuration matrix C can be defined in different forms. Two well-known configurations are the nearest neighbor coupling and the global coupling [7]. We illustrate our results considering the clusters linked in a ring format with the two nearest neighbors clusters (two-nearest neighbors coupling), whose configuration matrix, C 1 , is given by

0

1=2

0

...

0

1=2

1

B 1=2 0 1=2 0 ... 0 C C B B . . . .. C C B . . . B 0 1=2 . . . . C C; C1 ¼ B C B . .. .. B .. . . 0 1=2 0 C C B C B @ 0 ... 0 1=2 0 1=2 A 1=2

0

...

0

1=ðn  1Þ

B B 1=ðn  1Þ B C2 ¼ B .. B @ .

0 .. .

... .. . .. .

1=ðn  1Þ

...

1=ðn  1Þ

1 1=ðn  1Þ C .. C . C C: C 1=ðn  1Þ A

ð23Þ

0

Both configuration matrices are circulant and, thus, diagonalizable. Moreover, the eigenvalues of C 1 are given by k0 ¼ 1 and kj ¼ cosð2npjÞ, while the eigenvalues of the matrix C 2 are given by 1 k0 ¼ 1 and kj ¼  n1 , for all j ¼ 1; . . . ; n  1. Analogously, we can define the configuration matrix C that represents the connectivity between the patches that compose a cluster. In our simulations we consider clusters composed by two patches and two different ways that individuals can distribute among these patches after the long-range dispersal takes place. In the first case individuals of all patches in neighbor clusters are equally connected, the uniform distribution. In a second example we assume that are different preferences in the way individuals distribute, the triangular distribution. In Sections 4.1 and 4.2 we describe the model considering the case of uniform distribution and triangular distribution, respectively. In both cases, there is no short-distance dispersal (m ¼ 0). Section 4.3 presents a more general case where we assume short-and long-distance dispersal.

ð20Þ

where ¼ Synchronization is detected when et ! 0. Moreover, the maximum parallel Lyapunov number (h) and the maximum transversal Lyapunov number (K) are calculated numerically using an algorithm described in [31] that calculates the average separation of trajectories in orthogonal directions. If K < 1, it means that the synchronized state is stable, consequently perturbations near to this state tend to zero and the synchronization state is reached. We consider that the local dynamics of each patch is given by the following Ricker function

0

0

4.1. Uniform distribution with two patches in each cluster

d X n 1 X jxj;i  xjþ1;i j t nd i¼1 j¼1 t

xnþ1;i t

and all clusters equally coupled (globally connected), whose matrix, C 2 , is given by

0

1=2

0

ð22Þ

We start describing the case where individuals that reach a new cluster choose between patch 1 and 2 with the same probability. Thus, the distribution matrix is given by

 W¼

0:5 0:5 0:5 0:5

 :

ð24Þ

We call this case as uniform distribution and we can write the system of equations (5) that generates the synchronized attractors by

x1tþ1

! ¼

x2tþ1

y1t y2t

!

!





!

0:5 0:5 l1 y1t l1 y1t þ ;  2 0:5 0:5 l2 yt l2 y2t



ð25Þ

where li is the fraction of individuals that leave patch i and migrate i to neighboring clusters, and yit is given by yit ¼ xit erð1xt Þ ; i ¼ 1; 2. Thus, the above equations can be written as

x1tþ1 x2tþ1

! ¼

1

1

2

2

2

1

x1t erð1xt Þ  0:5l1 x1t erð1xt Þ þ 0:5l2 x2t erð1xt Þ x2t erð1xt Þ  0:5l2 x2t erð1xt Þ þ 0:5l1 x1t erð1xt Þ

! :

ð26Þ

Observe that Eq. (26) synthesize the behavior of two populations coupled by migration [6]. In the case of l :¼ l1 ¼ l2 , the equations (26) reduce to the case of two populations with symmetric dispersal [12]. Fig. 3 shows the synchronized attractors and its behavior for the case of l :¼ l1 ¼ l2 . We can observe that the synchronized attractors are mainly characterized by three different behaviors. For small values of l (l < 0:09), they are chaotic. For intermediate values of l ð0:09 6 l 6 0:25Þ, they are periodic and out of phase (x1tþ1 ¼ x2t ), and for values of l bigger than 0.25, they are over the phase diagonal (x1t ¼ x2t ). The stability of the synchronized attractor is given calculating the maximum transversal Lyapunov number. The term Pj;t in the limit (13) is given by

6

V. Manica, J.A.L. Silva / Mathematical Biosciences 250 (2014) 1–9

(a)

We note the perfect agreement between the criterion and the synchronization error. Moreover, the regions where the synchronized attractors are periodic (0:09 6 l 6 0; 25), the system is stable and the K values do not vary with the increase of the number of clusters or changes in the network topology. For l > 0:25, the attractors are fully synchronized and the stability region depends on the number of clusters and changes in the network topology. While for l < 0:09, the synchronized attractors are unstable and the clusters do not synchronize. Observe that l :¼ l1 ¼ l2 , we have M ¼ lI and ~ 1 is an eigenvector of A ¼ ð1  lÞI þ lW, since the lines of W add up to 1. For this case, as shown in Fig. 6, an increase in the number of clusters has a negative effect in the stability of the synchronized mode using the two-nearest neighbors topology.

(b)

2

2

xt

x2t

3

1

3

1

1

0 0

0.5 μ

0 0

1

(c)

0.5 μ

1

(d)

3

2

xt

h

2

2 1

1 1

2

0 0

3

1 xt

0.5 μ

1

Fig. 3. Behavior of the synchronized attractor for two patches and uniform distribution. Local patch dynamics is given by the Ricker map with r ¼ 3:1; m ¼ 0 and l :¼ l1 ¼ l2 . (a) x1t vs l. (b) x2t vs l. (c) x2t vs x1t for l ¼ 0:05 (black.), l ¼ 0:2 (blue+), and l ¼ 0:5 (red.). (d) Maximum parallel Lyapunov number vs l. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

P j;t ¼

ð1  l1 þ 0:5kj l1 Þf 0 ðx1t Þ

0:5kj l1 f 0 ðx2t Þ

0:5kj l1 f 0 ðx1t Þ

ð1  l2 þ 0:5kj l2 Þf 0 ðx2t Þ

! ;

where f 0 is the derivative of the Ricker map (21), x1t and x2t belong to the synchronized attractor given by the equations (26). Thus, the maximum transversal Lyapunov number (K) is calculated numerically considering the average separation of orbits in two orthogonal directions. Figs. 4 and 5 show the maximum transversal Lyapunov number and the synchronization error for two different network topologies (two-nearest neighbors coupling and globally connected network).

2

1 0

0

K

K

(a)

0.5

μ

(b)

1 0

1

0

0.5

μ

W¼



1

Fig. 4. Maximum transversal Lyapunov number vs l for two patches and uniform distribution. The local patch dynamics is given by the Ricker map with r ¼ 3:1. (a) 5 clusters. (b) 10 clusters. The thickness line represents the globally connected network, while the other one represents the two-nearest neighbors coupling. In the region where K < 1, the clusters have a synchronized trajectory.

x2tþ1

et

t

e

! ¼

1

2

2

2

x1t erð1xt Þ þ 0:5lx2t erð1xt Þ

! ð29Þ

:

x2t erð1xt Þ  0:5lx2t erð1xt Þ

4.3. More general cases In the previous sections, we did not consider the dispersal effect between the patches that compose a cluster. This local dispersal changes the behavior of the synchronized attractors as shown in

1

1

LΛ

2

0.5 μ

ð28Þ

(a)

2

(b)

(b)

2

0 0

0 0:5

 :

Observe that the above equations synthesize two populations coupled by dispersal where just one of them migrates [10]. We note that changing the distribution form between patches generate synchronized attractors totally different when compared with the case of uniform distribution, as shown in Fig. 7. Synchronized attractors corresponding to full synchronization are not observed and high values of l induce the synchronized attractors to a stabilizing effect, that is, it goes from chaos to periodic trajectories due to the dispersal process. Moreover, periodic synchronized attractors are stable, while the stability of chaotic synchronized attractors depends on factors like number of clusters and network topology (Figs. 8 and 9). In this case of triangular distribution, ~ 1 is not an eigenvector of A ¼ I  M þ WM since the lines of W do not add up to 1, consequently the phase diagonal is not an attractor and there is no full synchronization.

2

(a)

1 0:5

We call this case as triangular distribution and we can write the system of equations (5) that generates the synchronized attractors by

x1tþ1 ð27Þ

2

Now we assume that individuals in patch 1 move to patch number 1 in the neighboring cluster while individuals in patch 2 chose between patch 1 and 2 in the neighboring cluster with the same probability. The distribution matrix is given by

1

0 0

1

0 0 0.5 μ

1

Fig. 5. Synchronization error vs l for uniform distribution considering 5 clusters and 2 patches within each cluster. (a) two-nearest neighbors coupling. (b) globally connected network.

LΛ

0 0

4.2. Triangular distribution with two patches in each cluster

0.5 μ

1

1

0 0

0.5 μ

1

Fig. 6. Maximum transversal Lyapunov number for full synchronization vs l for two patches and uniform distribution. The local patch dynamics is given by the Ricker map with r ¼ 3:1. (a) 5 clusters. (b) 10 clusters. The thickness and black line represents the globally connected network, while the other one represents the two-nearest neighbors coupling.

V. Manica, J.A.L. Silva / Mathematical Biosciences 250 (2014) 1–9

(a)

(b) 4 x2t

xt

1

4 2 0 0

7

0.5 μ

2 0 0

1

(c)

0.5 μ

1

(d)

3

2

xt

h

2

2 1

1 0 0

2 1 xt

0 0

4

0.5 μ

1

Fig. 7. Behavior of the synchronized attractor for two patches and triangular distribution. Local patch dynamics is given by the Ricker map with r ¼ 3:1; m ¼ 0 and l :¼ l1 ¼ l2 . (a) x1t vs l. (b) x2t vs l. (c) x2t vs x1t for l ¼ 0:2 (black.), l ¼ 0:5 (red.), and l ¼ 0:8 (blue+). (d) Maximum parallel Lyapunov number vs l. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

(a)

2

1 0 0

K

K

2

0.5 μ

(b)

1 0 0

1

Fig. 10. Maximum parallel Lyapunov number ((a) and (c)), and correspondent maximum transversal Lyapunov number ((b) and (d)) vs m and l, l :¼ l1 ¼ l2 . (a) Uniform distribution. (c) Triangular distribution. The system has 5 clusters with 2 patches in each cluster. Local patch dynamics is given by the Ricker map with r ¼ 3:1 and two-nearest neighbors coupling.

0.5 μ

1

Fig. 8. Maximum transversal Lyapunov number vs l for two patches and triangular distribution. The local patch dynamics is given by the Ricker map with r ¼ 3:1. (a) 5 clusters, (b) 10 clusters. The thickness and black line represents the globally connected network, while the other one represents the two-nearest neighbors coupling. In the region where K < 1, the clusters have a synchronized trajectory.

(a)

(b)

2

2

0 0

t

e

e

t

Fig. 11. Same as Fig. 10, but 10 patches within each cluster.

1

0.5 μ

1

1

0 0

0.5 μ

1

Fig. 9. Synchronization error vs l for triangular distribution considering 5 clusters and 2 patches within each cluster. (a) two-nearest neighbors coupling, (b) globally connected network.

reach a steady or periodic state (Fig. 12). Moreover, while an increase in the number of clusters favors the instability of synchronized dynamics with chaotic behavior (see Figs. 4 and 8), an increase in the number of patches that compose a cluster leads to an increase in the region where the synchronized attractors are periodic (Figs. 10–12). 5. Discussion

Fig. 10. For the case of uniform distribution, we can identify two regions where the synchronized attractors are periodic (Fig. 10(a)). Moreover, in those regions, the synchronized attractors are locally asymptotically stable (Fig. 10(b)). For the triangular distribution, the synchronized attractors are located in different regions (Fig. 10(c)). We observe that high dispersal fractions between the clusters lead the model to a stable periodic synchronized attractor, except when the dispersal fractions between the patches also are high (Fig. 10(d)). An increase to ten patches in each cluster changed significantly the values of h and K in comparison with the two-patch case (see Fig. 11). We notice that the maximum parallel Lyapunov number is closer to the value 1, especially for the triangular migration case. It means that the separation of the trajectories get smaller, and

In this paper we develop a model of a network of equal patches linked by dispersal in two spatial scales. In a local scale, nearby patches are connected by short-range dispersal forming clusters. We assume that these clusters are too far away from each other to be linked by short-range dispersal. Thus, in a regional scale, we assume that clusters of patches are linked by long-distance dispersal (Fig. 1). The time evolution of the system involves three processes: local patch dynamics, dispersal of individuals between patches that compose a cluster, and dispersal between the clusters. The first two processes can be described by well-known models in the literature [1,7,37], while the model in the regional scale incorporate the dispersal process between the clusters of patches. This process considers a distribution matrix, W, that tell us how the

8

V. Manica, J.A.L. Silva / Mathematical Biosciences 250 (2014) 1–9

(a)

(b)

3

4

xt

x1t

1

2 2

1 0 0

0.5 μ

0 0

1

(c)

1

(d)

3

8 6 x1t

xt

1

2 1 0 0

0.5 μ

4 2

0.5 μ

1

0 0

0.5 μ

1

Fig. 12. Synchronized attractors of patch indexed as 1 vs l, considering m ¼ 0:1 and two-nearest neighbors coupling. (a) 2 patches and (c) 10 patches with uniform distribution. (b) 2 patches and (d) 10 patches with triangular distribution.

migrants in the departure cluster are distributed as well as how they distribute in the destination cluster after the long-distance dispersal. For example, individuals can migrate preferentially to certain patches or they can distribute equally in all patches that compose the neighboring cluster. We then analyze the phenomenon of synchronization between the clusters of 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 decomposition of the Jacobian matrix evaluated at the synchronized mode in blocks (Theorem 1) allows us to describe the synchronized trajectories behavior and its local stability. The key point is that one of the blocks corresponds to the variational matrix of the equations that generate the synchronized attractors. From this block, we calculate the maximum parallel Lyapunov number h that provides information about the synchronized attractors, i.e., periodic (h < 1) or chaotic dynamics (h > 1). The other blocks correspond to transversal directions of the synchronous invariant manifold and it allows us to calculate the maximum transversal Lyapunov number (K) that inform whether orbits that initiate close to the synchronized attractor are attracted (K < 1) or repelled (K > 1) to it (Theorem 2). A case of full synchronization where all patches evolve with the same density is also investigated. In this case, the criterion is given by the product of two quantifiers: the separation rate of two nearby orbits in the single patch measured by the Lyapunov number (L), and a term that depends on the whole dispersal process in both scales (Theorem 3). Our observation based on theoretical results and on numerical simulations reveals the importance of analyzing a metapopulation model composed by clusters of patches. In a situation where each cluster has two patches with chaotic local dynamics, the numerical simulations for two different distributions give very versatile results. The case where individuals that reach the neighboring cluster distribute equally in all patches, uniform distribution, generates different synchronized attractors (Fig. 3). We observe that weak interaction generate unstable synchronized attractors with chaotic behavior, intermediate dispersal fractions induce the model to stable synchronized attractors with periodic behavior, while high dispersal fractions induce to a fully synchronized attractor and its stability depends on factors like dispersal rate, number of

clusters and network topology (Fig. 4). When we change the distribution, the synchronized attractors are totally different. If we assume that individuals in patch number 1 move only to patch number 1 in the neighboring cluster while individuals in patch number 2 choose between patch 1 and 2 in the neighboring clusters with the same probability, triangular distribution, fully synchronized trajectories are not observed and as the dispersal rate is increased, the synchronized trajectories has a stabilizing effect, that is, it goes from chaos to periodic trajectories due to the dispersal process (Fig. 7). Moreover, periodic synchronized attractors are stable, while the stability of synchronized attractors with a chaotic behavior, as in the uniform distribution case, depends on factors like dispersal fraction, number of clusters and network topology (Fig. 8). The fact that dispersal has a stabilizing effect on the dynamics of populations is shown in studies in the literature [6,12,21], and it depends on the dispersal rule [37]. Moreover, stable synchronous dynamics (all patches evolving with the same density) depend on factors like environmental disturbance [1], distance between populations [13], density-dependent dispersal [33]. In our model in two spatial scales, we do not consider factors like environmental disturbance, distance between populations, density-dependent dispersal. On the other hand, our metapopulation model considers two dispersal processes for each generation and a preference for individuals to settle in certain patches. Moreover, our model can be used for testing different hypotheses and different dispersal rules in order to describe metapopulations patterns in ecology and epidemiology. Acknowledgments VM thanks CAPES for supporting his PhD at UFRGS, Brazil, A⁄STAR Research Attachment Programme (ARAP) and Dr. Fu Xiuju for supporting and supervising during his internship at IHPC, Singapore. Also thanks to Rogerio Manica and Erickson Tjoa for helpful discussions. References [1] J.C. Allen, W.M. Schauffer, D. Rosko, Chaos reduces species extinction by amplifying local population noise, Nature 364 (1993) 229. [2] B. Blasius, A. Huppert, L. Stone, Complex dynamics and phase synchronization in spatially extended ecological systems, Nature 399 (1999) 353. [3] B. Cazelles, S. Bottani, L. Stone, Unexpected coherence and conservation, Proc. R. Soc. London B 268 (2001) 2595. [4] P. Checco, M. Biey, L. Kocarev, Synchronization in random networks with given expected degree sequences, Chaos Solitons Fract. 35 (2008) 562. [5] D.A.T. Cummings, R.A. Irizarry, N.E. Huang, T.P. Endy, A. Nisalak, K. Ungchusak, D.S. Burke, Travelling waves in the occurrence of dengue haemorrhagic fever in Thailand, Nature 427 (2004) 344. [6] M. Doebeli, Dispersal and Dynamics, Theor. Popul. Biol. 47 (1995) 82. [7] D.J.D. Earn, S.A. Levin, P. Rohani, Coherence and conservation, Science 290 (2000) 1360. [8] D.J.D. Earn, R. Pejman, B.T. Grenfell, Persistence, chaos and synchrony in ecology and epidemiology, Proc. R. Soc. London B 265 (1998) 7–10. [9] J.P. Eckmann, D. Ruelle, Ergodic theory of chaos and strange attractors, Rev. Mod. Phys. 57 (1985) 617. [10] J.E. Franke, A. Yakubu, Periodic dynamical systems in unidirectional metapopulation models, J. Differ. Equ. Appl. 11 (2005) 687. [11] A.J. Green, J. Figuerola, Recent advances in the study of long-distance dispersal of aquatic invertebrates via birds, Divers. Distrib. 11 (2005) 149. [12] A. Hastings, Complex interactions between dispersal and dynamics: lessons from coupled logistic equations, Ecology 74 (1993) 1362. [13] 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. [14] M. Karsten, B.J. van Vuuren, A. Barnaud, J.S. Terblanche, Population genetics of Ceratitis capitata in South Africa: implications for dispersal and pest management, PLoS ONE 8 (1) (2013) e54281, http://dx.doi.org/10.1371/ journal.pone.0054281. [15] D. Keller, R. Holderegger, J. Marten, M.J. van Strien, Spatial scale affects landscape genetic analysis of a wetland grasshopper, Mol. Ecol. 22 (9) (2013) 2467.

V. Manica, J.A.L. Silva / Mathematical Biosciences 250 (2014) 1–9 [16] B.E. Kendall, Spatial structure, environmental, heterogeneity, and population dynamics: analysis of the coupled logistic map, Theor. Popul. Biol. 54 (1998) 11. [17] B.P. Kinlan, S.D. Gaines, S.E. Lester, Propagule dispersal and the scales of marine community process, Divers. Distrib. 11 (2005) 139. [18] D. Kleinhans, P.R. Jonsson, On the impact of dispersal asymmetry on metapopulation persistence, J. Theor. Biol. 290 (2011) 37. [19] P. Lancaster, M. Tismenetsky, The Theory of Matrices, Academic Press, London, 1985. [20] R. Levins, Some demographic and genetic consequences of environmental heterogeneity for biological control, Bull. Entomol. Soc. Am. 15 (1969) 237. [21] A.L. Lloyd, The coupled logistic map: a single model for the effects of spatial heterogeneity on population dynamics, J. Theor. Biol. 173 (1995) 217. [22] H.D. Loxdale, J. Hardie, S. Halbert, R. Foottit, N.A.C. Kidd, C.I. Carter, The relative importance of short- and long-range movement of flying aphids, Biol. Rev. 68 (1993) 291. [23] R.M. May, Biological populations with nonoverlapping generations: stable points, stable cycles, and chaos, Science 186 (1974) 645. [24] R.M. May, A.L. Lloyd, Synchronicity, chaos and population cycles: spatial coherence in an uncertain world, Trends Ecol. Evol. 14 (1999) 417. [25] R.M. May, A.L. Lloyd, Simple mathematical models with very complicated dynamics, Nature 261 (1976) 1. [26] A. Meats, J.E. Edgerton, Short- and long-range dispersal of the Queensland fruit fly, Bactrocera tryoni and its relevance to invasive potential, sterile insect technique and surveillance trapping, Aust. J. Exp. Agric. 48 (9) (2008) 1237.

9

[27] A. Meats, C.J. Smallridge, Short- and long-range dispersal of medfly, Ceratitis capitata (Dipt., Tephritidae), and its invasive potential, J. Appl. Entomol. 131 (2007) 518. [28] L.M. Pecora, T.L. Carroll, Master stability functions for synchronized coupled systems, Phys. Rev. Lett. 80 (1998) 2109. [29] G.D. Ruxton, Density-Dependent migration and stability in a system of linked populations, Bull. Math. Biol. 58 (1996) 643. [30] N. Shigesada, K. Kawasaki, Invasion and range expansion of species: effects of long-distance dispersal, Dispersal Ecol. (2002) 350. [31] I. Shimada, T. Nagashima, A numerical approach to ergodic problem of dissipative dynamical systems, Prog. Theor. Phys. 61 (1979) 1605. [32] J.A.L. Silva, M.L. De Castro, D.A.R. Justo, Stability in a metapopulation model with density-dependent dispersal, Bull. Math. Biol. 63 (2001) 485. [33] J.A.L. Silva, J.A. Barrionuevo, F.T. Giordani, Synchronism in population networks with non linear coupling, Nonlinear Anal. Real World Appl. 11 (2010) 1005. [34] J.A.L. Silva, F.T. Giordani, Density-dependent migration and synchronism in metapopulations, Bull. Math. Biol. 68 (2006) 451. [35] H. Thunberg, Periodicity versus chaos in one-dimensional dynamics, SIAM Rev. 43 (2001) 3. [36] A. Trakhtenbrot, R. Nathan, G. Perry, D.M. Richardson, The importance of longdistance dispersal in biodiversity conservation, Divers. Distrib. 11 (2005) 173. [37] J. Ylikarjula, S. Alaja, J. Laakso, D. Tesar, Effects of patch number and dispersal patterns on population dynamics and synchrony, J. Theor. Biol. 207 (2000) 377. [38] J.B. Wenger, E.N. Naumova, Seasonal synchronization of influenza in the United States older adult population, PLoS ONE 5 (4) (2010) e10187, http:// dx.doi.org/10.1371/journal.pone.0010187.