Edited by Ran Nathan, The Hebrew University of Jerusalem, Jerusalem, Israel, and accepted by the Editorial Board June 2, 2008 (received for review February 26, 2008)
The movement paths of individuals over landscapes are basically represented by sequences of points (xi, yi) occurring at times ti. Theoretically, these points can be viewed as being generated by stochastic processes that in the simplest cases are Gaussian random walks on featureless landscapes. Generalizations have been made of walks that (i) take place on landscapes with features, (ii) have correlated distributions of velocity and direction of movement in each time interval, (iii) are Lévy processes in which distance or waiting-time (time-between steps) distributions have infinite moments, or (iv) have paths bounded in space and time. We begin by demonstrating that rather mild truncations of fat-tailed step-size distributions have a dramatic effect on dispersion of organisms, where such truncations naturally arise in real walks of organisms bounded by space and, more generally, influenced by the interactions of physiological, behavioral, and ecological factors with landscape features. These generalizations permit not only increased realism and hence greater accuracy in constructing movement pathways, but also provide a biogeographically detailed epistemological framework for interpreting movement patterns in all organisms, whether tossed in the wind or willfully driven. We illustrate the utility of our framework by demonstrating how fission–fusion herding behavior arises among individuals endeavoring to satisfy both nutritional and safety demands in heterogeneous environments. We conclude with a brief discussion of potential methods that can be used to solve the inverse problem of identifying putative causal factors driving movement behavior on known landscapes, leaving details to references in the literature.
The movement of all organism, but especially sentient animals, is a complex process that depends on both an individual's ability to perform various tasks and the nature of the landscape through which it moves (1⇓⇓⇓⇓–6). These tasks include the individual's intrinsic ability to move in different ways (e.g., a horse walks, trots, canters, and gallops), the individual's internal state to perform certain activities (e.g., forage, head home, flee, and seek a mate), and the individual's ability to sense its environment, remember landmarks, construct mental maps, and process information (7). Landscape variables that influence movement include topography, abiotic variables (8, 9), location of resources (10), conspecifics by gender and age, and heterospecific competitors and predators.
Emerging digital and communications technologies have refined our ability to measure movement at the resolution of fractions of seconds with concomitant spatial precision (11), while kinematical (e.g., acceleration), physiological (e.g., heart beat and temperature), and behavioral (e.g., vocalizations) information are simultaneously recorded. The internal state driving movement, however, remains largely hidden: in animals, for example, states of hunger, thirst, and fear are either inaccessible or, at best, only indirectly inferable.
The sampling frequency of movement data affects our ability to detect short-duration fundamental movement elements (FME) (6), such as a lunge versus a step taken at normal speed. Such fine-scale events are generally unrecoverable from sampling movement at intervals coarser than the duration of these events (12) [see supporting information (SI) Fig. S1]. Furthermore, activities such as foraging or heading home involve a mix of FMEs such as being stationary, ambling, and walking; and these activities may differ only in the way the FMEs are strung together. If a string of normal steps interspersed with stationary periods, for example, is on the order of minutes for both foraging and heading to a target, then movement paths sampled every 10 min during either of these two activities can be distinguished only if they produce different characteristic “distance moved in each sampling interval” distributions. This suggests that to appropriately characterize movement components of an individual's path over time, we should endeavor to identify canonical activity mode (CAM) distributions that emerge from the mix of FMEs that characterize the activity in question: i.e., CAMS are composites of the FMEs (Fig. S1), and their characteristic step size and direction of heading distributions will depend on the length of sampling intervals (Fig. S2) and scale of analysis (13).
Ideally, if one assumes that FMEs are characterized purely by a fixed speed (because they relate to biomechanical traits of individuals), then one could mechanistically construct a movement path by specifying a sequence of FMEs with a direction of heading according to know distributions of sequence lengths and correlated heading directions for particular activities. Alternatively, at fixed points in time one could specify the next location of an individual by drawing “distance moved” and “heading direction” from empirical step-size and heading distributions (1, 6, 14⇓⇓–17). To date, such distributions are invariably derived from sampling intervals considerably longer than the shortest FME. Thus, with most current data, it is not possible to construct distributions in terms of strings of FMEs, but only in terms of longer-lasting CAMs. Consequently, CAMs are currently the preferred place to start developing a framework for movement analysis, despite the fact that any set of CAMs is unlikely to account for all of an individual's time. In many cases, however, a reasonable tradeoff may exist in defining several CAMS that account for most of an individual's time.
The most basic set of data that can be collected on the movement path of an individual is a sequence of positions (xi, yi) at points ti, i = 0, 1, 2, …, n: that is, the set 𝒟 = {u0, …, un} with ui = (xi, yi), where xi = x(ti) and yi = y(ti). For convenience, we set t0 = 0 but do not necessarily require all of the time intervals τi = [ti−1, ti] to be of equal length. In fact, in some analyses, the focus is on the distribution of waiting times τi associated with the events of first being in position (xi−1, yi−1) and then next in position (xi, yi) (14). These data can be transformed into polar coordinates (6) to obtain a set of vectors zi = (ti, di, θi) with
and
A considerable body of diffusion and stochastic process theory exists to analyze such data in the context of featureless landscapes, including uncorrelated and correlated random walks (14⇓⇓–17) and super- and subdiffusive Lévy walks (18⇓⇓–21) and Lévy modulated correlated walks (22).
Consider the distribution of waiting times τi, velocities vi, and absolute displacements di associated with a set of displacement data 𝒟. If the distributions of both waiting times and velocities have finite mean and variance, then the stochastic process associated with the data is said to be diffusive. From theory (16⇓⇓⇓–20), this implies that over an ensemble of K sets 𝒟k, k = 1, …, K, where each set is one realization of the same stochastic process, the mean-square displacement (msd) ψi of the series dik averaged over all K sets is asymptotically linear, that is ψi = (1/K)Σk=1Kdik2 ∼ t. By contrast, if upon plotting this relationship we find that ψi ∼ tp for p > 1 or < 1, then the walk is respectively referred to as superdiffusion or subdiffusion. As we demonstrate below, using simulated data drawn from a modified Pareto power law distribution (23) (Fig. S3) and plotted on a log–log scale, such plots can be misleading if time is not sufficiently large. The reason is that, initially, the relationship is affected by the fact that the distance moved in the first sampling interval is controlled by the actual step-size distribution, but from the second step onwards the distance moved from the origin is now also affected by turning angles.
Over the past decade, analyses of the movement paths of several organisms, including albatrosses (24) and spider monkeys (14), have concluded that the associated movement processes are superdiffusive, although a re-analysis of these data refute this finding (ref. 1, but see ref. 25). A possible source of error in estimating p in the msd relationship ψi ∼ tp (Fig. 1) is t must be sufficiently large for the asymptotic value to emerge. Another source of error is that superdiffusion predicts fractal looking movement paths (14, 20, 24). Using one of several different methods to estimate the fractal dimension (26) of such paths, some organisms have been pronounced as superdiffusive (14). In these organisms, however, repeated fractal patterns occur at no more than two or three particular scales and are more a feature of the way resources are distributed across the landscape than of a genuine superdiffusive process. Thus, it is important to understand how animal movement is influenced by landscape features and to assess the extent to which step-size distributions are modified by landscape heterogeneity.
The logarithm of mean-square displacements (msd) ψ versus the logarithm of time t average over 10,000 simulations are plotted for a random walk with step lengths drawn from modified Pareto distributions (Upper Left, q = 2; Upper Right, q = 1.5; Lower Left, q = 1) and directions for each step completely random. From the lines inserted by “eye” (red, small t; blue, large t), Upper Left represents diffusion (P = 1 for large t), and Upper Right and Lower Left represent super diffusion (respectively, P = 1.15 > 1 and P = 2 for large t). In Lower Right, for the case q = 1, the length of these excursions are truncated at a step size of 100 (biologically, an upper bound is set by the maximum velocity of the organism multiplied by the length of the time interval), which is far out in the tails of the distribution (two orders of magnitude beyond the mode; see Fig. S3). In this case, the noisy superdiffusive behavior is completely tamed even though, initially, it looks superdiffusive (ψ ∼ t1.4 for t ≤ 2).
We note here that because all bounded step-size distributions produce Gaussian random walks when turning angles are uncorrelated, critical information in the step-size distribution, such as multimodality arising from mixed distributions of FMEs, is lost when using the statistics of emergent global characteristics such as msd as a function of time. This stresses the importance of knowing the actual step-size distributions when deconstructing movement paths and trying to understand the causal processes creating local path structures.
The framework presented here is formulated in the context of a group of N known individuals indexed by k. This specificity allows us to account for the following: (i) species, gender, and age-specific differences; (ii) unique memory and knowledge of landscape; and (iii) cues and vectors individuals use to select a new position on the landscape. Further, we assume that each individual has nmk FMEs (Table 1), each with its own characteristic speed sjk, j = 1, …, nmk, from which its movement track is generated. The set of FMEs constitute the basic motion capacity Ωk (7). Furthermore, in any segment of the movement path (Fig. S1), these FMEs are mixed in various proportions to constitute a set of Ark CAMs (Table 1), r = 1, …, nak best characterized by distributions of speeds (equivalently distances) with means and standard deviations s̄r and στr for each r (Fig. S2). We note, however, that when the sample intervals τik = τk are fixed over all intervals i, for reasons discussed below only the standard deviation (and not the mean) depends on sampling frequency 1/τk.
Definition of parameters and variables
As an example, Ωk in horses has been defined in terms of five FMEs or gaits (27) of increasing speeds: stationarity (s1 = 0), walking (a four beat gait with s2 ≈ 4 mph), trotting (a two beat gait with s3 ≈ 8 mph on average), cantoring (a three beat gait for which s4 > s3), and galloping (a four beat gait for which s5 ≈ 25–30 mph, varying across horses). Beyond these natural gaits, some horses have been bred to implement so-called “ambling gaits” that provide a smoother ride, but not all horses can execute these, thereby providing an example of selection at work on the motion capacity component Ω of our guiding conceptual model (7).
For each of the nak, CAMs associated with individual k, the distribution of speeds (equivalently distances or “step sizes”) associated with a particular CAM is affected by the value of τ: if τ is small an individual will generally be in one movement mode or another in the mix that constitutes the activity, whereas if τ is large an individual will more likely be executing some mix of modes over one interval (Fig. S2). Thus, even though the proportion of different movement mode events used to construct the movement tracks may be quite stable with regard to a specific activity, the directions of heading from one event to the next make the relationship between the length of a track and displacement quite complicated. Only for very simple cases, such as tracks generated from one type of FME, is the relationship between length-of-track and displacement easily cracked by analytical methods. For the rest, the easiest way to generate the distribution of speeds (distances) Aτr(s) as a function of sample interval size τ is to use Monte Carlo simulation.
The proportions of FMEs in a particular CAM is likely to vary with changes in landscape. For example, the mix of stationary, walking, and running modes used by a foraging African antelope will differ in an open savannah compared with thick bushveld. This problem can be dealt with by assuming that the parameters of a CAM distribution depend on external landscape factors in addition to the frequency of sampling along a movement path.
One of the goals in an analysis of empirical data is to see how cleanly a set of CAM distribution of step-size FMEs can be extracted from the movement track to explain the actual activity producing different segments of the movement track. Some segments may reflect a pure activity (e.g., foraging) while others are a mix of activities (e.g., foraging and resting) or even a compromise between competing activities such as an individual foraging as it heads to a known water source or home. To be able to account for such mixed activities, as well as assess factors that may lead to an individuals switching from one CAM to another, we must infer that the individual has a multidimensional internal state that drives the behavior (7). The current state of this driver can in turn be associated with a goal emerging from an individual's internal state, which in general will vary with time. The relationship between an individual's internal state (i.e., the vector w—see Fig. S4) and its current goal state can be treated in various ways. One way is to construct a mapping of a continuous n-dimensional vector space to a discrete space of g goals. Another is to consider the internal state as a weighting vector
that produces a goal-modified “ideal” distribution Gτik(s) of speeds (distances) s for individual k at time ti = iτ with the weighted sum of its nak CAMs: i.e.,
In the case of organisms that have no internal mechanism for generating goals (e.g., plants), the internal state may, for example, represent elements controlling dehiscence of seed dispersal structures.
The spatial information embedded in the distribution Gτik(s) needs to be made explicit before other relevant spatially explicit geographical and biological information can be incorporated into the movement process. The distribution Gτik(s) is defined on s ∈ [0, smax]. On a flat structureless landscape—that is, in the absence of all cues and other directional biasing factors beyond the context of a correlated walk—an individual is equally likely to move in any direction. Thus, the ideal (i.e., featureless landscape) distribution Gτik(s) can be given an explicit spatial dimension by rotating it around the point [0, 0] on a plane parameterized by coordinates (α, β) to obtain a radially symmetric distribution 𝒢τik(α, β) that has a top-half-of-a-donut-like structure (Fig. S5). The distribution can now be relocated so that its center of symmetry is the current point of location uik = (xik, yik) of organism k.
We incorporate landscape features into each individual's goal as follows. First, we cover the landscape with a raster of rectangular cells Cαβ. Second, we locate the position uik = (xik, yik) of each individual within the cell containing this point: in general, this allows several individuals of one or more types to be located in each cell. Also, we associate an nc-dimensional external state to account for all factors relevant to the movement of each of the N individuals on the landscape, including other individuals. Third, the one-dimensional CAM distributions Aτrk(s) are used, as described above, to generate two-dimensional radially symmetric, but discretized, distributions 𝒜τrk(α, β) ≥ 0 with Σ(α,β)𝒜τrk(α, β) = 1 for all k and r = 1, …, nak. Fourth, we incorporate the landscape effects through a set of nak landscape modifier matrices (LMM) Lirk (i.e., specific to individual k, their activity r, and time iτ) with elements ℓir,αβk constructed from those elements of the external state vector that are applicable for the activity in question. These elements ℓir,αβk are used to modify the CAM distributions to reflect the preference that each individual has for each of the landscape cells while involved in one of its CAMs (e.g., one cell may be the most desirable from a foraging point of view, whereas another is desirable as a target when heading for water). These LMMs are used to modify the movement distribution 𝒢τik(α, β) to account for the state of the landscape, as well as individual specific information that inter alia relates to the location of places [remembered, sensed, and inferred through landscape cues, including the earth's magnetic field (28)] associated with heading activities, the distribution of resources across the landscape, and the location of other individuals on the landscape; and hence the elements of these matrices must be updated each time step.
The simplest way to incorporate this landscape information into the movement process is to use the activity-specific LMM elements ℓir,αβk (Table 1) as a way to weight the terms in the idealized goal determined movement distribution Gτik(s) to produce a composite realized movement distribution, with spatial matrix elements
for all α, β, i, and k = 1, where
normalizes the discrete elements of the probability distribution over the cells Cαβ for each individual k and all time i.
The final component of the movement process is how the organism selects the particular cell (α, β) that will determine its next position ui+1k = (xi+1k, yi+1k) at time t = (i + 1)τ. The simplest rule is for an organism to move to cell (α, β) with largest value Mτik(α, β), effectively without making any interim decision on its way to the new cell. For the case of an organism driven purely by stochastic landscape processes (such as wind), movement to the next cell can be regarded as random with the probability of selecting a particular cell (α, β) equal to the values Mτik(α, β). A third possibility is to select cells with a probability that is proportional to some power of Mτik(α, β), a solution that is intermediate between the first two if this power is >1. If this power is <1, the solution is intermediate between the second mechanism and a purely random solution that gives no weight to the relative values Mτik(α, β).
The following example illustrates the application of our framework to simulating the movement of a herd of social ungulates foraging on a heterogeneous landscape with conflicting needs to both assuage hunger and remain safe by staying close to other individuals.
In social ungulates, the existence of groups is presumably the result of safety offered by the group (29). However, membership in a group comes at a cost, namely competition for limited resources (usually food). These conflicting needs presumably drive the observed fission–fusion dynamics typical of many social ungulates (30). Thus, the goal for each individual is feeding while remaining near conspecifics, and the main internal states for this goal are levels of hunger and safety, with each individual seeking to assuage hunger while remaining safe. The CAM in this example is pure foraging, and the FMEs allow either moving between patches constrained by a maximum distance traversed in each time set, or remaining in the current patch while feeding (see SI Text). Navigation in our example is relatively simple: individuals move directly using visual information to a cell selected within a fixed radius interpreted as the observable range. The distribution of food patches and the location of other herd members represent the external states that a given herd-member responds to in its choice of direction. The external states are combined by using LMMs to obtain the realized movement distribution from which an individual selects its next target.
Our model is spatially explicit with each grid cell, in this particular case roughly the size of an individual. For each time step an animal may perform one of three activities: (i) moving toward a selected patch (one cell per time step), (ii) feeding in a patch (consuming one unit per time step), or (iii) resting in a patch when it has a full gut. There are two matrices, each representing one of the external states: a food (vegetation) matrix and a safety matrix (Fig. 2Upper). The vegetation matrix is static and is generated at the start of each simulation by using a combination of random procedures to create a patchy heterogeneous landscape (see Materials and Methods and SI Text). The safety matrix is a function of the spatial location of all group members (Fig. 2), risk does not vary with landscape features, and there are no visible predators. In each time step, a safety score is calculated to each cell as the sum of the inverse of the distance of all herd members to that cell (see Materials and Methods). The LMM is a weighted sum of the vegetation and safety matrices, where the weightings α and (1 − α) depend on the internal driver (state) α ∈ [0, 1] reflecting the current priorities to the individual in trading safety against hunger. When an animal stays within a patch to feed, its value of α increases and when it moves its value of α decreases so that the relative importance of food to safety oscillates up and down as the animal approaches satiation or becomes increasingly hungry (see Materials and Methods and SI Text).
Each of the five panels is an extract from a much larger mapping of the values of elements in the landscape modifier matrices (LMMs) for the vegetation and safety landscapes of the realized discretized movement distributions constructed from these matrices. The focal individual is represented by the small red squares in each of the five panels, with the positions of its conspecifics represented by other small squares in Upper Left and Lower. In the distributions represented in Lower, the most attractive areas are the lighter areas (but ignoring the small dark blue squares, which are just conspecific position markers for reference). The relative weighting of safety over resources ranges from safety being the only consideration (Lower Left), safety and resources being equally important (Lower Center), and resources being the only consideration (Lower Right). Imposed upon Lower is a circle representing the maximum possible movement displacement in one time step.
By playing around with simulation parameters, it is possible to test how various movement-related questions—such as fission–fusion dynamics, subgroup structure, and trajectory patterns—are affected by patch size, gut size, and feeding strategies. For example, we explored how patch density affects fission–fusion and subgroup structure as follows: 20 animals with equal gut sizes were initially spread randomly over a 20 × 20 cell section of a 220 × 1,100 cell grid, with movement rules detailed in SI Text. We assumed that the vegetation in each cell is not renewable, so that depleted patches act as repulsive regions. The narrow landscape channels the herd into a directional movement along the long axis, but the short axis is sufficiently wide to enable fission events (Fig. 3).
Simulation of fission–fusion behavior as a function of vegetation quality. Open squares, high quality; filled diamonds, low quality. Population defined to be in a two-herd (one herd) state when the fission index is <0.5 (>0.5). See main text and SI Text for details.
The simulation was repeated with the same starting locations with two different vegetation matrices, each with a mean patch size of 38 cells. The first, representing “good” habitat, produced a landscape where ≈75% of the cells were positive (Fig. 3), whereas the second, representing “poor” habitat, produced a landscape with ≈65% positive cells. We ran each simulation for 10,000 time steps. Group structure was analyzed every 500 steps, starting with the 1,000th step, using cluster analysis with a two-group restriction (see SI Text). We then calculated the mean individual distance for the entire herd and for each subgroup. We considered a ratio of 0.5, between the sum of the two subgroup mean distances to the mean distance of the entire herd, as an indicator of fission. Simulating the herd's movement by using the “poor” habitat matrix produced three fission events compared with none for the “good” habitat (Fig. 3). All fissions were followed by fusion after 1,000–1,500 time steps. We note that the herd was generally cohesive, with fission and fusion events emerging rather than explicitly constructed, despite that fact that the movements of the individuals themselves were completely deterministic. Our approach highlights the importance of quantifying the internal drivers in understanding animal movements. In the context of our specific simulation, empirical data on giving-up densities (31) appear to be a good tool for quantifying internal drives and understanding movement patterns.
The framework presented here provides a mathematically detailed exposition of the conceptual model developed by Nathan et al. (7) that can be used for both the construction and deconstruction of the pathway of organisms on structured landscape, the former through simulations and the latter through state-space methods (32) for fitting parameters. We have not dealt in any detail with the problem of movement mode identification that is the key to the deconstruction component other than stressing that the frequency with which data are collected limits our ability to identify CAMs and their underlying FMEs (12). To undertake such an analysis is not a trivial problem: it requires computationally complex methods that can only be successfully applied to high-resolution data, but a start has been made (33⇓–35).
Before automated GPS data collection, VHF telemetry position data were typically collected too infrequently and were also not sufficiently accurate to be useful for reconstructing pathways; but these data were suitable for constructing home ranges or types of utilization distributions at a seasonal scale, using both parametric (36⇓–38) and nonparametric (39, 40) kernel methods as the preferred methodology. Since the mid-1990s (24) movement data have been fitted to Lévy models to evaluate the extent to which this movement is superdiffusive with, as we have discussed, mixed success (1) and also to assess the degree to which path characteristics can be fitted by correlated random walks (17) or mixtures of random walks (4), but these analyses do not explicitly incorporate landscape structure or factors.
In terms of general methods for deconstructing movement pathways, various time series and frequency domain techniques can be brought to bear on the problem under the rubric of exploratory data analysis (EDA) (41). Furthermore, stochastic differential equation methods (42) can be used to construct vector fields from data on the contemporaneous movement of many individuals, and then thin plate splines can be used to fit potential fields to these vector fields to identify regions on the landscape that are either repealing or attracting the individuals at a particular time of day (2). Recently, techniques new to the field of movement ecology, such as wavelet analysis (33) and artificial neural networks (34), are being applied to obtain insights into the effects that the internal and environmental states of a system have on movement paths. Beyond these, as the resolution of movement and landscape data improves dramatically over the next decade, we should expect to see the application of state-space estimation methods (32) that can take advantage of formulations such as ours, because our formulation permits the inclusion of detailed landscape information.
Details are elaborated on in SI Text.
This matrix is the only component in the model that has stochastic elements. It consists of a series of patches set up using Monte Carlo methods, with a parameter controlling patch density and a beta distribution controlling patch size. The quality of resources in the center of each patch was then assigned a number at random between 1 and 10 with values declining to the edge of the patch. The resource value of each cell was reduced by a set amount in each time step for which the patch was occupied by a feeding individual. Fig. 2 Upper Left depicts the result of one such construction.
For each of the cells containing an individual (i.e., focal individual), a matrix of values for all of the remaining cells was constructed. The values associated with each of these remaining cells is based on the sum of the inverse distances of all of the remaining organisms to these cells. Thus, cells of highest values are those closest in an integrative sense to the organisms as a group that excludes the focal individual (see Fig. 2 Upper Right).
All points within a fixed distance of an individual (and only these points) were regarded as selectable targets to move to next (circles in Fig. 2). An individual then moves toward the cell that has the highest value of a weighted sum of the vegetation and safety matrix values for that cell.
Each animal has an energy bank that determines its level of hunger and, hence, its weighting of vegetation and safety matrix values. During each time step, an animal can either feed and increase its energy bank or move and decrease its energy bank. Once an individual has consumed the vegetation locally, it moves to a patch within navigation range that maximizes its current tradeoff for resources versus safety.
We thank the members of the Institute of Advanced Studies movement ecology group for ideas that we have absorbed into our text during the course of our weekly meetings over many months, and Leo Polansky for comments on the manuscript. This work was supported by the Institute of Advanced Studies at the Hebrew University of Jerusalem and by a James S. McDonnell 21st Century Science Initiative Award.
Author contributions: W.M.G. and D.S. designed research, performed research, contributed new reagents/analytic tools, analyzed data, and wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission. R.N. is a guest editor invited by the Editorial Board.
This article contains supporting information online at www.pnas.org/cgi/content/full/0801732105/DCSupplemental.
