New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Locating the transition from periodic oscillations to spatiotemporal chaos in the wake of invasion

Edited by John Guckenheimer, Cornell University, Ithaca, NY, and accepted by the Editorial Board April 22, 2009 (received for review January 7, 2009)
Abstract
In systems with cyclic dynamics, invasions often generate periodic spatiotemporal oscillations, which undergo a subsequent transition to chaos. The periodic oscillations have the form of a wavetrain and occur in a band of constant width. In applications, a key question is whether one expects spatiotemporal data to be dominated by regular or irregular oscillations or to involve a significant proportion of both. This depends on the width of the wavetrain band. Here, we present mathematical theory that enables the direct calculation of this width. Our method synthesizes recent developments in stability theory and computation. It is developed for only 1 equation system, but because this is a normal form close to a Hopf bifurcation, the results can be applied directly to a wide range of models. We illustrate this by considering a classic example from ecology: wavetrains in the wake of the invasion of a prey population by predators.
 periodic travelling wave
 absolute stability
 reactiondiffusion
 complex Ginzburg–Landau equation
 predator–prey
Wavetrains are a fundamental solution type in spatially extended oscillatory systems. They were observed in the Belousov–Zhabotinskii reaction >30 years ago (1), and experimental and theoretical demonstrations of their importance are now widespread, in applications including hydrodynamics (2–4), solar cycles (5), chemical reactions (6, 7), cell biology (8), and ecology (9, 10). Moreover, wavetrains provide the background state for many more complicated behaviors, including spatiotemporal chaos. A striking illustration of this is provided by the invasion of an unstable equilibrium, where a band of wavetrains can occur between the invasion front and a region of spatiotemporal chaos (Fig. 1 A). This behavior is associated with the wavetrain being unstable. Intuitively, the wavetrain is generated by the invasion process, but at the same time, unstable linear modes are excited. Immediately behind the front, the amplitude of these modes is negligible, and the wavetrain is visible; further back, the instabilities have grown larger, eventually dominating the solution at some distance behind the invasion front. Behavior of this type occurs in a range of model types including reaction–diffusion equations (11–13), the complex Ginzburg–Landau equation and generalizations (14, 15), integrodifferential equations (16), integrodifference equations (17), and cellular automata (18). The beststudied application is the invasion of a prey population by predators (9, 19). Here, the wavetrain would correspond to spatially organised multiyear population cycles, which have been observed in a number of field studies (10, 19). Simulations such as Fig. 1 B show that after an initial increase, the width of the wavetrain band remains constant as the invasion progresses. The key issue of whether wavetrains will be observed in applications depends on the value of the band width, which we calculate. We focus on 1 particular equation system, but its generic nature enables our results to be applied directly to a wide range of models spanning many application areas.
Wavetrains Behind Invasion in a Prototype Model
We develop our method for equations of “λ – ω” reaction–diffusion form:
This is the complex Ginzburg–Landau equation with a real diffusion coefficient (20). Its kinetics are the normal form of any oscillatory system of differential equations close to a standard supercritical Hopf bifurcation, making our calculations directly relevant to many models with lowamplitude cycles. Eqs. 1 are most conveniently studied in terms of amplitude r = (u ^{2} + v ^{2})^{1/2} and phase θ = tan^{−1}(v/u); wavetrain solutions have r = R, a constant, and θ = (ω_{0} − ω_{1} R ^{2})t ± (1 − R ^{2})^{1/2} x, and exist for 0 ≤ R ≤ 1. The scenario we consider is local perturbation (near the lefthand boundary) of u = v = 0, on a large finite domain with zero Neumann boundary conditions. The resulting dynamics depend only on the parameter ω_{1}; ω_{0} simply determines the rotation rate of the phase, and does not affect r or θ_{x}. For suitable values of ω_{1} (detailed below), the solution has the form of an invasive transition wave in r, with spatiotemporal irregularities further back (Fig. 1 B). Careful numerical studies suggest that these irregular oscillations are a genuine example of spatiotemporal chaos (21). Substitution of the travelling wave ansatz
Stability in Moving Frames of Reference
The key to understanding the width of the wavetrain band is to consider absolute stability (26) when viewed in a frame of reference moving with a fixed, arbitrary velocity V. That is, we consider whether perturbations to the wavetrain grow or decay over time when viewed at a fixed point traveling with velocity V. Perturbations that grow but simultaneously propagate at other velocities may contribute to convective instability in this frame of reference, but these are not relevant to our calculation. We denote by λ_{max}(V) the growth rate λ of the most unstable linear mode, with ν_{max}(V) being the corresponding spatial eigenvalue. Using general theory of absolute stability (26, 27), we reduce the calculation of these quantities to the numerical tracking of solutions of a quartic polynomial as parameters in the coefficients vary (see Methods). Fig. 2 A illustrates the typical form of λ_{max}(V). There is a range of velocities (V _{L},V _{R}) for which λ_{max} > 0; all perturbations decay in frames of reference moving with velocities below V _{L} or above V _{R}. Note that in this figure, V _{R} < c _{inv}, so that all growing perturbations move away from the invasion front, resulting in a nonzero band width. One advantage of our method of calculation is that it is straightforward to monitor V _{L} and V _{R} while varying ω_{1} (Fig. 2 B). V _{L} = V _{R} at the onset of wave instability, which can be shown analytically to occur at ω_{1} = ω_{1} ^{stab} = 1.0714 (for c _{inv} = 2) (23, 28). As ω_{1} increases through ω_{1} ^{stab}, the behavior behind invasion changes from an uninterrupted expanse of wavetrain, to spatiotemporal irregularity behind a wavetrain band. The change in sign of V _{R} corresponds to the onset of absolute instability in a stationary frame of reference; we compute this point as ω_{1} = 1.512658 (for c _{inv} = 2).
A Formula for the Band Width
Having calculated λ_{max}(V), we can address the width of the wavetrain band. For this, we require a precise definition of its lefthand edge, which we take as the first point at which the perturbations to the wavetrain that are present immediately behind the invasion front become amplified by a factor ℱ. The value of ℱ is arbitrary, but we will show that the dependence of the wavetrain band width on ℱ and ω_{1} decouples, having the form log(ℱ)𝒲(ω_{1}). We refer to 𝒲(ω_{1}) as the “band width coefficient.” Let (x*,t*) be a point on the invasion front. We make the generic assumption that the most unstable linear modes are present in the perturbation given to the wavetrain at such a point. As t increases above t*, these perturbations will spread out in space and time, growing along all rays x = x* + (t − t*)V with V ∈ (V _{L},V _{R}) (a selection of such rays is sketched in Fig. 3). On any such ray, the initial perturbation becomes amplified by the factor ℱ at time t _{crit}(V) = t* + log(ℱ)/Re[λ_{max}(V)], at location x _{crit}(V) = x* + V log(ℱ)/Re[λ_{max}(V)]. The curve (x _{crit}(V),t _{crit}(V)) (V _{L} < V < V _{R}) is indicated by the thick solid line in Fig. 3. The lefthand edge of the wavetrain band occurs at the point on this curve that is closest to the invasion front x = x* + (t − t*)c _{inv}, which occurs at V = V _{band}, where The width itself is given by We show in Methods that λ_{max}′(V) = ν_{max}(V), so that Eq. 2 is with the band width coefficient 𝒲 = 1/Re[ν_{max}(V _{band})].
It is of course entirely expected that 𝒲 is the reciprocal of the real part of the spatial eigenvalue; our key result is Eq. 3 for the frame velocity V _{band} at which this eigenvalue should be calculated. The form of V _{band}(ω_{1}) is plotted in Fig. 2 B, and the corresponding width coefficient 𝒲 is shown in Fig. 4 A for a limited range of ω_{1}. At the onset of instability (ω_{1} = 1.0714), V _{band} = V _{L} = V _{R} and 𝒲 = ∞. The width decreases with ω_{1}, through the onset of absolute stability (ω_{1} = 1.512658). This transition has no particular implications for the width of the wavetrain band, though it does alter the behavior behind the band, from a pattern of sources and sinks to more comprehensive spatiotemporal disorder (29, 30). Finally, at ω_{1} = ω_{1} ^{max} = 18.72751, V _{band} reaches the invasion speed c _{inv} = 2. At this point, V _{R} = 2 also, and the V _{band}(ω_{1}) curve folds; there is another branch of solutions of Eq. 3 with values greater than V _{R}, which has no practical relevance. For ω_{1} > ω_{1} ^{max}, perturbations to the wavetrain can actually outrun the invasion, so that there is no wavetrain band.
To test our predictions, we performed numerical simulations of Eqs. 1, defining the wavetrain band by the condition that ∂r/∂x is below a small, arbitrary threshold. The comparison with our theoretical prediction 𝒲 is extremely good (Fig. 4). Moreover, the slope of the regression line in this figure is a suitable value for log(ℱ), which can be used to convert predictions of band width coefficient to actual band widths in applications. As a further test, we repeated this comparison using initial conditions in which u and v decay slowly to zero, rather than being exactly zero away from the lefthand boundary. This can generate a “pushed” rather than “pulled” invasive wave (4) with speed c _{inv} > 2, which changes both the wavetrain amplitude and the band width. Fig. 4 shows that for c _{inv} = 4 the comparison between theory and simulations is again extremely good.
Example of Application: Predator–Prey Invasion
The wider significance of our results hinges on the fact that the kinetics of the λ – ω system (Eqs. 1) are not only the standard prototype oscillator but also the normal form of any oscillatory system near a standard supercritical Hopf bifurcation. This makes our results applicable to a wide range of models with lowamplitude cycles across many disciplines. To illustrate this, we consider a classic ecological example: the invasion of a prey population by predators, using the standard Rosenzweig–MacArthur model (see Methods). We do not attempt a detailed ecological investigation—rather, our objective is to illustrate the ease with which our method can be applied. To simulate invasion, we use initial conditions in which a small, localized predator density is introduced into a prey population at its carrying capacity. This generates an advancing front of predators, and a corresponding receding front of prey. A large body of theoretical work has shown that in the case of cyclic population dynamics, the behavior behind the invasion front consists of spatiotemporal oscillations (13, 31, 32). In some cases, these are a stable wavetrain; when the appropriate wavetrain is unstable, more complicated dynamics occur, with a typical example shown in Fig. 5. Immediately behind the invasion front, the population densities are almost constant, at their (unstable) coexistence steady state. Behind this, there is a wavetrain band whose width remains constant as the invasion progresses, and which is followed by spatiotemporal chaos (31, 32).
The predator–prey invasion process is different from that studied for the λ – ω equations; in particular, the preyonly state has no analogue in Eqs. 1, with u = v = 0 corresponding to a constant coexistence of predators and prey. However, detailed calculations show that the behavior behind the initial invasion front corresponds exactly to the λ – ω invasion scenario (33). This is a general result, not restricted to the particular case of Eqs. 4; effectively, the leading invasion front leaves the system at the coexistence steady state and also initiates an invasion of that state in the opposite direction. This means that close to Hopf bifurcation in the kinetics, the width of the wavetrain band is given by our calculations for the λ – ω system (Eqs. 1).
The key step in applying our results, both to predator–prey invasion and more generally, is thus to determine ω_{1} as a function of the model parameters. This is done by the standard method of reduction to normal form (34); detailed presentations of this calculation for the Rosenzweig–MacArthur equations are given in refs. 35 and 36. By combining this with our previous calculation of the band width coefficient as a function of ω_{1}, we can calculate the dependency of the width of the wavetrain band on the underlying ecological parameters. For example, Fig. 6 shows contours of band width (with exponential spacing) in a key parameter plane. Note that because they are based on a normal form reduction, predictions such as those in Fig. 6 are valid only as approximations close to Hopf bifurcation in the local dynamics.
One natural population in which wavetrains have been widely reported is the field vole (Microtus agrestis) (10, 37, 38). Possible causes of the cyclic dynamics of this population include their predation by weasels (Mustela nivalis) (39, 40). The Rosenzweig–MacArthur model omits a number of factors that are likely to influence the behavior of any real weasel/field vole interaction, but nevertheless it represents a useful caricature model, and typical parameter estimates (35) yield a band width coefficient of 156. It is convenient to express this as the number of wavelengths in the wavetrain band, which is given by multiplying by our estimate for log(ℱ) and dividing by the wavelength = 2π/(1 − R*)^{1/2}. This predicts 295 wavelengths in the band, which is much longer than any practical ecological habitat, suggesting that even though the appropriate wavetrain is unstable, spatiotemporal disorder will not be observed after invasion. This conclusion is valid even after the whole domain has been invaded, for typical (zero flux) boundary conditions, because wavetrain solutions then undergo a gradual transition to spatially uniform oscillations. In contrast, systems in which the prey birth rate is smaller relative to that of the predator would have a smaller band width, or in extreme cases, no wavetrain band at all (Fig. 6). Irregular oscillations would then develop more quickly after invasion, and would persist after the whole domain had been invaded (25). This illustrates one important practical benefit of band width calculations for conducting experiments and surveys: estimating whether or not the domain is sufficiently large for irregular spatiotemporal dynamics to be observed. It is important to note that our results only apply to the behavior after the invasion of prey by predators. Other scenarios may generate different behaviors, including spatiotemporal disorder, in the model Eqs. 4—for instance, landscape obstacles can also generate chaos via unstable wavetrains (41). Moreover, some alternative predator–prey models predict an even wider range of spatiotemporal phenomena, including “patchy invasion” in which the invasion process and the transition to chaos are simultaneous (9, 42). For the specific case of field voles, field data reveal a variety of population dynamics, including both regular cycles and quasichaos (40). A natural extension to our study is to apply the methods to more realistic models, potentially including the stochasticity that is implicit in both the underlying biology and the process of data collection.
An important implication of the contours in Fig. 6 is that the width of the wavetrain band is much more sensitive to variations in some parameters than others. For instance, in the field vole–weasel case, a 5% increase in predator mortality reduces the band width by just 0.5%, whereas a corresponding increase in prey birth rate increases the band width by 22%. An even more striking example of this nonlinearity is provided by the interaction between the zooplankton Daphnia pulex and the phytoplankton Chlamydomonas reinhardii. To the best of our knowledge, there are no spatiotemporal data from which wavetrains could be inferred in this case. However, using estimates of the relevant parameters (43), we predict that the width of a wavetrain band behind invasion would be extremely sensitive to the estimate of zooplankton birth rate, with a reduction of only 5.2% sufficient to double the band width. Such sensitivity has an important general implication: Major differences in spatiotemporal dynamics may result from even slight changes in parameter values. For example, a large body of evidence now indicates that as well as increasing the frequency of ecological invasions (44), climate change is having a significant effect on the demographic parameters underlying oscillatory ecological systems (45, 46). Our results suggest that the consequent changes in some spatiotemporal behaviors may be even more dramatic.
Methods
Calculation of λ_{max}(V).
The key quantity in our calculation of wavetrain band width is λ_{max}(V), the maximum growth rate of a linear mode in a frame of reference moving with velocity V. To determine this, we replace x by the moving coordinate x − Vt, and then compute “branch points” in the absolute spectrum (26). We do this via the numerical continuation approach of Rademacher et al. (27). Briefly, we linearize the equations satisfied by r(x,t) and θ(x,t) about the wavetrain, and denote by 𝒟(λ,ν;V) = 0 the dispersion relation for linear modes with temporal and spatial eigenvalues λ and ν, respectively. For given λ, 𝒟 is a quartic polynomial in ν, and we label the roots ν_{1},…,ν_{4} so that Reν_{i} ≥ Reν_{i+1}. The absolute spectrum defined in ref. 26 is the set of λ for which Reν_{2} = Reν_{3}, and the branch points are the 4 values of λ for which 𝒟 has repeated roots for ν. Note that for a branch point, the condition ν_{2} = ν_{3} corresponds to the “pinching condition” of refs. 47, 48, and its instability implies pointwise growth of perturbations in the comoving frame. In the language of ref. 4, absolute instability corresponds to a positive value of the linear spreading speed.
More generally, the absolute spectrum is the accumulation set of eigenvalues in the limit of large domain length under generic separated boundary conditions. Branch points typically form endpoints of curves of the absolute spectrum, which can therefore contain λs with a larger real part than any branch point in the absolute spectrum. However, these would lead to “remnant instabilities,” associated with perturbations reflecting repeatedly off the domain boundaries (26). Hence, we expect that such an instability would be irrelevant in the present context. In fact, computations of the full absolute spectrum (via continuation from branch points, as in ref. 27), for a wide range of ω_{1} and V, indicate that for our system, the most unstable point in the absolute spectrum is always a branch point.
Therefore λ_{max}(V) can be found simply by calculating the four branch points, which are the solutions of 𝒟(λ,ν;V) = 𝒟 _{ν}(λ,ν;V) = 0, and identifying which of these are in the absolute spectrum. We numerically continue these solutions in V, giving λ_{max} and the corresponding spatial eigenvalue ν_{max} as functions of V. A complication is that the appropriate branch point can change at “triple points,” when 3 roots of 𝒟(λ,ν) have equal real parts; therefore, and to check Reν_{2} = Reν_{3}, it is necessary to continue all 4 branch points and monitor their ν values. We calculate V _{band} simply by monitoring the 2 sides of Eq. 3 during the continuation. Moreover, having found V _{band} for one value of ω_{1}, continuation in ω_{1} is straightforward, enabling the plot in Fig. 2 B.
The above procedure is significantly simpler, and more portable, than that used in the few previous articles we are aware of that compute λ_{max}(V) for other systems (49–51); these studies involve continuation of a saddle point of λ(ν) − Vν in the complex ν plane.
Calculation of λ_{max}′(V).
The derivative λ_{max}′(V) is required for calculation of the band width, and this can be found very easily. We have 𝒟(λ_{max},ν_{max};V) = 0 for all V, so that at λ = λ_{max} and ν = ν_{max} (subscripts denote partial derivatives). But since λ_{max} and ν_{max} occur at a branch point, 𝒟 _{ν} = 0. Therefore λ_{max}′(V) = −𝒟 _{V}/𝒟 _{λ}. For any reaction–diffusion system, the dispersion relation satisfies Therefore λ_{max}′(V) = ν_{max}(V),
Predator–Prey Example.
The application of our results to predator–prey invasion uses the Rosenzweig–MacArthur model (52–54). Here p and h denote predator and prey densities, which depend on space x and time t. The populations are assumed to have the same dispersal coefficient, δ; however, we have shown previously that the “band width phenomenon” also occurs in this model when the dispersal coefficients are different (41) or even densitydependent (55). The local population dynamics depend on the positive parameters a, b, c, r, h _{0}, and k. The prey consumption rate per predator is an increasing saturating function of the prey density: c represents the maximal consumption rate, and k reflects how quickly the consumption rate saturates as prey density increases. Parameters a and r denote maximal percapita predator and prey birth rates; for predators, this is the birth rate when the prey density is very high, whereas for prey it is the birth rate at very low prey density. The percapita predator death rate is denoted by b, and h _{0} is the preycarrying capacity. For suitable parameter values, the local dynamics of Eqs. 4 are oscillatory: specifically, there is a standard supercritical Hopf bifurcation as kh _{0}(a − b) − a − b increases through zero.
Acknowledgments
We thank Björn Sandstede and Leonid Brevdo for helpful discussions and Andrew White and Kevin Painter for comments on the manuscript. J.A.S. was supported in part by a Leverhulme Trust Research Fellowship. J.D.M.R. acknowledges support from the Nonlinear Dynamics of Natural Systems program of the Nederlandse Organisatie voor Wetenschappelijk Onderzoek.
Footnotes
 ^{1}To whom correspondence should be addressed. Email: jas{at}ma.hw.ac.uk

