Synchrony in Patch Models: How Topology can Distinguished Spatial Scales.

Timothy Reluga

Date: October 30, 2001
Begun Monday, 30 April, 2001

Abstract:

One of the issues in bridging the gap in spatial modeling between local and global dynamics is understanding the scale transition. To this end, I present a general discuss of sufficient conditions for the existence of synchronous dynamics in single-species coupled map lattices, and apply this analysis to lattices of spatial homogeneity and spatial heterogeneity. The results provide an intuitive basis for the minimal spatial scales at which mean-field approximations are biologically applicable.

1 Introduction

One of the key problems in modern theoretical ecology is elucidating the relationship between space and population dynamics. The two standard mathematical approaches to space are the homogeneity and the mean-field postulates. The homogeneity postulate assumes that populations are well-mixed and obey mass-action laws analogous to those of chemical systems. The mean-field postulate assumes independent behavior among spatially allocated subpopulations, allowing the population to be described in terms of global average. These postulates apply at distinct spatial scales. The mean-field postulate is useful at the largest, global spatial scales, where distances are so large that migration is insignificant. The homogeneity assumption is useful at an intermediate ``patch'' scale, where distances are great enough that the motion of individual organisms appears random, but distances are not so great as to limit mixing. There are two other spatial scales which are believe to play an important role in ecological systems. At the smallest ``local'' spatial scales, much work has been done on individual-based models where each organism and its interactions are modeled in detail. The other, and one which to be addressed here, is the ``metapopulation'' scale. The ``metapopulation'' scale is a spatial scale between the ``patch'' scale and the ``global'' scale, where distances are great enough to prevent the population from being well-mixed, but not great enough to make patches behavie independently.

The metapopulation scale has been an area of recent interest. One reason for the interest is an appealling explicit appearance of topology and geometry of source and sink patches in the models. This provides an excellent opportunity to explore the transition in behavior from the homogenous to the mean-field postulates. Using recent progreess in coupled lattice map theory, I compare the spatial dynamics in 3 common topologies and discuss the relationship between migration and emergent synchronous behavior. Sections 2 and 3 summarize and expand the results of Silva et al.[SCJ00] concerning synchrony in a coupled lattice map. Section 4 summarizes standard results concerning synchrony in 2-patch, uniformly coupled, and ring lattices. Section 5 performs a similar analysis on the heterogenous topology of a binary tree. Section 6 compares analytic and numeric results concerning synchrony in each topology. Section 7 discusses the significance for these results.

2 General Model and Theory

Silva, De Castro, and Justo first investigated synchrony in a coupled map lattice of arbitrary size[SCJ00] under density-independent migration. Here I present a more general analysis which encompasses density dependent migration.

A single species coupled lattice framework consists of an initial lattice configuration $\vec{x}_{0}\in \Re^{n}$, and a propagating equation

$$\vec{x}_{t+1}= F(G(\vec{x}_{t})),$$ (1)

where $F$ is some migration operator, and $G$ is some local growth operator. For notational convience, the vector notation will henceforth be omitted. Either of these operators may be non-linear. The miminal constraints we usually impose are that $F():\Re ^{n+} \times \Re^{n+}$ is dissipative( $\vert\vert F(x)\vert\vert _{1} \leq \vert\vert x\vert\vert _{1}$) and that $G():\Re ^{n+} \times \Re^{n+}$ is diagonal( $\partial G_{i} / \partial x_{j} = 0$ if $i \neq j$).

Our intuitive definition of a synchronous lattice configuration needs a formalization. Let $e$ be the vector of all $1$'s, and let $\Lambda_{1} = \{ a \vec{e} : a \in \Re^{+} \}$. A lattice configuration $x$ is then synchronous if $x \in \Lambda_{1}$. This a stronger definition than would be produced by the idea of coherence, and coherence may well be a more significant biological concept. Nevertheless, we need this full strength to procede mathematically.

