Detecting nonlinear dynamics in spatio-temporal systems
Ph'YSIOA ELSEVIER
Physica D 96 (1996) 321-333
Detecting nonlinear dynamics in spatio-temporal systems, examples from ecological models Sarah Little a,*, Stephen Ellner b, Mercedes Pascual a, Michael Neubert a, Daniel Kaplan c, Timothy Sauer d, Hal Caswell a, Andy Solow a a Woods Hole Oceanographic Institution, Wood Hole, MA 02543, USA b North Carolina State University c McGill University d George Mason University
Abstract
Mathematical models of marine populations exhibit chaotic dynamics. However, we hypothesize that in moving water, Eulerian sampling of spatially heterogeneous populations may obscure any deterministic signal beyond the resolving capabilities of presently available nonlinear signal processing techniques. To examine this hypothesis we created two spatio-temporal models of population dynamics. To caricature actual ocean sampling limitations, we sampled the model output in two ways, random walks to simulate Eulerian sampling, and spatial averages to simulate population measurements from finite volumes. Results indicate that the ability to identify underlying nonlinear dynamics quickly degrades as the step size of a random walk sampling increases. On the other hand, the analysis techniques used are more robust in the face of spatial averaging. Keywords: Nonlinear dynamics; Marine ecology; Spatio-temporal chaos; Sampling
1. Introduction
In this paper we examine how noise introduced by limitations in sampling of a chaotic spatio-temporal system affects the ability to decipher underlying nonlinear dynamics. This study arose out of work on sampling of marine ecosystems, but the results are extendable to spatio-temporal systems in general. In the ocean, complex spatio-temporal interactions between predators, prey, and nutrients, coupled with ocean physics, govern the population dynamics *Corresponding author. Fax: +1-508-289-3245; e-mail: [email protected] (508)289-3245.
of everything from phytoplankton to whales. Understanding the behavior of populations of oceanic organisms has great practical significance ranging from managing commercial fisheries such as cod and salmon, to controlling blooms of toxic plankton. However, the simple act of taking meaningful, long-term, spatio-temporal measurements of populations in a moving fluid is a formidable task itself. First, all point samples of a population must be spatial averages of some sort. Second, long-term ocean sampling is most easily accomplished in the Eulerian sense, that is, in a reference frame fixed with respect to the moving fluid, such as from a moored buoy or a pier. For an externally advected population, this effectively means
0167-2789/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved PII S0167-2789(96)00030-9
322
S. Little et al./Physica D 96 (1996) 321-333
that temporally successive samples are not taken from the same spatial location with respect to the spatial pattern of the population. The effects of these two sampling limitations on the ability to detect nonlinear deterministic behavior are addressed in this paper by applying several nonlinear time series analysis techniques to the output of two different spatio-temporal population models. In theoretical ecology, there are many examples of temporal population models which exhibit chaos. The interaction of three variables in a predatorprey-nutrient system [1] is a well studied chaotic system, as is the simple logistic map [2]. Examples of spatio-temporal models are less common, but include population models in discrete time and space [3,4], and discrete time and continuous space [5]. The tests conducted in this paper use a onedimensional continuous time and space model of a predator-prey-resource interaction [6], and a onedimensional discrete time and space model of a single species; both models exhibit chaotic dynamics. Time series representing spatial averaging are generated by taking spatial averages of constant length and location from the output of the modes. Time series which caricature Eulerian sampling are generated by sampling at points determined by a random walk in one-dimensional space. These series are compared to a time series of samples from a single spatial location. The tests used in this paper utilize several recently developed methods [7-9] which are intended for situations where the data are scarce, the functional form of the underlying process is unknown, and noise may be present. The methods assay different aspects of nonlinear dynamics - predictability, continuity, and sensitive dependence on initial conditions. The purpose of this paper is not to compare methods, but to determine the relative effects of Eulerian sampling vs. spatial averaging on our ability to detect that the underlying dynamics are deterministic spatio-temporal chaos. Because the methods are new and their performance on spatio-temporal data has not been examined, we chose several methods so that consistent answers from all methods will give us confidence in the results.
These methods, and most other nonlinear dynamics signal processing ones, are based on phase space reconstruction using delay coordinate embedding V(ti) ---- (xi, xi+r, x i + 2 r . . . . .
Xi+(d-l)r),
(1)
where V is the embedded vector, x(t) the measured variable, d the embedding dimension, r the time lag, and i the index for time t. Here the methods diverge, [8,7] tests approach the problem of detecting nonlinear determinism with a two step process. First, they develop a measure of the data which will be sensitive to determinism. Their techniques capitalize on the geometric structure of the attractor in phase space which is the hallmark of deterministic systems. On the one hand, the complexity of this geometric structure is determined by the dynamics, and on the other, noise will tend to destroy the smoothness of the geometric structures. The techniques for measuring determinism here work by looking at the smoothness of local neighborhoods in phase space of the reconstructed attractor, the idea being that the smoother the structure, the more likely that it arose from a deterministic process and not a stochastic one. Secondly, they test whether the determinism arises linearly or nonlinearly. They use the determinism measure to compare the data with linear surrogates which have the same power spectrum as the data but have a random phase spectrum [10]. The surrogates are representative of a linear random process with the same autocorrelation as the data. If a difference is found between the data and their linear surrogates, then the conclusion is drawn that the data do not arise from a linear stochastic system. Importantly though, the test does not rule out the situation where a deterministic process underlies the data but the data have a large noise component, nor does it rule out that the data arose from a nonlinear stochastic process. The Nychka, Ellner, Gallant, McCaffrey test ([9,11]) (hereafter NEGM) approaches the question of detecting nonlinear determinism by measuring the Lyapunov exponent of a model fitted to the data. Because ecological time series typically have fewer than 1000 points, it is impossible to reliably estimate Lyapunov exponents from data directly. The approach generally taken, then, is to fit a nonlinear model to
S. Little et al./Physica D 96 (1996) 321-333 the data and calculate the exponents directly from the model. The main pitfall of this approach is that the results can depend entirely on the choice of model. The issues of fitting an appropriate model to a data set are not trivial, and have been addressed in [9]. Analyses of population dynamics data have often been based on fitting the data with simple models having a few adjustable parameters (e.g. stock-recruitment models of fish population). The methods used in this paper are all "model-free" (i.e. "nonparametric") in the sense of not presuming to know the underlying mechanisms or the form of the underlying equations. The former can be criticized because the choice of a specific model may strongly bias the outcome. Modelfree or nonparametric methods have fewer assumptions, and their assumptions are usually qualitative rather than quantitative. Implemented carefully, such methods are flexible enough to avoid biases without requiring excessive amounts of data. However, modelfree and nonparametric methods are nonmechanistic, so it is difficult to incorporate any information apart from the time-series population sizes. Deriving hybrid semi-mechanistic models, which use a mechanistic "skeleton" insofar as possible, incorporating what is known about the biology generating the dynamics, and using lags and nonparametrics to deal with areas of ignorance, is an area of current research.
323
gradient r, and the population dynamics of the prey p, and predator h, are given by the following dimensionless reaction-diffusion model with logistic growth of the prey and a type II [13,14] functional response of the predator:
Op Ot
.
Oh Ot
.
--
r(x)p(1 .
.
aph 1 + bp
-
-
p)
aph + D -02p 1 + bp Ox2' OZh
m h
(2)
(3)
+ D OX2,
r(x) = 2 - 1.4x,
(4)
where the boundary conditions are zero flux at x = 0 and x = 1. The model is started with predator and prey uniformly distributed along the spatial domain. The coefficients are a = 5, b --- 5, m = 0.6, and D = 10 -4. With these parameter values, the resulting dynamics are chaotic in time [6]. The weak coupling of oscillators with different amplitudes and frequencies along the gradient leads to chaos in this model. The spatial heterogeneity of prey resource induces an additional effect on the dynamics: the temporal history of the model at the high end of the gradient is more regular than at the low end. This can be seen in Fig. 1 where the time axis shows regular cycles near x = 0 (high resource), and irregular oscillations near x = 1 (low resource).
2. Spatio-temporai population models 2.1. Continuous, spatially heterogeneous model Pascual [6] studied the dynamics of a diffusive predator-prey model in one-dimensional space. Spatial heterogeneity is introduced by letting one of the parameters, the prey's growth rate, vary linearly with space. In the absence of diffusion, the model is a standard predator-prey system, which exhibits limit cycles or equilibria [12]. The fixed (in time) spatial heterogeneity in prey resource (such as light or nutrients) induces chaos in this model. In the model, both species diffuse at the same constant rate D. At any point x and time t, the resource
1 St~ca
150
Time
Fig. l. Output from the continuous, heterogeneous, spatit-temporal predator-prey-nutrient model described in Eqs. (2)-(4) (cf. [6]). The plot shows the evolution of the prey population. Note that temporal complexity increases along the line from spatial point 0 where the nutrient level is high, to point 1 where it is low.
324
S. Little et al./Physica D 96 (1996) 321-333
The model is run for 1000 time units (after transients have died out) with the numerical scheme used in [6]. This scheme combines a fully implicit method for the diffusion terms, with a fourth-order RungeKutta method for the predator-prey interaction terms. The domain is on a 100 point spatial grid (from x = 0 to x = 1), and the model generates an output data matrix of 100 spatial points by 1000 time steps (Fig. 1). A control time series is generated by extracting all the time points from one spatial point in the matrix (a vector 1000 points long). Test time series are generated by point sampling in time the spatial axis of this matrix with random walks, and by taking time series of spatial averages at fixed locations. Each resulting test time series is a vector 1000 steps long. The coherence length of the spatial pattern is calculated from the distance at which the autocorrelation of the spatial pattern first crosses zero. The mean coherence length, an average of the coherence length of the spatial pattern at each time step, is 20 points long. The random walk's step sizes are all less than one quarter of this length, while the spatial averages start at approximately one quarter and increase to five times this length.
Five time series of randomly sampled locations are generated with random walk step sizes of 1, 2, 3, 4, and 5 spatial points. The random walk sequence is generated with a given step size, and this time series of spatial locations is used to pick points out of the model output matrix. Note that even for the smallest random walk step size of l point, the entire 100 point domain may be sampled by the time 1000 time steps have been completed. Seven time series of spatial averages are generated by averaging lengths of 6, 10, 20, 30, 50, 75 and the entire domain of 100 points. Selected examples of the time series generated by this model are shown in Fig. 2.
2.2. Discrete, spatially homogeneous model The second model is a one-dimensional coupled map lattice which describes the growth and dispersal of a single population over a finite number of habitat patches (see e.g. [15,16]). In each patch the population behaves as a map with density dependent growth, that is, the growth rate is dependent on the population size. At each time step, the population at each point is allowed to disperse over a finite spatial kernel which
Sampled Continous Model /
'Singles;atial poi'nt
.
I
Randomwalk: step size 2
.
.
.
.
.
o.~ -0.2
Random walk: step size 4
~ -0.2
o. E
Spatial average:
length 10
0.2/
Spatial average:
length aO
/
_0.2/
Spatial average:
length 100
/
0
100
200
300
400
500 600 Time step
700
800
900
1000
Fig. 2. Examples of the control and test time series generated from the continuous mode]. Note that as spatia! averaging length increases, the amplitude of the resultant time series decreases.
325
S. Little et al./Physica D 96 (1996) 321-333
Density Dependent Growth
Dispersal Kernel 6,0
0,25
0.20 4.0 0.15
O.lO 2.0 0.05
0.00
-7 -6 -5 -4 -3 -2 -I
0 I i-)
2
3
4
5
6
0.0
7
0.0
0~
1.0 n
1.5
2.0
A Snapshot of the CML 6.0 5.0 4.0 ~3.0 2.0
1.0 0.0
100
0
200
300
400
500
i Fig. 3. Pictorial description of the discrete, homogeneousspatio- temporal single species population model given in Eqs. (5) and (6). decreases geometrically with distance from the point. The population size n, in patch s, and time i, is given by the equation: 512
ni+l(s) = ~
(5)
K(s - j ) f [ n i ( j ) ] .
j=0 The dispersal kernel K ( s - j ) (graphed in Fig. 3) gives the fraction of the population in patch j which disperses to patch s: K ( s - j ) = 0.2236
therefore Z
K ( s ) = 1.
,
Is - Jl _< 5,
(6)
s
The density dependent growth in each patch is given by the Ricker curve [17] f(n): f ( n ) ----ne 4"l~l-n)
(7)
This model exhibits chaos. Unlike the first model, this model is spatially homogeneous; the parameters are independent of position. As for the first model, this model is run for 1000 time steps (after transients have died out), but with a spatial grid of 512 points, this generates a matrix of 512 spatial points by 1000 time steps. A control time series is generated by sampling at one spatial point in the matrix. Test data are generated, as above, by sampling this matrix with random walks, and with spatial averages. The mean coherence length of the spatial pattern generated by this model is about 10 points long. The random walks' step sizes are all less than this length, while the spatial averages start at approximately one quarter and increase to fifty times this length. Four time series of randomly sampled locations are generated with random walk step sizes of 1, 2, 4 and 8 spatial points. The random walk sequence is generated
326
S. Little et al./Physica D 96 (1996) 321-333
Sampled Discrete Model
<E -~'I
Spatial average: length 17
!
_21-
Spatial average: length 512
_~1
0
I
I
i
i
100
200
300
400
i
[
500 600 Time step
i
i
i
700
800
900
1000
Fig. 4. Examples of the control and test time series generated from the discrete model. Note that the amplitude of the time series decreases with increasing spatial averaging length.
with a given step size, and this time series of spatial locations is used to pick points out of the model output matrix. Note that even for the smallest random walk step size of 1 point, the entire 512 points domain may be sampled by the time 1000 time steps have been completed. Five time series of spatial averages are generated by averaging lengths of 3, 5, 9, 17 and the entire domain of 512 points. Selected examples of the time series generated by this model are shown in Fig. 4.
3. Nonlinear analysis methods 3.1. A measure of determinism To determine whether the spatial sampling schemes obscured the underlying determinism in our models, we used several techniques, three based on the concept of prediction, and one on the concept of continuity. 3.1.1. Predictability Predictability is a comparison between the output of a prediction algorithm, and "out of sample" data points, the data not used by the algorithm. The general
prescription for the prediction algorithms utilizes local neighborhoods of points in phase space and a variety of methods based on this have been published e.g. [ 1821,7]. The time evolution of points on the attractor will be smooth in local neighborhoods if no noise is present in the system. Close trajectories will stay close in local regions. This property is used for prediction in the following ways. For every point V/in phase space, find its k-nearest neighbors, where k is a tunable parameter which usually is about 1-10% of the length of the data set. These neighbors are used in various ways to generate a prediction of one time step ahead. One of the simplest prediction methods is to find the center of mass of the neighbors after one time step and compare it to the point Vi after one time step, vi+l, as described in Fig. 5. In other words, for each time i in the embedded time series V/, we calculate a predicted value at time i + l (or, more generally, i + h, where h is a "prediction horizon") in the following way. First, find the k-nearest neighbors to the embedded point V at time i. Move each neighbor ahead one time step and find the center of mass of these images of the neighbors. This gives us a " predicted value" called Yi+I,
327
S. Little et al./Physica D 96 (1996) 321-333 Prediction
TuneSeries
Prediction DelayEmbedding
~ \ \ t II centerof mass ~ \ J./F ofneatrest n e i g h b o r s \ ~ I timet+l /~ ~'x ~nearest / ,,~ ~ I I \ neighbors / ~ / I~ k yfi+l)
V(i)=(xi,xi+x)
~ n,.a,=, n~ight~s
[ t
/Oil t I II
ar map ived from arest neighbors
I /
Fig. 5. Center of mass prediction method used for Kaplan's predictability in two dimensions. A discretely sampled time series is embedded using delay coordinates. For each point V(i), a prediction is made using the center of mass of the images of its nearest neighbors one time step ahead. The difference between this predicted point, y(i + 1), and V(i + 1), is the prediction error for that point. which we compare to the actual value, Vi+l. If the geometry is smooth, the distance between Vi+l and yi+l will be small, and the prediction will be good. If the geometry is noisy, or highly chaotic, then the prediction will be worse. We evaluate the overall quality of prediction of the whole data set and get a measure of prediction error by looking at the root mean square (RMS) value of Vi+l - Y i + l for all i. This analysis was performed using an embedding dimension of 3, and time lag of 1, using 5 nearest neighbors. This method is referred to as Kaplan predictability and the results are presented as the log of the RMS value of the prediction error. In a slightly more complicated version of this [7], a linear mapping is fit in a least squares sense from the neighbors to their future values (typically one time step ahead) (see Fig. 6). This mapping then is ap-
Fig. 6. Local linear map prediction method on which Sauer's method is based. A linear map is calculated, using least squares approximation, from the nearest neighbors to point V(i), to their images one time step ahead. This map is then applied to the point V(i) to get a predicted point y(i + 1). The prediction error for that point is the difference between y(i + 1) and V(i + 1).
plied to the point Vi, and the resulting prediction yi-+-I is compared to the actual Vi+j. This method uses further refinements for noise reduction by performing local singular value decomposition (SVD) on the neighbors before calculating the mapping. In this process, the d-dimensional neighbors are projected onto a smaller (m + 1)-dimensional subdomain of the embedded space. The purpose of this is to reduce the unwanted effects of noise by keeping only the (m + 1) highest singular values. By taking more modes, or principal components, the method allows for greater flexibility in reducing noise and keeping signal. For very noisy data, few modes usually produce better predictions, while for clean data, retaining more modes produces better predictions. In this paper, all modes were calculated, but the results did not differ significantly so only m = 3 is presented. This analysis is
328
S. Little et aL/Physica D 96 (1996) 321-333
performed using embedding dimension 4, and time lag of 1, and 10-35 nearest neighbors (depending on the data set). A further variation is used, NEGM, where a weighted average is calculated of the neighbors one time step ahead, using weights proportional to the distances of the original neighbors to the point to be predicted. Yi+l : ~ Wj,i Vj+I, J
(8)
where Yi+l is the predicted value, wj,i are the weights calculated from distances between neighbors Vj to the point Vi, and ~+1 are the neighbors at one time step ahead. The advantage of this method is that nearby neighbors, which should be better predictors, are weighted more heavily. This method is called a kernel regression model [22]. Rather than comparing the prediction to linear surrogates, however, NEGM uses ordinary cross-validation [22] to find the optimum prediction while varying the parameters of embedding dimension, time delay, and the weighting function. NEGM uses this method to arrive at an optimum embedding dimension which is then used in a program to estimate the maximum Lyapunov exponent of a time series [9]. This analysis was performed with embedding dimension 1-9, depending on the data set, and time lag of 1.
3.1.2. Continuity Continuity of trajectories in phase space can be used as a test for smoothness. The idea behind this measure is that in a deterministic system, the smooth geometry is created by trajectories which stay close to each other for (at least) short periods of time. In a noisy system, the trajectories will diverge more rapidly. So smoothness is quantified by looking at how rapidly nearby points in phase space diverge. Continuity is measured as follows [8]: For a deterministic system, if two points, at time i and j , are close together in the embedding space, then their future values one time step ahead, V/+I and Vj+I should also be close together. By examining a plot of IVi+I - V j + l l (called Ei,j) VS. the distance between V/and ~ in the embed-
ding space (called ~i,j), w e can assess the amount of determinism. On a plot of ei,j v s . ~i,j, determinism appears when Ei, j tends to zero for those i, j where ~i,j is small.
3.2. Linear surrogates Surrogate data are realizations of a random process with the same histogram and similar autocorrelation (and power spectrum) as the original data. Theiler et al. [10] describe Fourier transform methods for generating surrogate data. We used "amplitude adjusted" surrogates [10, algorithm II] which tests the hypothesis that the data are a nonlinear transformation of a Gaussian linear stochastic process. We constructed ten realization of the surrogates. If the prediction error of the data is at least two standard deviations below the mean prediction error of the surrogates, we conclude that a nonlinear component is present in the underlying dynamics.
3.3. Lyapunov exponent estimation A nonmechanistic model for calculating the largest Lyapunov exponent is used in this paper: Xi+h : f ( x i , X i - r , Xi-2z . . . . .
X i - ( d - l ) r ) q- ei,
(9)
where h is the prediction horizon, ei is a noise term, and f is a statistical regression model to be fit to the data x [9]. The analyses here used the feedforward neural network model for f [23]. This method is often quite successful for modeling nonlinear and chaotic time series data e.g. [24]. In particular, the feedforward neural network is robust against errors in choosing the embedding dimension [23]. The model used explicitly allows for a random component in the data. In such systems with exogenous noise, Lyapunov exponents are calculated from the derivatives of the estimated map f , essentially by proceeding as if f were the true map but evaluating the derivatives along the observed trajectory (note that this is not the same as computing the Lyapunov exponent of the model (Eq. (9)).
329
S. Little et al./Physica D 96 (1996) 321-333
4. Model results The different methods give extremely consistent results when applied to time series from the spatiotemporal models. In the controls, they all find high predictability and detect nonlinear deterministic behavior. In the random walk sampled data, they all find that the predictability and the ability to detect underlying nonlinear determinism falls rapidly as step size increases. In the spatially averaged data they all find that predictability and the ability to detect nonlinear determinism remains high even as averaging length increases beyond the spatial coherence length. To compare consistency among techniques, the application of all the techniques to a subset of time series from only the continuous spatio-temporal model is qualitatively summarized in Table 1. The first three rows simply report the relative level of predictability measured. For each of three methods, the controls and spatial averages are highly predictable, while the random walks degrade rapidly. The next three rows report the comparison of predictability between the time series and their linear random surrogates. Again, for all methods the controls and the spatial averages are distinguishable from their linear surrogates, meaning that the underlying nonlinear deterministic model has been detected. The random
walks series, however, are indistinguishable from their linear surrogates, except for Kaplan's predictability method which is able to detect nonlinear deterministic behavior in the smaller of the two random walk step sizes. The last row, the results of NEGM Lyapunov calculations, confirm the previous results, and show positive Lyapunov exponents for the control and spatial averages, but not for the random walk samples. Secondly, to compare consistency between models, the application of the simple center of mass prediction technique to all the time series from both models is presented in Figs. 7 and 8. On both models, the control time series from a single spatial location is easily distinguishable from its linear surrogates, as evidenced by the large separation between predictability estimates for controls vs. surrogates. On both models, the random walk sequences tend towards lower predictability and smaller separation between the predictability estimates for the times series and their surrogates. The details of this trend are different between models, with separation remaining out to the largest random walk step sizes for the discrete model but not for the continuous one. On both models the prediction estimates from spatial averages are distinguishable from their linear surrogates, except for the spatial average of the entire domain for the discrete model.
Table 1 Summary of results of methods applied continuous model
Predictability, Kaplan Predictability,. Sauer Predictability, Ellner Linear surrogates, Kaplan continuity Linear surrogates, Sauer predictability Linear surrogates, Kaplan predictability Lyapunov exponent, NEGM
Single spatial point
Random walk step size 2
Random walk step size 4
Spatial average length 10
Spatial average length 30
Spatial average length 100
High High High NL/DE
Med Med Med Indist.
Low Low Low Indist.
High High High NL/DE
High High High NL/DE
High High High NL/DE
NL/DE
Indist.
Indist.
NL/DE
NL/DE
NL/DE
NL/DE
NL/DE
Indist.
NL/DE
NL/DE
NL/DE
POS
neg.
neg.
POS
POS
POS
330
S. Little et al./Physica D 96 (1996) 321-333
Kaplan Prediction on Continous Model -2.5
i
.-3
0 -
Surrogate prediction error
I I
~-3.5
ILl
SINGLE i POINT i
¢-
0
I
._o ._o -4 -lo
I I
SPATIAL AVERAGES
I ,
Q. o~-4.5
I I
RANDOM WALKS
=~;
i
..J
I I
I ,
-5
I
I
I
I
t -5.!
•I 0
i
!e I
I
[
i
1
2
3
4 5 6
Step Size
• [
i
I
L
10
20
30
50
I
I
75 100
Length
Fig. 7. Results of the application of Kaplan's center of mass prediction method to the sampled output of the continuous population model. The x-axis is a measure of distance: for the random walk it is step size, for the spatial averages it is averaging length. The prediction error for the control sample from a single point location is statistically distinguishable from that of its linear surrogates. For the random walk samples, the difference between data and their surrogates drops as step size increases, although not monotonically. Since the random walks are random samples, some may be less scrambled than others, and it is not surprising to see this kind of nonmonotonic behavior. For the spatially averaged samples, the overall prediction error is much lower for both data and surrogates, however, the difference between them remains statistically significant even out to a spatial average of the entire domain of 100 points. It is clear that random walks through these spatiotemporal systems generate time series that quickly approach stochastic noise, while spatial averages preserve the character of the dynamics. Since the methods used in this paper are searching for smoothness in phase space, it is instructive to look at phase portraits of the data. Fig. 9 presents phase portraits for the sampled model output. It is quite clear from these simple plots that (1) structure is present for spatial averages, and (2) structure is quickly lost for random walks. 5. Discussion
It is clear from the results of the model test that random walk sampling degrades the ability to detect nonlinear determinism much more rapidly than spatial averaging does. All the methods used in this paper agree with this finding.
The spatial average sample is a simple linear combination of spatial elements from the spatio-temporal system. This transformation should then, in theory, retain the dynamical information of the system. In some sense the spatial averaging is just another embedding, but in the spatial domain of a spatio-temporal system the effects of this embedding are not well understood. On the one hand, if there is noise present, this averaging will reduce noise levels and make detection of dynamics easier. On the other hand, spatial averaging also can reduce the signal level, as out-of-phase portions of space cancel each other (this can be seen in both Figs. 1 and 3, where in the noise-free case of spatial averages, the signal level drops as the spatial average length increases). In this study, an average of the entire spatial domain for the continuous model retained its underlying nonlinear deterministic structure, while an average of the entire spatial domain of the discrete model did not. The competing
S. Little et al./Physica D 96 (1996) 321-333
331
Kaplan Prediction on Discrete Model J
i
,
i
L
e.- Data prediction error
C~ Surrogate prediction error •
-0.5
•
UJ tO
._o-1.s
SINGLE POINT
SPATIAL AVERAGES RANDOM WALKS
"(:3
I:~ 0 -2 .-J -2.5
-3
i
i
i
t
L
0
1
2
4
8
Step Size
L I
I
I
35 9 17
512
Length
Fig. 8. Results of the application of Kaplan's center of mass prediction method to the sampled output of the discrete population model. The x-axis is a measure of distance: for the random walk it is step size, for the spatial averages it is averaging length. The overall predictability is lower than for the continuous model. Nevertheless, the prediction error for the control sample from a single point location is statistically distinguishable from that of its linear surrogates. For the random walk samples, the difference between data and their surrogates drops as step size increases, although not monotonically, however, it is not surprising to see this kind of nonmonotonic behavior in random samples. For the spatially averaged samples, the overall prediction error is lower for data and slightly lower for the surrogates. The difference between data and surrogates remains statistically significant for spatial averages up to 17 points. When the entire spatial domain is averaged, the time series is indistinguishable from its linear stochastic surrogates.
effect on noise and signal due to spatial averaging is an interesting area for further study in spatio-temporal systems. The random walk sampling is quite different. In the cases studied here, the individual step size is typically small compared to the correlation length of the spatial pattern. At each step, the random walk sampling adds a random component to the data whose time correlation is dependent on the correlation of spatial patterns. The derivative of the random walk sampled time series is a combination of the local time derivative and spatial derivative. This added complexity quickly degrades the ability of the techniques used here to detect the underlying deterministic dynamics. For the heterogeneous model, the random walk adds additional complexity, because over many time steps the spatial location of the sample can become quite
far from the original sample point, far enough away to be in dynamic regimes that are apparently different. For the spatially heterogeneous case, although the form of the equations for the system is the same everywhere, certain locations are more regular than others (see for example Fig. 1, where the time history at spatial point 0 is nearly periodic, while that of spatial point 1 is quite aperiodic). The random walk sample goes between these, and over the long term the resultant signal is a scrambled combination of the different regimes. So the random walk sampling can add two levels of complexity to the data, local mixing of time and space derivatives, and slow drift into different dynamic regimes. The ramifications of these results to uncovering nonlinear dynamics in the ocean are severe. They imply that the usual Eulerian sampling of the ocean may
332
S. Little et al./Physica D 96 (1996) 321-333
Continous Model
o:O4 oO 1, /D oOllal Single point sample
RW step 2
-0.2
RW step 4
-0.2
. _o,-0.4-0.2
0 0.2 0.4
-0-0.4"0.2 0 0.2 0.4
[
-0.4
-0.4-0.2 0 0.2 0.4
LO + v X
oOlO1
SA length 30
SA length 10
O.
0.
:Oo::[
-0.2 I
-0.4
. . . .
-0.4-0.2 0 0.20.4
SA length 100
-0.4-0.2 0 0.2 0.4
-0.4-0.2 0 0.2 0.4
x(t) Discrete Model RW step 2
Single point sample 4
RW step 4
,
2"
~.~!.':':"':'.';.
i:
r~-.".¢ : "o.."
~.~*.,~..:".'..
P~:.~;~,~.. ' -2 -2
0
2
-2
-2
0
2
-2 -2
0
2
+ v x SA length 17
SA length 9 4
_,
-2
f:%_ 0
2
SA length 512 4
!~£ '~:~:,.. -2
-2
0
2
-2 -2
0
2
x(t) Fig. 9. Phase portraits of selected time series from the models. The top six are from the continuous model, the bottom six are from the discrete model. The degradation of smooth geometry in phase space due to the random walk sampling is clearly illustrated in the first and third rows. The maintenance of smooth structure and the decrease in signal amplitude in the face of spatial averaging is illustrated in rows 2 and 4. be a very poor method to collect data for detecting underlying nonlinear deterministic population dynamics. A more general conclusion may also be drawn for experiments in other fields. When it is not possible to compensate for the movement of a spatio-temporal system in relation to the measurement point, it is bet-
ter to take fixed spatial averages rather than moving individual point samples. For oceanographic data, this amounts to obtaining spatial averaging of large sections of populations, for example using satellite images of plankton populations. Alternatively, drifting (Lagrangian)
S. Little et al./Physica D 96 (1996) 321-333
m e a s u r e m e n t s should be obtained, where the sample is taken from approximately the same point in the population each time step.
Acknowledgements This work arose out of a workshop on N o n l i n e a r Data A n a l y s i s in M a r i n e Ecology held at the Woods Hole O c e a n o g r a p h i c Institution in April 1994, sponsored by the Office of Naval Research. We thank the participants for discussions and ideas generated during the workshop. This work was f u n d e d in part by O N R grant n u m b e r N00014-92-J- 1527. W H O I contribution n u m b e r 9095.
References [1] M. Kot, G.S. Sayler and T.W. Schultz, Complex dynamics in a model microbial system, Bull. Math. Biol. 54 (1992) 619-648. [2] R.M. May, Simple mathematical models with very complicated dynamics, Nature 261 (1976) 459-467. [3] R.V. Sole and J. Vails, On structural stability and chaos in biological systems, J. Theoret. Biol. 155 (1992) 87-102. [4] M.P. Hassell, H.N. Comins and R.M. May, Spatial structure and chaos in insect population dynamics, Nature 353 (1991) 255-258. [5] M. Kot, Diffusion-driven period-doubling bifurcations, Biosystems 22 (1989) 279-287. [6] M. Pascual, Diffusion-induced chaos in a spatial predatorprey system, Proc. Roy. Soc. Lond. Ser. B 251 (1993) 1-7. [7] T. Sauer, Time series prediction by using delay coordinate embedding, in: Time Series Prediction: Forecasting the Future and Understanding the Past, eds. A.S. Weigend and N.A. Gershenfeld, SFI Studies in the Sciences of Complexity, Proc. Vol. XV (Addison-Wesley, Reading, MA, 1993). [8] D.T. Kaplan, Exceptional events as evidence for determinism, Physica D 73 (1994) 38-48.
333
[9] S. Ellner and E Turchin, Chaos in a "noisy" world: New methods and evidence from time series analysis, Amer. Naturalist 145 (1995) 343-375. [10] J. Theiler, J. Galdrikian, A. Longtin, S. Eubank and J.D. Farmer, Using surrogate data to detect nonlinearity in time series, in: Nonlinear Modeling and Forecasting, eds. M. Casdagli and S. Eubank (Addison-Wesley, Reading, MA, 1992). [11] D.W. Nychka, S. Ellner, D. McCaffrey and A.R. Gallant, Finding chaos in noisy systems, J. Roy. Statist. Soc. Ser. B 54 (1992) 399-426. [12] R.M. May, Stability and complexity in model ecosystems, Monographs in Population Biology (Princeton University Press, Princeton, NJ, 1973). [13] L. Michaelis and M.I. Menten, Die Kinetick der Invertinwirkung, Biochem. Z. 49 (1913) 333-369. [14] C.S. Holling, Some characteristics of simple types of predation and parasitism, Can. Ent. 91 (1959) 385-395. [15] K. Kaneko, Overview of coupled map lattices, Chaos 2(3) (1992) 279. [16] I. Waller and R. Kapral, Spatial and temporal structure in systems of coupled nonlinear oscillators, Phys. Rev. A 30 (1984) 2047-2055. [17] W.E. Ricker, Stock and recruitment, J. Fish. Res. Bd. Can. 11 (1954) 559-623. [18] J.D. Farmer and J.J. Sidorowich, Predicting chaotic time series, Phys. Rev. Lett. 59 (1987) 845-848. [19] M. Casdagli, Nonlinear prediction of chaotic time series, Physica D 35 (1989) 335-356. [20] G. Sugihara and R.M. May, Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series, Nature 344 (1990) 734-741. [21] R. Cawley and G.-H. Hsu, Local-geometric-projection method for noise reduction in chaotic maps and flows, Phys. Rev. A 46 (1992) 3057-3082. [22] B. Cheng and H. Tong, On consistent nonparametric order determination and chaos, J. Roy. Statist. Soc. Ser. B 54(2) (1992) 427-449. [23] D. McCaffrey, S. Ellner, D.W. Nychka and A.R. Gallant, Estimating the Lyapunov exponent of a chaotic system with nonlinear regression, J. Amer. Statist. Assoc. 87 (1992) 682-695. [24] M. Casdagli and S. Eubank, Nonlinear modeling and forecasting, Santa Fe Institute Studies in the Sciences of Complexity, Proc. Vol. XII (Addison-Wesley, New York, 1992).
Physica D 96 (1996) 321-333
Detecting nonlinear dynamics in spatio-temporal systems, examples from ecological models Sarah Little a,*, Stephen Ellner b, Mercedes Pascual a, Michael Neubert a, Daniel Kaplan c, Timothy Sauer d, Hal Caswell a, Andy Solow a a Woods Hole Oceanographic Institution, Wood Hole, MA 02543, USA b North Carolina State University c McGill University d George Mason University
Abstract
Mathematical models of marine populations exhibit chaotic dynamics. However, we hypothesize that in moving water, Eulerian sampling of spatially heterogeneous populations may obscure any deterministic signal beyond the resolving capabilities of presently available nonlinear signal processing techniques. To examine this hypothesis we created two spatio-temporal models of population dynamics. To caricature actual ocean sampling limitations, we sampled the model output in two ways, random walks to simulate Eulerian sampling, and spatial averages to simulate population measurements from finite volumes. Results indicate that the ability to identify underlying nonlinear dynamics quickly degrades as the step size of a random walk sampling increases. On the other hand, the analysis techniques used are more robust in the face of spatial averaging. Keywords: Nonlinear dynamics; Marine ecology; Spatio-temporal chaos; Sampling
1. Introduction
In this paper we examine how noise introduced by limitations in sampling of a chaotic spatio-temporal system affects the ability to decipher underlying nonlinear dynamics. This study arose out of work on sampling of marine ecosystems, but the results are extendable to spatio-temporal systems in general. In the ocean, complex spatio-temporal interactions between predators, prey, and nutrients, coupled with ocean physics, govern the population dynamics *Corresponding author. Fax: +1-508-289-3245; e-mail: [email protected] (508)289-3245.
of everything from phytoplankton to whales. Understanding the behavior of populations of oceanic organisms has great practical significance ranging from managing commercial fisheries such as cod and salmon, to controlling blooms of toxic plankton. However, the simple act of taking meaningful, long-term, spatio-temporal measurements of populations in a moving fluid is a formidable task itself. First, all point samples of a population must be spatial averages of some sort. Second, long-term ocean sampling is most easily accomplished in the Eulerian sense, that is, in a reference frame fixed with respect to the moving fluid, such as from a moored buoy or a pier. For an externally advected population, this effectively means
0167-2789/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved PII S0167-2789(96)00030-9
322
S. Little et al./Physica D 96 (1996) 321-333
that temporally successive samples are not taken from the same spatial location with respect to the spatial pattern of the population. The effects of these two sampling limitations on the ability to detect nonlinear deterministic behavior are addressed in this paper by applying several nonlinear time series analysis techniques to the output of two different spatio-temporal population models. In theoretical ecology, there are many examples of temporal population models which exhibit chaos. The interaction of three variables in a predatorprey-nutrient system [1] is a well studied chaotic system, as is the simple logistic map [2]. Examples of spatio-temporal models are less common, but include population models in discrete time and space [3,4], and discrete time and continuous space [5]. The tests conducted in this paper use a onedimensional continuous time and space model of a predator-prey-resource interaction [6], and a onedimensional discrete time and space model of a single species; both models exhibit chaotic dynamics. Time series representing spatial averaging are generated by taking spatial averages of constant length and location from the output of the modes. Time series which caricature Eulerian sampling are generated by sampling at points determined by a random walk in one-dimensional space. These series are compared to a time series of samples from a single spatial location. The tests used in this paper utilize several recently developed methods [7-9] which are intended for situations where the data are scarce, the functional form of the underlying process is unknown, and noise may be present. The methods assay different aspects of nonlinear dynamics - predictability, continuity, and sensitive dependence on initial conditions. The purpose of this paper is not to compare methods, but to determine the relative effects of Eulerian sampling vs. spatial averaging on our ability to detect that the underlying dynamics are deterministic spatio-temporal chaos. Because the methods are new and their performance on spatio-temporal data has not been examined, we chose several methods so that consistent answers from all methods will give us confidence in the results.
These methods, and most other nonlinear dynamics signal processing ones, are based on phase space reconstruction using delay coordinate embedding V(ti) ---- (xi, xi+r, x i + 2 r . . . . .
Xi+(d-l)r),
(1)
where V is the embedded vector, x(t) the measured variable, d the embedding dimension, r the time lag, and i the index for time t. Here the methods diverge, [8,7] tests approach the problem of detecting nonlinear determinism with a two step process. First, they develop a measure of the data which will be sensitive to determinism. Their techniques capitalize on the geometric structure of the attractor in phase space which is the hallmark of deterministic systems. On the one hand, the complexity of this geometric structure is determined by the dynamics, and on the other, noise will tend to destroy the smoothness of the geometric structures. The techniques for measuring determinism here work by looking at the smoothness of local neighborhoods in phase space of the reconstructed attractor, the idea being that the smoother the structure, the more likely that it arose from a deterministic process and not a stochastic one. Secondly, they test whether the determinism arises linearly or nonlinearly. They use the determinism measure to compare the data with linear surrogates which have the same power spectrum as the data but have a random phase spectrum [10]. The surrogates are representative of a linear random process with the same autocorrelation as the data. If a difference is found between the data and their linear surrogates, then the conclusion is drawn that the data do not arise from a linear stochastic system. Importantly though, the test does not rule out the situation where a deterministic process underlies the data but the data have a large noise component, nor does it rule out that the data arose from a nonlinear stochastic process. The Nychka, Ellner, Gallant, McCaffrey test ([9,11]) (hereafter NEGM) approaches the question of detecting nonlinear determinism by measuring the Lyapunov exponent of a model fitted to the data. Because ecological time series typically have fewer than 1000 points, it is impossible to reliably estimate Lyapunov exponents from data directly. The approach generally taken, then, is to fit a nonlinear model to
S. Little et al./Physica D 96 (1996) 321-333 the data and calculate the exponents directly from the model. The main pitfall of this approach is that the results can depend entirely on the choice of model. The issues of fitting an appropriate model to a data set are not trivial, and have been addressed in [9]. Analyses of population dynamics data have often been based on fitting the data with simple models having a few adjustable parameters (e.g. stock-recruitment models of fish population). The methods used in this paper are all "model-free" (i.e. "nonparametric") in the sense of not presuming to know the underlying mechanisms or the form of the underlying equations. The former can be criticized because the choice of a specific model may strongly bias the outcome. Modelfree or nonparametric methods have fewer assumptions, and their assumptions are usually qualitative rather than quantitative. Implemented carefully, such methods are flexible enough to avoid biases without requiring excessive amounts of data. However, modelfree and nonparametric methods are nonmechanistic, so it is difficult to incorporate any information apart from the time-series population sizes. Deriving hybrid semi-mechanistic models, which use a mechanistic "skeleton" insofar as possible, incorporating what is known about the biology generating the dynamics, and using lags and nonparametrics to deal with areas of ignorance, is an area of current research.
323
gradient r, and the population dynamics of the prey p, and predator h, are given by the following dimensionless reaction-diffusion model with logistic growth of the prey and a type II [13,14] functional response of the predator:
Op Ot
.
Oh Ot
.
--
r(x)p(1 .
.
aph 1 + bp
-
-
p)
aph + D -02p 1 + bp Ox2' OZh
m h
(2)
(3)
+ D OX2,
r(x) = 2 - 1.4x,
(4)
where the boundary conditions are zero flux at x = 0 and x = 1. The model is started with predator and prey uniformly distributed along the spatial domain. The coefficients are a = 5, b --- 5, m = 0.6, and D = 10 -4. With these parameter values, the resulting dynamics are chaotic in time [6]. The weak coupling of oscillators with different amplitudes and frequencies along the gradient leads to chaos in this model. The spatial heterogeneity of prey resource induces an additional effect on the dynamics: the temporal history of the model at the high end of the gradient is more regular than at the low end. This can be seen in Fig. 1 where the time axis shows regular cycles near x = 0 (high resource), and irregular oscillations near x = 1 (low resource).
2. Spatio-temporai population models 2.1. Continuous, spatially heterogeneous model Pascual [6] studied the dynamics of a diffusive predator-prey model in one-dimensional space. Spatial heterogeneity is introduced by letting one of the parameters, the prey's growth rate, vary linearly with space. In the absence of diffusion, the model is a standard predator-prey system, which exhibits limit cycles or equilibria [12]. The fixed (in time) spatial heterogeneity in prey resource (such as light or nutrients) induces chaos in this model. In the model, both species diffuse at the same constant rate D. At any point x and time t, the resource
1 St~ca
150
Time
Fig. l. Output from the continuous, heterogeneous, spatit-temporal predator-prey-nutrient model described in Eqs. (2)-(4) (cf. [6]). The plot shows the evolution of the prey population. Note that temporal complexity increases along the line from spatial point 0 where the nutrient level is high, to point 1 where it is low.
324
S. Little et al./Physica D 96 (1996) 321-333
The model is run for 1000 time units (after transients have died out) with the numerical scheme used in [6]. This scheme combines a fully implicit method for the diffusion terms, with a fourth-order RungeKutta method for the predator-prey interaction terms. The domain is on a 100 point spatial grid (from x = 0 to x = 1), and the model generates an output data matrix of 100 spatial points by 1000 time steps (Fig. 1). A control time series is generated by extracting all the time points from one spatial point in the matrix (a vector 1000 points long). Test time series are generated by point sampling in time the spatial axis of this matrix with random walks, and by taking time series of spatial averages at fixed locations. Each resulting test time series is a vector 1000 steps long. The coherence length of the spatial pattern is calculated from the distance at which the autocorrelation of the spatial pattern first crosses zero. The mean coherence length, an average of the coherence length of the spatial pattern at each time step, is 20 points long. The random walk's step sizes are all less than one quarter of this length, while the spatial averages start at approximately one quarter and increase to five times this length.
Five time series of randomly sampled locations are generated with random walk step sizes of 1, 2, 3, 4, and 5 spatial points. The random walk sequence is generated with a given step size, and this time series of spatial locations is used to pick points out of the model output matrix. Note that even for the smallest random walk step size of l point, the entire 100 point domain may be sampled by the time 1000 time steps have been completed. Seven time series of spatial averages are generated by averaging lengths of 6, 10, 20, 30, 50, 75 and the entire domain of 100 points. Selected examples of the time series generated by this model are shown in Fig. 2.
2.2. Discrete, spatially homogeneous model The second model is a one-dimensional coupled map lattice which describes the growth and dispersal of a single population over a finite number of habitat patches (see e.g. [15,16]). In each patch the population behaves as a map with density dependent growth, that is, the growth rate is dependent on the population size. At each time step, the population at each point is allowed to disperse over a finite spatial kernel which
Sampled Continous Model /
'Singles;atial poi'nt
.
I
Randomwalk: step size 2
.
.
.
.
.
o.~ -0.2
Random walk: step size 4
~ -0.2
o. E
Spatial average:
length 10
0.2/
Spatial average:
length aO
/
_0.2/
Spatial average:
length 100
/
0
100
200
300
400
500 600 Time step
700
800
900
1000
Fig. 2. Examples of the control and test time series generated from the continuous mode]. Note that as spatia! averaging length increases, the amplitude of the resultant time series decreases.
325
S. Little et al./Physica D 96 (1996) 321-333
Density Dependent Growth
Dispersal Kernel 6,0
0,25
0.20 4.0 0.15
O.lO 2.0 0.05
0.00
-7 -6 -5 -4 -3 -2 -I
0 I i-)
2
3
4
5
6
0.0
7
0.0
0~
1.0 n
1.5
2.0
A Snapshot of the CML 6.0 5.0 4.0 ~3.0 2.0
1.0 0.0
100
0
200
300
400
500
i Fig. 3. Pictorial description of the discrete, homogeneousspatio- temporal single species population model given in Eqs. (5) and (6). decreases geometrically with distance from the point. The population size n, in patch s, and time i, is given by the equation: 512
ni+l(s) = ~
(5)
K(s - j ) f [ n i ( j ) ] .
j=0 The dispersal kernel K ( s - j ) (graphed in Fig. 3) gives the fraction of the population in patch j which disperses to patch s: K ( s - j ) = 0.2236
therefore Z
K ( s ) = 1.
,
Is - Jl _< 5,
(6)
s
The density dependent growth in each patch is given by the Ricker curve [17] f(n): f ( n ) ----ne 4"l~l-n)
(7)
This model exhibits chaos. Unlike the first model, this model is spatially homogeneous; the parameters are independent of position. As for the first model, this model is run for 1000 time steps (after transients have died out), but with a spatial grid of 512 points, this generates a matrix of 512 spatial points by 1000 time steps. A control time series is generated by sampling at one spatial point in the matrix. Test data are generated, as above, by sampling this matrix with random walks, and with spatial averages. The mean coherence length of the spatial pattern generated by this model is about 10 points long. The random walks' step sizes are all less than this length, while the spatial averages start at approximately one quarter and increase to fifty times this length. Four time series of randomly sampled locations are generated with random walk step sizes of 1, 2, 4 and 8 spatial points. The random walk sequence is generated
326
S. Little et al./Physica D 96 (1996) 321-333
Sampled Discrete Model
<E -~'I
Spatial average: length 17
!
_21-
Spatial average: length 512
_~1
0
I
I
i
i
100
200
300
400
i
[
500 600 Time step
i
i
i
700
800
900
1000
Fig. 4. Examples of the control and test time series generated from the discrete model. Note that the amplitude of the time series decreases with increasing spatial averaging length.
with a given step size, and this time series of spatial locations is used to pick points out of the model output matrix. Note that even for the smallest random walk step size of 1 point, the entire 512 points domain may be sampled by the time 1000 time steps have been completed. Five time series of spatial averages are generated by averaging lengths of 3, 5, 9, 17 and the entire domain of 512 points. Selected examples of the time series generated by this model are shown in Fig. 4.
3. Nonlinear analysis methods 3.1. A measure of determinism To determine whether the spatial sampling schemes obscured the underlying determinism in our models, we used several techniques, three based on the concept of prediction, and one on the concept of continuity. 3.1.1. Predictability Predictability is a comparison between the output of a prediction algorithm, and "out of sample" data points, the data not used by the algorithm. The general
prescription for the prediction algorithms utilizes local neighborhoods of points in phase space and a variety of methods based on this have been published e.g. [ 1821,7]. The time evolution of points on the attractor will be smooth in local neighborhoods if no noise is present in the system. Close trajectories will stay close in local regions. This property is used for prediction in the following ways. For every point V/in phase space, find its k-nearest neighbors, where k is a tunable parameter which usually is about 1-10% of the length of the data set. These neighbors are used in various ways to generate a prediction of one time step ahead. One of the simplest prediction methods is to find the center of mass of the neighbors after one time step and compare it to the point Vi after one time step, vi+l, as described in Fig. 5. In other words, for each time i in the embedded time series V/, we calculate a predicted value at time i + l (or, more generally, i + h, where h is a "prediction horizon") in the following way. First, find the k-nearest neighbors to the embedded point V at time i. Move each neighbor ahead one time step and find the center of mass of these images of the neighbors. This gives us a " predicted value" called Yi+I,
327
S. Little et al./Physica D 96 (1996) 321-333 Prediction
TuneSeries
Prediction DelayEmbedding
~ \ \ t II centerof mass ~ \ J./F ofneatrest n e i g h b o r s \ ~ I timet+l /~ ~'x ~nearest / ,,~ ~ I I \ neighbors / ~ / I~ k yfi+l)
V(i)=(xi,xi+x)
~ n,.a,=, n~ight~s
[ t
/Oil t I II
ar map ived from arest neighbors
I /
Fig. 5. Center of mass prediction method used for Kaplan's predictability in two dimensions. A discretely sampled time series is embedded using delay coordinates. For each point V(i), a prediction is made using the center of mass of the images of its nearest neighbors one time step ahead. The difference between this predicted point, y(i + 1), and V(i + 1), is the prediction error for that point. which we compare to the actual value, Vi+l. If the geometry is smooth, the distance between Vi+l and yi+l will be small, and the prediction will be good. If the geometry is noisy, or highly chaotic, then the prediction will be worse. We evaluate the overall quality of prediction of the whole data set and get a measure of prediction error by looking at the root mean square (RMS) value of Vi+l - Y i + l for all i. This analysis was performed using an embedding dimension of 3, and time lag of 1, using 5 nearest neighbors. This method is referred to as Kaplan predictability and the results are presented as the log of the RMS value of the prediction error. In a slightly more complicated version of this [7], a linear mapping is fit in a least squares sense from the neighbors to their future values (typically one time step ahead) (see Fig. 6). This mapping then is ap-
Fig. 6. Local linear map prediction method on which Sauer's method is based. A linear map is calculated, using least squares approximation, from the nearest neighbors to point V(i), to their images one time step ahead. This map is then applied to the point V(i) to get a predicted point y(i + 1). The prediction error for that point is the difference between y(i + 1) and V(i + 1).
plied to the point Vi, and the resulting prediction yi-+-I is compared to the actual Vi+j. This method uses further refinements for noise reduction by performing local singular value decomposition (SVD) on the neighbors before calculating the mapping. In this process, the d-dimensional neighbors are projected onto a smaller (m + 1)-dimensional subdomain of the embedded space. The purpose of this is to reduce the unwanted effects of noise by keeping only the (m + 1) highest singular values. By taking more modes, or principal components, the method allows for greater flexibility in reducing noise and keeping signal. For very noisy data, few modes usually produce better predictions, while for clean data, retaining more modes produces better predictions. In this paper, all modes were calculated, but the results did not differ significantly so only m = 3 is presented. This analysis is
328
S. Little et aL/Physica D 96 (1996) 321-333
performed using embedding dimension 4, and time lag of 1, and 10-35 nearest neighbors (depending on the data set). A further variation is used, NEGM, where a weighted average is calculated of the neighbors one time step ahead, using weights proportional to the distances of the original neighbors to the point to be predicted. Yi+l : ~ Wj,i Vj+I, J
(8)
where Yi+l is the predicted value, wj,i are the weights calculated from distances between neighbors Vj to the point Vi, and ~+1 are the neighbors at one time step ahead. The advantage of this method is that nearby neighbors, which should be better predictors, are weighted more heavily. This method is called a kernel regression model [22]. Rather than comparing the prediction to linear surrogates, however, NEGM uses ordinary cross-validation [22] to find the optimum prediction while varying the parameters of embedding dimension, time delay, and the weighting function. NEGM uses this method to arrive at an optimum embedding dimension which is then used in a program to estimate the maximum Lyapunov exponent of a time series [9]. This analysis was performed with embedding dimension 1-9, depending on the data set, and time lag of 1.
3.1.2. Continuity Continuity of trajectories in phase space can be used as a test for smoothness. The idea behind this measure is that in a deterministic system, the smooth geometry is created by trajectories which stay close to each other for (at least) short periods of time. In a noisy system, the trajectories will diverge more rapidly. So smoothness is quantified by looking at how rapidly nearby points in phase space diverge. Continuity is measured as follows [8]: For a deterministic system, if two points, at time i and j , are close together in the embedding space, then their future values one time step ahead, V/+I and Vj+I should also be close together. By examining a plot of IVi+I - V j + l l (called Ei,j) VS. the distance between V/and ~ in the embed-
ding space (called ~i,j), w e can assess the amount of determinism. On a plot of ei,j v s . ~i,j, determinism appears when Ei, j tends to zero for those i, j where ~i,j is small.
3.2. Linear surrogates Surrogate data are realizations of a random process with the same histogram and similar autocorrelation (and power spectrum) as the original data. Theiler et al. [10] describe Fourier transform methods for generating surrogate data. We used "amplitude adjusted" surrogates [10, algorithm II] which tests the hypothesis that the data are a nonlinear transformation of a Gaussian linear stochastic process. We constructed ten realization of the surrogates. If the prediction error of the data is at least two standard deviations below the mean prediction error of the surrogates, we conclude that a nonlinear component is present in the underlying dynamics.
3.3. Lyapunov exponent estimation A nonmechanistic model for calculating the largest Lyapunov exponent is used in this paper: Xi+h : f ( x i , X i - r , Xi-2z . . . . .
X i - ( d - l ) r ) q- ei,
(9)
where h is the prediction horizon, ei is a noise term, and f is a statistical regression model to be fit to the data x [9]. The analyses here used the feedforward neural network model for f [23]. This method is often quite successful for modeling nonlinear and chaotic time series data e.g. [24]. In particular, the feedforward neural network is robust against errors in choosing the embedding dimension [23]. The model used explicitly allows for a random component in the data. In such systems with exogenous noise, Lyapunov exponents are calculated from the derivatives of the estimated map f , essentially by proceeding as if f were the true map but evaluating the derivatives along the observed trajectory (note that this is not the same as computing the Lyapunov exponent of the model (Eq. (9)).
329
S. Little et al./Physica D 96 (1996) 321-333
4. Model results The different methods give extremely consistent results when applied to time series from the spatiotemporal models. In the controls, they all find high predictability and detect nonlinear deterministic behavior. In the random walk sampled data, they all find that the predictability and the ability to detect underlying nonlinear determinism falls rapidly as step size increases. In the spatially averaged data they all find that predictability and the ability to detect nonlinear determinism remains high even as averaging length increases beyond the spatial coherence length. To compare consistency among techniques, the application of all the techniques to a subset of time series from only the continuous spatio-temporal model is qualitatively summarized in Table 1. The first three rows simply report the relative level of predictability measured. For each of three methods, the controls and spatial averages are highly predictable, while the random walks degrade rapidly. The next three rows report the comparison of predictability between the time series and their linear random surrogates. Again, for all methods the controls and the spatial averages are distinguishable from their linear surrogates, meaning that the underlying nonlinear deterministic model has been detected. The random
walks series, however, are indistinguishable from their linear surrogates, except for Kaplan's predictability method which is able to detect nonlinear deterministic behavior in the smaller of the two random walk step sizes. The last row, the results of NEGM Lyapunov calculations, confirm the previous results, and show positive Lyapunov exponents for the control and spatial averages, but not for the random walk samples. Secondly, to compare consistency between models, the application of the simple center of mass prediction technique to all the time series from both models is presented in Figs. 7 and 8. On both models, the control time series from a single spatial location is easily distinguishable from its linear surrogates, as evidenced by the large separation between predictability estimates for controls vs. surrogates. On both models, the random walk sequences tend towards lower predictability and smaller separation between the predictability estimates for the times series and their surrogates. The details of this trend are different between models, with separation remaining out to the largest random walk step sizes for the discrete model but not for the continuous one. On both models the prediction estimates from spatial averages are distinguishable from their linear surrogates, except for the spatial average of the entire domain for the discrete model.
Table 1 Summary of results of methods applied continuous model
Predictability, Kaplan Predictability,. Sauer Predictability, Ellner Linear surrogates, Kaplan continuity Linear surrogates, Sauer predictability Linear surrogates, Kaplan predictability Lyapunov exponent, NEGM
Single spatial point
Random walk step size 2
Random walk step size 4
Spatial average length 10
Spatial average length 30
Spatial average length 100
High High High NL/DE
Med Med Med Indist.
Low Low Low Indist.
High High High NL/DE
High High High NL/DE
High High High NL/DE
NL/DE
Indist.
Indist.
NL/DE
NL/DE
NL/DE
NL/DE
NL/DE
Indist.
NL/DE
NL/DE
NL/DE
POS
neg.
neg.
POS
POS
POS
330
S. Little et al./Physica D 96 (1996) 321-333
Kaplan Prediction on Continous Model -2.5
i
.-3
0 -
Surrogate prediction error
I I
~-3.5
ILl
SINGLE i POINT i
¢-
0
I
._o ._o -4 -lo
I I
SPATIAL AVERAGES
I ,
Q. o~-4.5
I I
RANDOM WALKS
=~;
i
..J
I I
I ,
-5
I
I
I
I
t -5.!
•I 0
i
!e I
I
[
i
1
2
3
4 5 6
Step Size
• [
i
I
L
10
20
30
50
I
I
75 100
Length
Fig. 7. Results of the application of Kaplan's center of mass prediction method to the sampled output of the continuous population model. The x-axis is a measure of distance: for the random walk it is step size, for the spatial averages it is averaging length. The prediction error for the control sample from a single point location is statistically distinguishable from that of its linear surrogates. For the random walk samples, the difference between data and their surrogates drops as step size increases, although not monotonically. Since the random walks are random samples, some may be less scrambled than others, and it is not surprising to see this kind of nonmonotonic behavior. For the spatially averaged samples, the overall prediction error is much lower for both data and surrogates, however, the difference between them remains statistically significant even out to a spatial average of the entire domain of 100 points. It is clear that random walks through these spatiotemporal systems generate time series that quickly approach stochastic noise, while spatial averages preserve the character of the dynamics. Since the methods used in this paper are searching for smoothness in phase space, it is instructive to look at phase portraits of the data. Fig. 9 presents phase portraits for the sampled model output. It is quite clear from these simple plots that (1) structure is present for spatial averages, and (2) structure is quickly lost for random walks. 5. Discussion
It is clear from the results of the model test that random walk sampling degrades the ability to detect nonlinear determinism much more rapidly than spatial averaging does. All the methods used in this paper agree with this finding.
The spatial average sample is a simple linear combination of spatial elements from the spatio-temporal system. This transformation should then, in theory, retain the dynamical information of the system. In some sense the spatial averaging is just another embedding, but in the spatial domain of a spatio-temporal system the effects of this embedding are not well understood. On the one hand, if there is noise present, this averaging will reduce noise levels and make detection of dynamics easier. On the other hand, spatial averaging also can reduce the signal level, as out-of-phase portions of space cancel each other (this can be seen in both Figs. 1 and 3, where in the noise-free case of spatial averages, the signal level drops as the spatial average length increases). In this study, an average of the entire spatial domain for the continuous model retained its underlying nonlinear deterministic structure, while an average of the entire spatial domain of the discrete model did not. The competing
S. Little et al./Physica D 96 (1996) 321-333
331
Kaplan Prediction on Discrete Model J
i
,
i
L
e.- Data prediction error
C~ Surrogate prediction error •
-0.5
•
UJ tO
._o-1.s
SINGLE POINT
SPATIAL AVERAGES RANDOM WALKS
"(:3
I:~ 0 -2 .-J -2.5
-3
i
i
i
t
L
0
1
2
4
8
Step Size
L I
I
I
35 9 17
512
Length
Fig. 8. Results of the application of Kaplan's center of mass prediction method to the sampled output of the discrete population model. The x-axis is a measure of distance: for the random walk it is step size, for the spatial averages it is averaging length. The overall predictability is lower than for the continuous model. Nevertheless, the prediction error for the control sample from a single point location is statistically distinguishable from that of its linear surrogates. For the random walk samples, the difference between data and their surrogates drops as step size increases, although not monotonically, however, it is not surprising to see this kind of nonmonotonic behavior in random samples. For the spatially averaged samples, the overall prediction error is lower for data and slightly lower for the surrogates. The difference between data and surrogates remains statistically significant for spatial averages up to 17 points. When the entire spatial domain is averaged, the time series is indistinguishable from its linear stochastic surrogates.
effect on noise and signal due to spatial averaging is an interesting area for further study in spatio-temporal systems. The random walk sampling is quite different. In the cases studied here, the individual step size is typically small compared to the correlation length of the spatial pattern. At each step, the random walk sampling adds a random component to the data whose time correlation is dependent on the correlation of spatial patterns. The derivative of the random walk sampled time series is a combination of the local time derivative and spatial derivative. This added complexity quickly degrades the ability of the techniques used here to detect the underlying deterministic dynamics. For the heterogeneous model, the random walk adds additional complexity, because over many time steps the spatial location of the sample can become quite
far from the original sample point, far enough away to be in dynamic regimes that are apparently different. For the spatially heterogeneous case, although the form of the equations for the system is the same everywhere, certain locations are more regular than others (see for example Fig. 1, where the time history at spatial point 0 is nearly periodic, while that of spatial point 1 is quite aperiodic). The random walk sample goes between these, and over the long term the resultant signal is a scrambled combination of the different regimes. So the random walk sampling can add two levels of complexity to the data, local mixing of time and space derivatives, and slow drift into different dynamic regimes. The ramifications of these results to uncovering nonlinear dynamics in the ocean are severe. They imply that the usual Eulerian sampling of the ocean may
332
S. Little et al./Physica D 96 (1996) 321-333
Continous Model
o:O4 oO 1, /D oOllal Single point sample
RW step 2
-0.2
RW step 4
-0.2
. _o,-0.4-0.2
0 0.2 0.4
-0-0.4"0.2 0 0.2 0.4
[
-0.4
-0.4-0.2 0 0.2 0.4
LO + v X
oOlO1
SA length 30
SA length 10
O.
0.
:Oo::[
-0.2 I
-0.4
. . . .
-0.4-0.2 0 0.20.4
SA length 100
-0.4-0.2 0 0.2 0.4
-0.4-0.2 0 0.2 0.4
x(t) Discrete Model RW step 2
Single point sample 4
RW step 4
,
2"
~.~!.':':"':'.';.
i:
r~-.".¢ : "o.."
~.~*.,~..:".'..
P~:.~;~,~.. ' -2 -2
0
2
-2
-2
0
2
-2 -2
0
2
+ v x SA length 17
SA length 9 4
_,
-2
f:%_ 0
2
SA length 512 4
!~£ '~:~:,.. -2
-2
0
2
-2 -2
0
2
x(t) Fig. 9. Phase portraits of selected time series from the models. The top six are from the continuous model, the bottom six are from the discrete model. The degradation of smooth geometry in phase space due to the random walk sampling is clearly illustrated in the first and third rows. The maintenance of smooth structure and the decrease in signal amplitude in the face of spatial averaging is illustrated in rows 2 and 4. be a very poor method to collect data for detecting underlying nonlinear deterministic population dynamics. A more general conclusion may also be drawn for experiments in other fields. When it is not possible to compensate for the movement of a spatio-temporal system in relation to the measurement point, it is bet-
ter to take fixed spatial averages rather than moving individual point samples. For oceanographic data, this amounts to obtaining spatial averaging of large sections of populations, for example using satellite images of plankton populations. Alternatively, drifting (Lagrangian)
S. Little et al./Physica D 96 (1996) 321-333
m e a s u r e m e n t s should be obtained, where the sample is taken from approximately the same point in the population each time step.
Acknowledgements This work arose out of a workshop on N o n l i n e a r Data A n a l y s i s in M a r i n e Ecology held at the Woods Hole O c e a n o g r a p h i c Institution in April 1994, sponsored by the Office of Naval Research. We thank the participants for discussions and ideas generated during the workshop. This work was f u n d e d in part by O N R grant n u m b e r N00014-92-J- 1527. W H O I contribution n u m b e r 9095.
References [1] M. Kot, G.S. Sayler and T.W. Schultz, Complex dynamics in a model microbial system, Bull. Math. Biol. 54 (1992) 619-648. [2] R.M. May, Simple mathematical models with very complicated dynamics, Nature 261 (1976) 459-467. [3] R.V. Sole and J. Vails, On structural stability and chaos in biological systems, J. Theoret. Biol. 155 (1992) 87-102. [4] M.P. Hassell, H.N. Comins and R.M. May, Spatial structure and chaos in insect population dynamics, Nature 353 (1991) 255-258. [5] M. Kot, Diffusion-driven period-doubling bifurcations, Biosystems 22 (1989) 279-287. [6] M. Pascual, Diffusion-induced chaos in a spatial predatorprey system, Proc. Roy. Soc. Lond. Ser. B 251 (1993) 1-7. [7] T. Sauer, Time series prediction by using delay coordinate embedding, in: Time Series Prediction: Forecasting the Future and Understanding the Past, eds. A.S. Weigend and N.A. Gershenfeld, SFI Studies in the Sciences of Complexity, Proc. Vol. XV (Addison-Wesley, Reading, MA, 1993). [8] D.T. Kaplan, Exceptional events as evidence for determinism, Physica D 73 (1994) 38-48.
333
[9] S. Ellner and E Turchin, Chaos in a "noisy" world: New methods and evidence from time series analysis, Amer. Naturalist 145 (1995) 343-375. [10] J. Theiler, J. Galdrikian, A. Longtin, S. Eubank and J.D. Farmer, Using surrogate data to detect nonlinearity in time series, in: Nonlinear Modeling and Forecasting, eds. M. Casdagli and S. Eubank (Addison-Wesley, Reading, MA, 1992). [11] D.W. Nychka, S. Ellner, D. McCaffrey and A.R. Gallant, Finding chaos in noisy systems, J. Roy. Statist. Soc. Ser. B 54 (1992) 399-426. [12] R.M. May, Stability and complexity in model ecosystems, Monographs in Population Biology (Princeton University Press, Princeton, NJ, 1973). [13] L. Michaelis and M.I. Menten, Die Kinetick der Invertinwirkung, Biochem. Z. 49 (1913) 333-369. [14] C.S. Holling, Some characteristics of simple types of predation and parasitism, Can. Ent. 91 (1959) 385-395. [15] K. Kaneko, Overview of coupled map lattices, Chaos 2(3) (1992) 279. [16] I. Waller and R. Kapral, Spatial and temporal structure in systems of coupled nonlinear oscillators, Phys. Rev. A 30 (1984) 2047-2055. [17] W.E. Ricker, Stock and recruitment, J. Fish. Res. Bd. Can. 11 (1954) 559-623. [18] J.D. Farmer and J.J. Sidorowich, Predicting chaotic time series, Phys. Rev. Lett. 59 (1987) 845-848. [19] M. Casdagli, Nonlinear prediction of chaotic time series, Physica D 35 (1989) 335-356. [20] G. Sugihara and R.M. May, Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series, Nature 344 (1990) 734-741. [21] R. Cawley and G.-H. Hsu, Local-geometric-projection method for noise reduction in chaotic maps and flows, Phys. Rev. A 46 (1992) 3057-3082. [22] B. Cheng and H. Tong, On consistent nonparametric order determination and chaos, J. Roy. Statist. Soc. Ser. B 54(2) (1992) 427-449. [23] D. McCaffrey, S. Ellner, D.W. Nychka and A.R. Gallant, Estimating the Lyapunov exponent of a chaotic system with nonlinear regression, J. Amer. Statist. Assoc. 87 (1992) 682-695. [24] M. Casdagli and S. Eubank, Nonlinear modeling and forecasting, Santa Fe Institute Studies in the Sciences of Complexity, Proc. Vol. XII (Addison-Wesley, New York, 1992).
