## Abstract

This work explores an information-theoretic approach to drawing inferences about coupling of spatially extended ecological populations based solely on time-series of abundances. The efficacy of the approach, time‐delayed mutual information, was explored using a spatially extended predator–prey model system in which populations at different patches were coupled via diffusive movement. The approach identified the relative magnitude and direction of information flow resulting from animal movement between populations, the change in information flow as a function of distance separating populations, and the diffusive nature of the information flow. In addition, when the diffusive movement was eliminated from the model, mutual information correctly provided no evidence of information flow, even when population synchrony was generated by a common environmental driving function. Thus, for this model system, time‐delayed mutual information was useful in discriminating between the Moran effect and animal movement as causes of population synchrony, as well as in characterizing dispersal in terms of direction, relative speed and diffusive nature.

## 1. Introduction

Ecology is the study of the distribution and abundance of organisms over time and space. A central goal of ecologists is to develop a mechanistic understanding of temporal and spatial variation that can be used for prediction (e.g. Rhodes *et al*. 1996; Lande *et al*. 1999; Kendall *et al*. 2000; Keeling & Rohani 2002). Synchronous fluctuations of single-species populations in different locations at both local and regional scales are a form of time–space variation that has been of special interest to ecologists (e.g. Moran 1953; Ranta *et al*. 1995, 1998; Grenfell *et al*. 1998; Bjørnstad *et al*. 1999; Hudson & Cattadori 1999; Koenig 1999; Swanson & Johnson 1999; Ripa 2000; Peltonen *et al*. 2002; Post & Forchhammer 2002). Mechanisms underlying such synchrony fall into two basic categories (e.g. Ranta *et al*. 1998; Bjørnstad *et al*. 1999; Koenig 1999; Swanson & Johnson 1999; Kendall *et al*. 2000): (i) those based on movement of animals (by either the focal species or their predators) and (ii) those based on spatially correlated environmental variation (often referred to as the Moran effect; see Moran 1953; Royama 1992). When evidence of synchrony is found, distinguishing between these alternative explanations is important not only with respect to mechanistic prediction, but also from a conservation perspective (e.g. Hudson & Cattadori 1999). For example, local extinction of an area is much more likely to be followed by recolonization if dispersal is a primary cause of synchronization.

Studies of single-species population dynamics over space focus on both the possible existence of synchronous dynamics and, if present, the responsible mechanism(s). Evidence for the existence of synchrony typically involves cross-correlation analyses of time-series of abundance data (e.g. see Bjørnstad *et al*. 1999; Koenig 1999). Evidence for mechanisms underlying synchrony can be based either on direct approaches involving study of the purported mechanisms themselves or on indirect approaches based on the time-series of abundance data. For example, investigations can be focused on the relationship between local weather conditions and large-scale climate variation on vital rates and even body weights of natural animal populations (e.g. Leirs *et al*. 1997; Sillett *et al*. 2000; Mysterud *et al*. 2001). Similarly, the dispersal process has been the subject of extensive direct investigation (e.g. Clobert *et al*. 2001). Indirect approaches attempt to use time-series of abundances at different locations to discriminate between environmental variation and dispersal as causes of synchrony. For example, synchrony produced by dispersal is expected to decrease as a function of distance between locations, whereas this tendency is not expected if different locations are being influenced by the same environmental driving variable operating at a regional or global scale (e.g. Koenig 1999; Swanson & Johnson 1999). Nevertheless, there are limitations to what inferences can be made about mechanisms inducing synchrony from population time-series data alone (Bjørnstad *et al*. 1999).