Determining when a metapopulation can behavie synchronously is equivalent to determining when $\exists \Omega \subseteq \Lambda_{1}$ in the neighborhood of which, $F(G())$ behavies like a contraction. This can be broken down into two requirements:

$$F(G(\Lambda_{1})) \subseteq \Lambda_{1},$$ (2)

$$\exists r : \forall \vert\vert\epsilon\vert\vert \sim 0, \epsilon... ...\nabla F(G(x))\nabla G(x) \epsilon, \vert\vert / \vert\vert \epsilon \vert\vert$$ (3)

Condition (2) is not true in general, and must be handled in our definitions of $F()$ and $G()$. It is sufficient that $G()$ be isotropic( $P^{-1}G(Px) = G(x)$ for any permutation matrix $P$) and $F(\Lambda_{1}) \subseteq \Lambda_{1}$. These are restrictive, but not unreasonable. The critical problem is to determine when synchronous configurations are stable, according to condition (3). Normally, one would look at the spectral radii of each of the linearized operators, and use the submultiplication property of matrix norms to provide an minimum upper bound on stability.

$$\max{\left(\rho(\nabla F(G(x))) \right)} \max{\left(\rho(\nabla G(x)) \right)} \leq 1$$ (4)

While this is still applicable, it is not satisfactory. In its weakest form, it says nothing more than that a spatially stable homogenous equilibria is synchronus. Since we are investigating the stability of a subspace, $\epsilon$ is choosen in a rank-deficient manner which can significantly weaken the requirements for stability. In addition, the spectral radii realized during iteration may be significantly lower than the theoretical upper bounds.

3 Density-Independent Migration Theory

Having run into a roadblock, we cannot precede without placing additional constraints on $F()$. Let us make the simplest reasonable assumption: conservative, density-independent migration.

$$F(x) = Mx,$$ (5)

where $M$ is a symetric stochastic matrix with a unique dominate eigenvalue $\lambda_{1}=1$ corresponding to subspace $\Lambda_{1}$. Denote the other eigenvalues of $M$ as $\lambda_{1} > \lambda_{2} > \lambda_{3} > \ldots \geq 0$. The symetry of $M$ means there is no directional bias in migration between connected patches. Additionally, let

$$\rho(G^{*}) = \sup{ \left\{ \frac{ \vert\vert\nabla G(x) \epsilon \vert\vert }{ \vert\vert\epsilon\vert\vert } : x \in \Omega \right\} }$$ (6)

Theorem 1   For any operators $M$ and $G$ meeting the above restrictions and

$$\vert\lambda_{2}\vert \vert\rho(G^{*})\vert < 1,$$

$\Omega \subseteq \Lambda_{1}$ is an attracting manifold.

Proof. This is a special case of more general spectral radius theorems. The key point is that, since the $\epsilon$ in Condition (3) is orthogonal to $\Lambda_{1}$, the dominate component of $M$ appears in the form of $\lambda_{2}$. $\qedsymbol$

The problem of synchrony is reduced to the production of practical bounds on $\lambda_{2}$, which depends on the lattic topology, and $\rho(G^{*})$. It has already been assumed that $G()$ is a diagonal isotropic operator, so it can be writen as $G(x): x_{i} \mapsto x_{i}g(x_{i})$ for some $g(x)$. The equations of Ricker and Hassell are often used as representations for the local dynamics in meta-population models. Either of these equations could be used here, but I will instead use a mathematically convient form of the Smith-Slatkin equation:

$$g(x) = \frac{2\beta}{\beta-1} \frac{1}{1+x^{\beta}} , \beta > 1.$$ (7)

This 1-parameter formulation has a bifurcation diagram that is qualitatively equivalent to that of the Ricker equation as $\beta \rightarrow \infty$. Standard analysis shows that there is a unique positive steady state, $x_{eq} = \left( \frac{\beta+1}{\beta-1} \right)^{1/\beta},$ which is attracting for $\beta < 3$. In addition it is a straight-forward manipulation to show that

$$\forall x \geq 0, f'(x_{eq}) \leq f'(x).$$