Author contributions: J.A.S., M.J.S., and J.D.M.R. designed research, performed research, and wrote the paper.

This article is a PNAS Direct Submission. J.G. is a guest editor invited by the Editorial Board.

The authors declare no conflict of interest.

Freely available online through the PNAS open access option.
References
 ↵
 Kopell N,
 Howard LN
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Scott SK
 ↵
 ↵
 Malchow H,
 Petrovskii SV,
 Venturino E
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Sherratt JA,
 Smith MJ
 ↵
 ↵
 ↵
 ↵
 Sherratt JA
 ↵
 Murray JD
 ↵
 Kay AL,
 Sherratt JA
 ↵
 ↵
 ↵
 Kopell N,
 Howard LN
 ↵
 ↵
 ↵
 Sherratt JA,
 Lewis MA,
 Fowler AC
 ↵
 Petrovskii SV,
 Malchow H
 ↵
 ↵
 Guckenheimer J,
 Holmes P
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Turchin P
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Briggs RJ
 ↵
 ↵
 Brevdo L
 ↵
 ↵
 ↵
 ↵
 May RM
 ↵
 Friedman A
 Cosner C
 ↵
Citation Manager Formats
More Articles of This Classification
Physical Sciences
Related Content
 No related articles found.
Cited by...
 Pattern selection and hysteresis in the Rietkerk model for banded vegetation in semiarid environments
 Speed of Invasion of an Expanding Population by a Horizontally Transmitted Trait
 Towards the rational design of synthetic cells with prescribed population dynamics
 Pattern solutions of the Klausmeier model for banded vegetation in semiarid environments II: patterns with the largest possible propagation speeds
 Identifying the dynamics of complex spatiotemporal systems by spatial recurrence properties