Recently, Nichols *et al*. (2005) noted that indirect approaches to the study of coupling based on linear correlation were not necessarily appropriate for populations exhibiting nonlinear dynamics, and such approaches provided no inferences about asymmetric directional movement. Thus, they used methods based on attractor reconstruction for nonlinear systems to assess the strength and asymmetry of dispersal-based coupling between spatially separated components of a predator–prey model with diffusive movement (Pascual 1993; Little *et al*. 1996). Specifically, they investigated dynamic interdependence (generalized synchrony) by drawing inferences about the nature of functions relating reconstructed attractors. Using mutual prediction (Schiff *et al*. 1996) and continuity statistics (Pecora *et al*. 1995; Moniz *et al*. 2004), Nichols *et al*. (2005) were able to detect not only the dynamic interdependence of the components of this system but also the asymmetric nature of the coupling.

In this paper, coupling is viewed in terms of information flow and an information-theoretic approach, time-delayed mutual information (Vastano & Swinney 1988), is used to investigate coupling using time-series data. Specifically, time-series of location-specific abundances is generated using the spatially extended predator–prey model of Pascual (1993; also see Little *et al*. 1996; Nichols *et al*. 2005). The time-series data are used to determine the direction of information transport among locations and to infer the diffusive nature of the coupling. The diffusion terms are then removed from the model to generate time-series of isolated patches exhibiting predator–prey dynamics in the presence and absence of a common environmental driving function. Although the common environmental driver does induce synchrony of the populations, the time‐delayed mutual information provides no evidence of information transport. Thus, this data-analytic approach is able to clearly distinguish between the Moran effect (Moran 1953; Royama 1992) and dispersal as causes of population synchrony for this model system.

## 2. Time‐delayed mutual information

Assume that there exist two systems, *Y* and *Z*, and that we are able to measure these systems at *N* discrete points in time, resulting in the time-series *y*_{i}, *z*_{i}, *i*=1,…,*N*. Each measurement is a discrete random variable with underlying probability density function *p*(*y*_{i}) and *p*(*z*_{i}), respectively. The joint probability density (assuming that *Y*, *Z* have been measured simultaneously) is given by *p*(*y*_{i}, *z*_{i}). The amount of information common to these two systems can be described by the mutual information(2.1)Equation (2.1) may be interpreted as a metric reflecting the amount of information the process *Z* has in common with the process *Y*. If the two systems are assumed to be independent, then one has *p*(*y*_{i},*z*_{i})=*p*(*y*_{i})*p*(*z*_{i}), resulting in *I*(*Y*; *Z*)=0. If the process *Z* contains some information about the process *Y*, then this assumption is violated and *I*(*Y*; *Z*) will record some non-zero value of bits (assuming base 2 logarithms are used) common to the two processes. The mutual information function is symmetric under exchange of its arguments *y*, *z*, and cannot immediately be used to detect the direction of transport of information between the two systems. However, if one instead looks at the *time‐delayed* mutual information,then a sense of directionality is acquired. If information present at system *Y* at discrete time *i* is transmitted to system *Z*, then there will be a peak in the curve *I*(*Y*; *Z*_{T}) at *T* >0 as the joint probability density increases and reaches its maximum. A peak that occurs for *T*<0 implies that the information is being transported from *Z* to *Y*. Time‐delayed mutual information has been used to detect the direction of information flow in neuron firings (Destexhe 1994), in a reaction diffusion system (Vastano & Swinney 1988) and a coupled map lattice (Ho & Shin 2003). By utilizing a more generic definition of coupling (shared information), mutual information can be used to detect information transfer in linear or nonlinear systems. It should be mentioned that a different metric, the transfer entropy, was recently introduced by Schreiber (2000). The transfer entropy metric was specifically designed to look at information transport, and has already been used to examine physiological coupling (Kaiser & Schreiber 2002) and financial time-series (Marschinski & Kantz 2002). Transfer entropy utilizes the notion of conditional probability and is naturally asymmetric (as opposed to requiring a time delay). This particular metric quantifies the additional information in one time-series (process) that is not contained in the other. Transfer entropy is not explored in this paper but is referenced for completeness.