The coincidence of the steady state and an inflection point will simplify the analysis which is to follow.

Let

$$\gamma = \sup{\left\{ \left\vert \frac{\partial}{\partial x} \left( x g(x) \right) \right\vert : x \in R^{+} \right\}}$$ (8)

$$\gamma^{*} = \left\vert \inf{ \left\{ \frac{\partial}{\partial x} \left( x g(x) \right) : x \in R^{+} \right\}} \right\vert$$ (9)

Using Eq. (7), $\gamma$ can be defined peicewise as

$$\gamma = \left \{ \begin{matrix} \frac{2\beta}{\beta-1} & \textno... ...\frac{\beta - 1}{2} & \textnormal{if} & \beta > 3+2\sqrt{2} \end{matrix}\right.$$ (10)

while $\gamma^{*}=\left\vert \frac{1-\beta}{2} \right\vert$

It is trivial to show $\rho(G^{*}) \leq \gamma$. I also contend

Conjecture 1   $\rho(G^{*}) \leq \gamma^{*}$

Proof. I don't know why this is so, but it does seem so, experimentally, for the Ricker, Hassell, and Smith-Slatkin models. $\qedsymbol$

4 Standard Examples

Gyllenberg et al.[GSE93] and Hastings[Has93] analyze the 2-patch coupled logistic model in depth. Sole and Gamarra specifically addressed the idea of synchrony in a coupled-patch model[SG98].

$$\left[ \begin{matrix}x_{t+1} y_{t+1} \end{matrix} \right] = \l... ...{matrix} \right] \left[ \begin{matrix}g(x_{t}) g(y_{t}) \end{matrix} \right]$$ (11)

with $g(x)$ as defined in Eq. 7. This fits the coupled lattice model described above. Since $0<m<1$, $\lambda_{1}=1, \lambda_{2}=1-2m$. Sole and Gamarra showed that the synchrony condition here is

$$\left\vert (1-2m) e^{\mathcal{L}(g)} \right\vert < 1.$$ (12)

This result is very similar to the standard homogeneous stability result derived in []. There is a horn of stability around $m=2$ outside of which synchronous configurations can not persist. This horn, of course, contains the domain of homogenous stability as a proper subset.(see Fig. 1).

Hastings[Has93] and Gyllenberg et al.[GSE93] discuss the stability of periodic orbits.

The two patch model is often generalized to encompass $N$ uniformly connected patches. In this case,