Computing the time‐delayed mutual information requires the estimation of single and joint probability densities as one system is time delayed with respect to the other. The method by which these quantities are estimated is predicated heavily on the amount and quality of the data being analysed. For a uni- or possibly multivariate time-series measurement **x**_{i} *i*=1,…,*N* (boldface type denotes a vector), one approach is to estimate the probability density at a given point *i* as(2.2)where the dependence of the estimate on the ‘bin size’ *R* is made explicit. Equation (2.2) represents a crude form of a kernel density estimation using the step kernel (Heaviside function):The integral quantifies the probability that a time-series point **x**_{j} will be found within some radius *R* of the point recorded at time *i*. Following the work of Liebert & Schuster (1989) and Prichard & Theiler (1995), the relevant quantities required by equation (2.1) can be given by(2.3)where *C*(* x*,

*R*) is the generalized correlation integral of order 1 (see Kantz & Schreiber 1999), obtained by averaging equation (2.2) over all points

*i*. Equation (2.1) is thereby altered to read(2.4)The biggest difficulty with using a non-parametric approach to analyse ecological data concerns the required time-series length. This general approach to density estimation requires that the time-series be of sufficient length to fully populate the partitioned data space. That is, the system must revisit a given state enough times to accurately quantify the local density. In this study and in previous work (see Nichols

*et al*. 2005), results were qualitatively the same for as few as 800–1000 point time-series even with this ‘brute-force’ approach to density estimation. For shorter datasets, the practitioner will probably have to alter the algorithm. One possibility is to make use of kernel density estimation techniques (see Silverman 1986), which use a more intelligent weighting scheme. Rather than simply counting points in an

*R*-ball, one may weight them according to their distance from the fiducial trajectory (

**x**_{i}in equation (2.2)). Parametric or semi-parametric approaches may also prove useful when there exists some

*a priori*knowledge of the distribution of the underlying populations. Also of interest is the situation where the practitioner has a model system (effectively unlimited data) and wishes to draw inferences from time-series generated by that model. Pascual & Levin (1999), for example, used time-series from a cellular automaton to assess a spatial scale for defining population densities while Wilson & Richards (2000) used a consumer-resource ‘patch’ model to analyse foraging in a spatial environment. Here, the potential utility of information theoretics is considered independent of the method of density estimation.

## 3. Numerical model

To examine coupling via information transport, a spatially one-dimensional predator–prey model was considered. The model was introduced by Pascual (1993) and explored further by Little *et al*. (1996), and describes the evolution of the prey *q* and predator *h* in both time *t* and space *x* according to(3.1)The dimensionless variables *q*, *h* and *x* represent the prey density, predator density and spatial coordinate, respectively. Reflective boundary conditions are considered at *x*=[0,1] asParameters are the predator death rate *m*, diffusion coefficient *d*, the predator/prey interaction *a*, the prey carrying capacity *b*, and the intrinsic growth rate *r*_{x}, of the prey population, which is (for non-zero *f*) a function of space. The resource term is a linearly declining function of space with value *r*_{0}=*e*=cost at the boundary. The parameter *f* governs the rate of resource decline and hence the degree of spatial asymmetry. As in Pascual (1993), the death rate, diffusion coefficient, boundary resource and carrying capacity are fixed at *m*=0.6, *d*=10^{−4}, *e*=5.0, *b*=2.0. Figure 1 shows the system response in both time and space. Oscillations near areas of high resource abundance tend to be periodic in nature while those near low resource areas exhibit chaotic dynamics (see Pascual 1993).

To simulate the Moran effect, an external driving term can be added in the form of a periodic function in the resource term resulting in(3.2)where *τ* is the period of oscillation and is taken to be roughly an order of magnitude longer than that associated with the system's natural period of oscillation. The diffusion coefficient is then set to 0 (*d*=0) to prohibit movement from one location to another. The system now represents a set of closed predator–prey communities influenced by a common environmental driving variable of amplitude *ϵ*.

## 4. Application to a predator–prey system

Using a finite difference scheme, equation (3.1) was integrated for *N*=10 000 time-steps (post-transient) with a dimensionless sampling time of 0.5. The dominant frequency for this system is approximately 0.065 cycles per unit time so that this sampling rate was deemed sufficient to capture the dynamics. Simulations were performed on a spatial grid consisting of *j*=1,…,100 sites, distributed evenly on *x*∈[0, 1]. Prior to the analysis, all time-series were normalized by subtracting their mean and dividing by the standard deviation. Simple changes in population mean and/or variance are therefore removed from consideration here.

In an earlier work by Nichols *et al*. (2005), it was hypothesized that information in this system flowed ‘downhill’ from areas of high resource abundance to areas of low resource. Put another way, the dynamics of the populations at the high resource end are influencing those at the low resource end but not vice versa. The main reason for this conjecture is that the populations at the higher resource end of the system are simply larger. Diffusive coupling links a given site on the lattice to its two neighbours (the numerical approximation to the second derivative includes the *j*−1, *j*, *j*+1 lattice site populations). Higher populations sizes for smaller *j* means that the term linked to the (*j*−1)th lattice site will more heavily influence the dynamics at the *j*th site than will the term at site *j*+1. Thus, the dynamics at the high resource end propagate towards the low resource end. This hypothesis may be directly tested using the concept of time‐delayed mutual information.

To this end the mutual information between the prey dynamics (time-series) at various spatial points at the high resource end and a ‘target’ point at the low resource end was computed. Specifically, the function *I* (*Y*^{(j)}, , *R*), *j*=70, 71,…,94, *k*=95 was computed for−350≤*T*≤350 and *R*=0.035 (3.5% of the signal variance). The notation *Y*^{(j)} simply denotes that the system *Y* was observed at spatial location *j*. Results are shown in figure 2 for select spatial sites. Each curve shows a single peak on which a number of periodic fluctuations are superimposed. These fluctuations correspond to the natural period of the dynamics. Essentially, the time-series are going in and out of phase as one is delayed with respect to the other. The amplitude of these fluctuations is larger than shown in figure 2. A 25-point ‘running smooth’ was used in which a sliding window of 25 points was passed through the curve and the mean of each window plotted. This was done purely for visualization purposes (otherwise, the fluctuations on the various curves overlap such that it becomes difficult to see the separate peaks). The important feature is the location of the dominant peak in time. All peaks occur at *T* >0, indicating that information flows from the high resource locations to location *x*=0.95. As expected, a larger spatial separation implies a longer time for the information to reach the target point. It should be noted that the time-scale of information transport is significantly longer than the natural period of the system.

Plots such as figure 2 may also be used to define critical length (distance) scales of influence. It is clear that once the spatial separation reaches Δ*x*=0.25 (e.g. lattice sites *x*=0.7 and 0.95), the prey populations have low mutual information (approximately 0.2 bits) and little information exchange (the curve is nearly flat).