$$M_{ij} = \left\{ \begin{matrix} 1 - (N-1) m & \textnormal{if}& i=j \\ m & \textnormal{if}& i \neq j \end{matrix}\right.$$ (13)

with the requirement that $0<m<1/(N-1)$. Then, $\lambda_{2} = 1-Nm$, with multiplicity $N-1$, $\Lambda_{2} \sim [1,0,0,..,0,-1,0,0,..]$. The synchrony condition

$$\left\vert (1-Nm) e^{\mathcal{L}(g)} \right\vert < 1,$$ (14)

applies. Notice that all the eigenvectors lose stability at the same time.

The simplest topology with non-trivial geometry is that of a circular ring of patches with nearest neighbor coupling. This is the case actually addressed by Sole and Gamarra[SG98]. Consider a set of $N$ patches, $N$ even, with indexes $i \in Z^{N}$. Define

$$M_{ij} = \left\{ \begin{matrix} 1 - 2 m & \textnormal{if}& i=j \\... ...al{if}& \vert i - j\vert = 1 \\ 0 & \textnormal{otherwise} \end{matrix}\right.$$ (15)

Note that under this definition, $M_{0,N-1} = M_{N-1,0} = m$. Since we are only considering nearest-neighbor migration, overdispersive artifacts can be avoided by restricting $0<m<1/4$.

This topology exhibits a full spectrum of eigenvalues.

$$\lambda_{i} = 1 - 2m + 2m\cos{ \left[ \frac{2\pi(i-1)}{N} \right] }$$ (16)

$$\Lambda_{i} = e^{2\pi(i-1)/N}$$ (17)

Sole and Gamarra show that a condition for synchrony here is

$$1-2m(1-\cos{\frac{2\pi}{N}}) < e^{-\mathcal{L}(g)},$$ (18)

but it should be noted that eigenvalues of higher frequency lose stability at smaller migration rates.

5 Trees: Heterogenous Topologies

Having presented the standard models, one is struct by the uniform distribution of their spatial structure. It is prudent to have some clumpier arrangements of patches for comparison. For instance, a river system will exhibit a very tree-like structure of suitable habitats through which migration occurs. None of the topologies discussed about can be reasonably applied to this situation.

Because of their key role in computer science, tree topologies have become a common site in modern science. Assume we have a some rooted tree with $N$ leaves where all the branches have been given lengths. A density-independent migration matrix can be derived from this tree in several ways. One simple method in keeping with the river-system analogy proposed above is as follows:

Imagine a river system feed by a set of small lakes. These lakes are the only reproductively suitable habitat for a species of fish. The lakes drain into a river system which provides good foraging habitat during the summer, but no reproductive opportunities. Every spring, after reproduction, individuals migrate down-stream towards the river mouth with some redistribution kernel $k(x)$, $x$ being the distance downstream from their lake of origin. The root branch(the ocean) is assumed to have infinite length, so there need be no explicit restrictions on $x$. The following fall, individuals migrate back up the tree(river system), making random decisions at every fork they encounter until they reach a leaf. The following spring, reproduction occurs in the lakes, and the cylce repeats.

This model can produce a wide range of migration matrii, conservative or dissipative, some of which violate the standard homogeneity stability criteria[SCJ01] and [JM00] as well as the assumptions of this paper. I shall only analyze one of the simplest cases.

Consider a binary tree of $2^{N}$ nodes, where all branches and leaves have length $\ell$, and an exponential redistribution kernel $k(x) = e^{-r x}/r$. Migration is riskless and unbiased upon the return. Applying the above model(see Appendix A),

$$M_{i,j}(N,p) = \frac{1}{1-p} p^{N+1} + \left( \frac{1-2p}{1-p} \right) p^{d(i,j)}$$ (19)

where $p = e^{-r\ell}/2$ and $d(i,j)$ is minimum tree-distance between $i$ and $j$. In this special case, $M$ meets all the necessary conditions for the application of Theorem 1. In the limit $p \rightarrow 1/2_{-}$, the uniform coupling model is recovered. Of special interest is a recursive block construction of $M$. Let $E_{n}$ be the $n \times n$ matrix of all $1$'s.

$$M = A(N) + \frac{p^{N+1}}{1-p} E_{2^{N}}$$ (20)

where $A(N)$ is defined recursively as

$$A(N) = \left[ \begin{matrix} A(N-1) & \mu(N) E_{2^{N-2}} \\ \mu(N) E_{2^{N-2}} & A(N-1) \end{matrix} \right]$$ (21)

$$A(0) = \mu(0)$$ (22)

$$\mu(x) = \frac{1-2p}{1-p} p^{x}$$ (23)

This leads to an elegant result.

Theorem 2   If $\Lambda$ is an eigenvector of $A(x)$ summing to 0, with eigenvalue $\lambda$. Then, $[\Lambda,\pm \Lambda]$ are both eigenvectors of $A(x+1)$ also with eigenvalue $\lambda$.

Proof. Do I have to? Not if I cite somebody. Notice that the proof is independent of $\mu(x)$. This means different levels of the tree could have different branch lengths. $\qedsymbol$

$[ 1, -1]$ and $[1,1]$ are eigenvectors of any symetric $2 \times 2$ matrix, specifically $A(1)$ and $E_{2}$. $e$ is an eigenvector of $A(x)$ for every $x$. Thus,

$$\lambda_{i} = \left\{ \begin{matrix}1 & \textnormal{if} & i = 1 1 - p^{N-i+2} & \textnormal{if} & 2 \leq i \leq N+2 \end{matrix} \right.$$ (24)

The eigenvectors can be explicitly constructed by concatination and negation.

6 Results

Experimental evidence suggests $\lambda_{\lceil N/2 \rceil}$ is actually the eigenvalue controlling synchrony. $\lambda_{2}$ has the low frequency of $2\pi/N$, while $\lambda_{\lceil N/2 \rceil}$ has the highest frequency eigenvector of all. If $\Lambda_{ \lceil N/2 \rceil }$ perturbations decay, all high-frequency components will eventually force out low-frequency perturbations as well. This is a concept beyond orthagonality, and plays a role in multigrid solutions of Laplace's equation. The experimentally observed requirement for synchrony is the requirement for local synchrony that

$$\left\vert 1 - 4m \right\vert < \frac{2}{\beta-1}$$ (25)

7 Discussion

Potentially disipative nature should increase stability further.

Despite this model's spatial nature, it still focuses very much on the homogenous space situation where all local patches have identical dynamics. Much of the analysis hinges on the unlikelyl isotropic nature of $G()$. Attempts to extend this theory must be able to address the concept of coherence as a generalization of synchrony to address a truer form of spatial heterogeneity.

Future work may include disipative migration matrixes, a generalization of lyaponov exponents(determinates?) to address coherence, and the effects of noise on stability.

Bibliography

GSE93

Mats Gyllenberg, Gunnar Soderbacka, and Stefan Ericsson.
Does migration stabilize local population dynamics? analysis of a discrete metapopulation model.
Mathematical Biosciences, 118:25-49, 1993.

Has93

Alan Hastings.
Complex interactions between dispersal and dynamics: Lessons from coupled logistic equations.
Ecology, 74(5):1362-1372, 1993.

JM00

S. R.-J. Jang and Arun K. Mitra.
Equilibrium stability of single-species metapopulations.
Bulletin of Mathematical Biology, 62:155-161, 2000.

SCJ00

Jacques A. L. Silva, Manuela L. De Castro, and Dagoberto A. R. Justo.
Synchronism in a metapopulation model.
Bulletin of Mathematical Biology, 62:337-349, 2000.

SCJ01

Jacques A. L. Silva, Manuela L. De Castro, and Dagoberto A. R. Justo.
Stability in a metapopulation model with density-dependent dispersal.
Bulletin of Mathematical Biology, 63:485-505, 2001.

SG98

Ricard V. Sole and Javier G. P. Gamarra.
Chaos, dispersal and extinction in coupled ecosystems.
Journal of Theoretical Biology, 193:539-541, 1998.

A. Derivation of Binary Tree migration

Consider a binary tree connecting $2^{N}$ patches, where all branches and leaves have length $\ell$. Let $\zeta(n)$ be the probability that an individual migrates a distance between $n\ell$ and $(n+1)\ell$, if $n<N$. Define $\zeta(N)$ as the probability of migrating a distance greater than $N\ell$.

$$\zeta(n) = \left\{ \begin{matrix} \int_{n\ell}^{(n+1)\ell} k(x)dx... ...<N \\ \int_{N\ell}^{\infty}k(x)dx & \textnormal{if} & n=N \end{matrix} \right.$$ (26)

By specifying an exponential migration kernel $k(x)=(\exp{-\lambda x})/\lambda$, and defining $p=e^{-\lambda \ell}/2$

$$\zeta(n) = \left\{ \begin{matrix} (2p)^{n}(1 - 2p) & \textnormal{if} & n<N \\ (2p)^{n} & \textnormal{if} & n=N \end{matrix} \right.$$ (27)

Now consider two patches $i$ and $j$ separated by distance $2d\ell$ on the binary tree. Then

$$M_{i,j}(N,p) = \sum_{k=d}^{N} 2^{-k} \zeta(k)$$ (28)

$$= 2^{-N}(2p)^{N} + \sum_{k=d}^{N-1} 2^{-k} \left( (2p)^{k}(1 - 2p) \right)$$ (29)

$$= p^{N}+\sum_{k=d}^{N-1} p^{k} (1 - 2p)$$ (30)

$$= \sum_{k=d}^{N} p^{k} - 2p\sum_{k=d}^{N-1} p^{k}$$ (31)

$$= p^{d} - \sum_{k=d+1}^{N} p^{k}$$ (32)

$$= p^{d} - p^{d+1} \frac{ 1 - p^{N-d}}{1-p}$$ (33)

$$= \frac{ p^{d} - 2 p^{d+1} + p^{N+1}}{1-p}$$ (34)

$$M_{i,j}(N,p) = \frac{1}{1-p} p^{N+1} + \left( \frac{1-2p}{1-p} \right) p^{d}$$ (35)

B. Minimum Coherence requirements

Patches	Minimum Requirement
2 Patch	$1-2m < 1 / \rho(G^{*})$
Complete connection of $N$ patches	$1 - N m < 1 / \rho(G^{*})$
Linear connection of $N$ patches	$1-4m < 1 / \rho(G^{*})$
Binary branching of $2^{N}$ patches	$1 - p^{N} < 1 / \rho(G^{*})$

C. To Do

• Phase in terms of ultrametric spaces, and Bethe lattices which are already well known objects in mathematical physics. The ultrametric space seems a little more appropriate, since all the patches in the tree model do satisfy a strong triangle inequality.

• Draw connection to phase-transitions in physics, where patchiness can be observed near threshold temperatures.

• Is the smallest eigenvalue of a doubly stochastic/conservative matrix greater than 0? Is there a better class of matri to choose?

• Verify that the second largest eigenvalue of a doubly stochastic matrix is less than 1. Ergoticity should be enough to show this. I may just have to assume the desired eigenstructure if I cannot show these properties.

• Deal with the problem of large positive derivatives. This complicates the analysis of the problem because I have to explain why I care about the negative derivative values, and not the positive derivative values.

• Incorporate Silva's ring model with lyaponov exponent formulation of synchrony.

• Cite Kaneko'93 for framework(refence in Silva00).

• Look at work on shuffling problem discussion of Trefethen and Diaconis.

• Get Keeling reference out of recent Theoretical ecology book for importance of metapopulation work. Chapter 3 in McGlade, "Spatial models of interacting populations", p 64-99. He specifically raises the question of what the size of homogenious cells should be in simulations.

• Address leakage to other fourier modes as the reason that the higher frequencies control global synchrony stability?

• Numerics summarizing the loss of synchrony in uniform, ring, and tree models.

• Address the alternate model formulation, where migration comes first, followed by survivorship. Are they equivalent? Are they different? Why use one instead of the other?

1 Migration in a potential

A standard result in coupled-lattice theory is that migration density independent migration does not destabilize homogenous steady-states. When migration is biased by some physical asymetry like river-flow, prevailing winds, or even gravity, this does not hold. In fact, such migration will even change the steady state.

$$x_{t+1} = MF(x_{t})$$ (36)

If $\Lambda_{1} \not \owns 1_{H}$ , then, even if $F()$ is diagonal isotropic, $F(x_{eq})=x_{eq} \not \rightarrow M F(x_{eq}) = x_{eq}$. Put another way, $F()$ is not necessarily isotropic, since different patches probably have different qualities. Define $x_{eq}: F(x_{eq}) = x_{eq}$. Under extension to a coupled lattice with density-independent migration, $MF(x_{eq}) = x_{eq} \iff x_{eq} \in \Lambda_{1}$

Figure 1: The coherence horn of a 2-patch system. Note that $m=1/2$ is the most stable coupling, and that synchrony disappears as $\beta \rightarrow \infty$.

About this document ...

Synchrony in Patch Models: How Topology can Distinguished Spatial Scales.

This document was generated using the LaTeX2HTML translator Version 2K.1beta (1.48)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -dir Patches -split 0 -show_section_numbers Paper.tex

The translation was initiated by You on 2001-10-30