In addition to the directionality information, one can also deduce the form of the coupling (i.e. the form of the dispersal functional) by examining the spread of the peaks. Following Vastano & Swinney (1988), the width of the curves can be quantified by the ‘full width half-max’ (width of the peak at half it's maximum value), denoted Δ*T*. Figure 3 shows the log–log plot of the Δ*T* values as a function of spatial separation Δ*x*. For diffusive coupling, it has been observed (Vastano & Swinney 1988; Destexhe 1994) that the width of the peaks tends to grow as the distance squared from the target point. As shown in figure 3 the slope is approximately 0.5 suggesting that indeed the coupling is diffusive in nature.

As a check, the simulations were repeated with the coupling removed (*d*=0.0 in the model). Figure 4 shows the resulting *I*(*Y*^{(90)}, , *R*) curves corresponding to spatial locations *x*=0.90, 0.95. The radius used in estimating equation (2.4) was fixed at *R*=0.035. As with figure 2, a 25-point running smooth was used in plotting in order to maintain consistency. The absence of any clear peak (asymmetry) illustrates that no information is passing from one spatial site to another. All that is being observed are fluctuations at the system's natural period as the delayed measurement goes in and out of phase with the fixed measurement *Y*^{(90)}. The amplitude of the fluctuations is also seen to be much lower than with the diffusive coupling added. Furthermore, the average number of bits *Y*^{(90)} and have in common is much lower (approximately 0.1 bits) than when diffusive coupling is present. However, because the magnitude of *I*(*Y*, *Z*_{T}, *R*) can be heavily dependent on *R*, care should be exercised when drawing conclusions based on scale alone (a discussion of the effects of *R* may be found in Kaiser & Schreiber (2002)).

The final test was to add in an external periodic driving term while the spatial coupling remained fixed at *d*=0.0. The goal here was to discriminate the synchronizing effects of external driving from dispersal. Simulations were repeated with the inclusion of the periodic forcing alluded to in §3 (equation (3.2)) using an amplitude of *ϵ*=0.25 at a frequency of 1/*τ*=0.01 cycles per unit time. Results are shown in figure 4 where the probabilities were again estimated using *R*=0.035. In this case, the dominant period of the driving is evident in the information content, yet there is no directionality associated with any information flow. Unlike figure 2, there is no dominant peak in either forward or reverse time indicating the absence of information transport. The conclusion is that although both systems (lattice sites) are ‘slaved’ to a common driver there is still no evidence of information transport between them. It may therefore be concluded that this approach may be used to distinguish the Moran effect from dispersal. However, if the drive signal is spatially varying (not global), then such a conclusion may be difficult to reach.

## 5. Conclusions

Drawing inferences about mechanisms underlying population synchrony from time-series of abundances is a difficult problem (e.g. Bjørnstad *et al*. 1999). In another work, Nichols *et al*. (2005) drew inferences about coupling based on general characteristics of the functions relating two reconstructed attractors. In the present work, inferences are based on an information-theoretic approach in which dynamic systems are viewed as generators of information (e.g. Prichard & Theiler 1995). The use of information theory in ecology is not new and ranges from species diversity metrics and inferences about diversity–stability relationships (e.g. see MacArthur 1955; Margalef 1958, 1968) to natural selection within age‐structured populations (Demetrius 1974, 1975), to model selection in the context of statistical inference (Burnham & Anderson 2002). The work in this paper uses the time‐delay mutual information approach of Vastano & Swinney (1988), but other information-theoretic approaches appear to hold promise as well (Schreiber 2000; Palûs *et al*. 2001; Kaiser & Schreiber 2002). All of these approaches require only that the practitioner be able to estimate the relevant probability densities from the data.

The approach used in this paper sought to identify the time delay at which the mutual information function between two time-series reaches a maximum. In the presence of diffusive movement within the model predator–prey system, the approach identified the relative magnitude and direction of information flow resulting from animal movement between populations, the change in information flow as a function of distance separating populations, and the diffusive nature of the information flow. In addition, when the diffusive movement was eliminated from the model, mutual information correctly provided no evidence of information flow, even when population synchrony was generated by a common environmental driving function. Thus, for this model system, time‐delayed mutual information was useful in discriminating between the Moran effect and animal movement as causes of population synchrony, as well as in characterizing dispersal in terms of direction, relative speed and diffusive nature.

## Acknowledgments

The author acknowledges the US Geological Survey for publication costs. The author also thanks Jim Nichols of the US Geological Survey for helpful discussions regarding the ecological interpretations of information transport and for helping to motivate this work. Thanks also to Evan Cooch of Cornell University for reading this manuscript and offering helpful suggestions.

## Footnotes

- Received May 28, 2004.
- Accepted August 2, 2004.

- © 2005 The Royal Society