## Abstract

Some of the strongest empirical support for Lévy search theory has come from telemetry data for the dive patterns of marine predators (sharks, bony fishes, sea turtles and penguins). The dive patterns of the unusually large jellyfish *Rhizostoma octopus* do, however, sit outside of current Lévy search theory which predicts that a single search strategy is optimal. When searching the water column, the movement patterns of these jellyfish change over time. Movement bouts can be approximated by a variety of Lévy and Brownian (exponential) walks. The adaptive value of this variation is not known. On some occasions movement pattern data are consistent with the jellyfish prospecting away from a preferred depth, not finding an improvement in conditions elsewhere and so returning to their original depth. This ‘bounce’ behaviour also sits outside of current Lévy walk search theory. Here, it is shown that the jellyfish movement patterns are consistent with their using optimized ‘fast simulated annealing’—a novel kind of Lévy walk search pattern—to locate the maximum prey concentration in the water column and/or to locate the strongest of many olfactory trails emanating from more distant prey. Fast simulated annealing is a powerful stochastic search algorithm for locating a global maximum that is hidden among many poorer local maxima in a large search space. This new finding shows that the notion of active optimized Lévy walk searching is not limited to the search for randomly and sparsely distributed resources, as previously thought, but can be extended to embrace other scenarios, including that of the jellyfish *R. octopus*. In the presence of convective currents, it could become energetically favourable to search the water column by riding the convective currents. Here, it is shown that these passive movements can be represented accurately by Lévy walks of the type occasionally seen in *R. octopus*. This result vividly illustrates that Lévy walks are not necessarily the result of selection pressures for advantageous searching behaviour but can instead arise freely and naturally from simple processes. It also shows that the family of Lévy walkers is vastly larger than previously thought and includes spores, pollens, seeds and minute wingless arthropods that on warm days disperse passively within the atmospheric boundary layer.

## 1. Introduction

A diverse range of organisms, including *Escherichia coli*, T cells, mussels, honeybees, some marine predators and even human hunter–gatherers have movement patterns that can be approximated by Lévy walks when searching for prey, conspecifics or home [1–9]. Lévy walks (also called Lévy flights in the biological literature) are random walks in which the step-lengths, *l*, have a probability distribution with a heavy power-law tail, *p*(*l*) ∼ *l ^{−μ}* with 1 <

*μ*≤ 3. The heavy power-law tail causes the measured mean-squared step-length to increase over time [10] and this implies the absence of a characteristic scale. This absence of a characteristic scale results in a fractal pattern of ‘walk clusters’ such that clusters of small-steps are interspersed with long steps to new locations, with this pattern repeating across all scales. It has long been realized that this fractal scaling can be advantageous when searching for sparsely and randomly distributed resources because it reduces oversampling [11]. This realization has given rise to the ‘Lévy flight foraging hypothesis’ [12]. This states that since Lévy walks can optimize search efficiencies, natural selection should have led to adaptations for Lévy walk foraging. The most convincing empirical support for this hypothesis comes from telemetry data for the dive patterns of a wide variety of marine predators, including sharks, bony fishes, sea turtles and penguins [4,5]. In accordance with theoretical expectations, these movement patterns are well represented by Lévy walks with

*μ*≈ 2.

*Rhizostoma octopus*, an unusually large jellyfish weighing 27 kg or more that is capable of actively swimming against localized currents, also perform extensive vertical movements that equate to moving up and down through the entire water column dozens of times a day [8]. These movements have been attributed to active feeding rather than to other activities such as reproduction and escape from predators. These dive patterns do, however, sit outside of current Lévy search theory [8]. This is because the jellyfish showed switching behaviour from exponential patterns to patterns better fitted by a variety of truncated Lévy walks (mean

*μ*= 1.96, range 1.2−2.9). It is also because in some periods, vertical excursions were followed by a vertical return to the depth occupied prior to the excursion. The jellyfish movement patterns appear to be consistent with their using optimized ‘fast simulated annealing’—a novel kind of Lévy search strategy—to locate the highest prey concentration in the water column and/or to locate the strongest olfactory cues emanating from more distant prey.

Simulated annealing is a powerful stochastic search algorithm for locating a global minimum (or maximum) that is hidden among many poorer local minima (maxima) in a large search space. It works by emulating physical annealing whereby a solid is slowly cooled so that when it eventually ‘freezes’ it does so at a minimum energy configuration. Kirkpatrick *et al*. [13] and Černý [14] were the first to establish a connection between this type of thermodynamic behaviour and the search for a global minimum (or maximum) of a cost function that may possess several local minima (or maxima). Simulated annealing has since been applied to numerous optimization problems [15]. The first applications used the so-called ‘Boltzmann’ annealing which is based on strictly local sampling, as step-lengths in the search patterns are Gaussian distributed. A major drawback with Boltzmann annealing is its slow rate of convergence. The problem was overcome by Szu & Hartley [16], who introduced ‘fast simulated annealing’ in which searching patterns are like Lévy walks with *μ* = 2. Fast simulated annealing will, in fact, outperform all Boltzmann-like simulated annealing schemes, i.e. all schemes that use thin-tailed step-length distributions.

Fast simulated annealing would be an effective searching strategy when *R. octopus* are presented with discrete prey layers rather than with a relatively uniform prey field with no effective prey maxima and minima, conditions that may suppress active searching behaviours. This speculation is consistent with ‘Lasker's stable ocean hypothesis’ [17]. This states that when plankton prey are localized in high density, discrete depth layers, then fish larvae focus on these layers and achieve high foraging success compared with when plankton prey are more dispersed, for example, after storm events.

Fast simulated annealing is, however, unlikely to be advantageous in the presence of convective currents, because its execution would require that the jellies actively resist being swept along by these currents and because the prey distribution will become dynamic. Turbulent convection in warm littoral waters occurs in response to a cold air outbreak, e.g. to the passage of a cold atmospheric front, and is well documented [18]. Taking advantage of convection currents would be an energetically favourably way to search the water column, even when some active swimming in the form of slow pulsing is required (*R. octopus* are nearly neutrally buoyant).

Here with the aid of numerical simulations I examine the extent to which the movement patterns of *R. octopus* are consistent their using fast simulated annealing and with their occasionally riding convective currents.

## 2. Material and methods

### 2.1. Fast simulated annealing

The fast simulated annealing algorithm of Szu & Hartley [16] is a simple adaptation of the Metropolis–Hastings algorithm [19], a Monte Carlo method used to generate sample states of a thermodynamic system. A cost function, e.g. the prey density, *c*(*x*), at the current position, *x*, is compared with the cost function, *c*(*y*), at another position, *y*, where the step-length, *l* = *|y* − *x|*, is drawn at random from a Cauchy distribution
2.1and where the ‘temperature’, *T*, is a measure of the size of the fluctuations in step-length. The new position is accepted with an ‘acceptance probability’ which often (and herein) takes the form
2.2where Δ*c* = *c*(*y*) − *c*(*x*) is the cost difference between the present and the previous positions. The new position is necessarily accepted if it is better than the previous one, i.e. if Δ*c* > 0 and is occasionally accepted if it is worse than the previous position. If the new position is not accepted, then the searcher returns to its previous position. The cycle is then repeated many times and at each step the temperature is gradually lowered. This annealing (cooling) corresponds to a slow decrease in the probability of accepting worse solutions as the searcher explores the solution space. Accepting worse solutions is a fundamental property of simulated annealing because it allows for a more extensive search for the optimal solution. Szu & Hartley [16] showed that this fast simulated annealing converges to the global maximum if annealing is scheduled according to *T*(*k*) = *T*_{0}/*k*, where *T*_{0} is the initial temperature and *k* is the step count. The efficiency of the algorithm can be increased by occasional ‘reheating’ which resets the iteration count. This ‘reheating’ allows the searcher to quickly escape from local maximum.

Szu & Hartley [16] focused their attention exclusively on reducing the expected number of steps in a search pattern because this is of crucial importance in computational settings. But in a foraging context, the overall length of a search path is crucial because it determines energy expenditure. There can only be selection for fast simulated annealing if the overall lengths of the resulting search paths are shorter than those produced by other processes. Here, this possibility was examined with the aid of numerical simulations. The effectiveness of the original fast simulated annealing algorithm of Szu & Hartley [16] was compared with other fast simulated annealing algorithms corresponding to different kinds of Lévy walk search pattern. Step-lengths in these search patterns were drawn from distributions
2.3with 1 < *μ* ≤ 3. The case with *μ* = 3 effectively corresponds to Boltzmann simulated annealing because the step-lengths have a characterized size. As an illustrative case, simulated prey concentration profiles in water columns were generated using a simple, Gaussian random walk model, i.e. incremental changes in prey concentration with changes in depth were Gaussian distributed but constrained to prevent concentrations from becoming negative (figure 1*a*). It is assumed that the prey distribution does not change significantly during the time course of the search. Later it will be shown that *R. octopus* could potentially execute several fast simulated annealing searches in a single day and so the assumption of a stationary prey distribution need only hold for a few hours*.* Robust predictions for the average length of search paths were obtained by ensemble averaging over 10^{5} searches of a given type.

### 2.2. Modelling turbulent convection

The objective of this modelling is to assess whether physical transport in the water column (independent of swimming) can yield movement patterns like those occasionally seen in *R. octopus. Rhizostoma octopus* could avail of convective currents but are not constrained to do so. Turbulence within convective flows is inhomogeneous and non-Gaussian [20]. The trajectories of passively advected, neutrally buoyant particles (a proxy for passively advected *R. octopus* and hereafter called ‘tracer-particles’) within such flows can be predicted accurately by Lagrangian stochastic (particle-tracking) models [21]. Thomson [21] showed how the evolution of the position, *x _{i}*(

*t*) , and velocity,

*u*(

_{i}*t*), of a tracer-particle in a turbulent flow with high Reynolds number can be represented in its most general form as an Ito equation, i.e. as a modified Langevin equation 2.4where

*C*

_{0}is the universal constant appearing in the Kolmogorov inertial subrange form for the Lagrangian velocity structure function,

*ɛ*(

*z*) is the mean rate of dissipation of turbulence kinetic energy, the indices denote Cartesian coordinates and d

*ξ*

_{i}(

*t*) is an incremental Wiener process with correlation property . The deterministic term,

*a*(

_{i}*u*,

*x*,

*t*), is determined by the Eulerian statistics of the flow via Thomson's [21] ‘well-mixed condition’. Thomson's [21] well-mixed condition is equivalent to or more stringent than all other criteria which have been used to distinguish between well and poorly formulated Lagrangian stochastic models.

Here, tracer-particle trajectories within a convective boundary layer were simulated using a one-dimensional Lagrangian stochastic model developed by Luhar & Britter [22]. Luhar & Britter [22] provide a full description of this model and its numerical implementation which is not repeated here. The model of Luhar & Britter [22] simulates vertical tracer-particle movements in stationary and horizontally homogeneous turbulence in a convective boundary layer with no mean flow. More elaborate simulations were also performed using a two-dimensional Lagrangian stochastic model for dispersion (stream-wise and vertical tracer-particle movements) in convective boundary layers with shear (horizontal flow) [23].

Simulated tracer-particle trajectories were partitioned into discrete bouts of upward and downward movements that take place within up- and downdrafts. The Akaike information criterion [24] was then used to test whether these simulation data are better represented as power-laws, as bi-exponentials, or as single exponentials, i.e. whether the simulated movement patterns can be regarded as being approximations to Lévy walks, bi-modal walks or as approximations to Brownian walks. Model distributions were fitted to the simulation data using log-maximum-likelihood methods [25]. The goodness of fit of the model distributions to the simulation data was determined by visual inspection of the survival functions (rank frequency distributions). This was sufficient because the model distributions differed significantly and because the best-fit distribution was clearly evident. Unlike probability distributions, survival functions can be constructed without first having to bin data; a procedure that can affect the shape of a distribution [26]. To construct the survival function, the simulation data for the step-lengths {*l _{i}*} were first ranked from largest to smallest {

*i*= 1, … ,

*n*}. The probability that a step-length is greater than or equal to

*l*(the survival function) and the probability that a step is greater than or equal to

_{i}*l*were then estimated as

_{i}*i*/

*n*.

## 3. Results

### 3.1. Fast simulated annealing

The initial temperature (initial sale of fluctuations in step-length) was taken to be *T*_{0} ≈ 0.1 *H*, where *H* is the depth of the water column and reheating (i.e. resets) was performed after every cycle of 50 iterations. These parameter values were chosen because they were found numerically to produce the most efficient searches of the simulated prey fields. These values are dependent upon the statistical structure. Different optimized parameter values were obtained when the simulated prey profile comprised discrete layers separated by voids (data not shown). Fast simulated annealing with Lévy walks having *μ* ≈ 2.5 are optimal for the location of the maximum prey concentration in the simulated Gaussian profiles (figure 1*b*). With this searching strategy, the theoretical maximum prey concentration is located after a searcher has traversed the water column around a dozen times. Fast simulated annealing with Lévy walks having *μ* ≈ 2.5 were also found to be optimal for some other kinds of random prey distributions, e.g. when incremental changes in prey density are Lévy distributed with *μ* = 2 rather than Gaussian distributed. Such universality with respect to the scaling properties of the targets is not unprecedented and has previously been seen in numerical simulations of Lévy walk searches in fractal environments [27]. The presence of Lévy walks with *μ* ≈ 2.5 in an optimal fast simulated annealing search may, however, not be evident at all times because the step-length distribution evolves over the time course of a search due to ‘cooling’ (figure 2). This switching behaviour is an artefact of keeping the observational range of spatial scales fixed and is compounded by the fact that the distribution of step-lengths for any time window is effectively an aggregation of many different distributions. Vertical excursions followed by returns to the depth occupied prior to the excursion are frequent and especially so during the later stages of a search when step-lengths are relatively short (figure 2).

### 3.2. Turbulent convection

Model predictions for the mean height and the mean spread of tracer-particles in convective flows are in close accordance with experimental data from classic controlled water tank experiments of turbulent dispersal in convective boundary layers [20,28,29] (electronic supplementary material, figure S1). This shows that the model of Luhar & Britter [22] can simulate accurately tracer-particle trajectories in convective boundary layers and has been *encoded* correctly. Luhar & Britter [22] obtained the same degree of model agreement with the water tank data. Simulation data for the lengths of upward and downward movements are best represented by truncated power-laws and as a consequence the simulated movement patterns can be regarded as being approximations to *asymmetric* Lévy walks (figure 3). The maximum-likelihood estimates for the Lévy exponents characterizing the upward and downward movements are *μ* = 1.29 and 1.21. The Akaike weights (weight of evidence) for truncated power-laws being the better of the trial model distributions are 1.00 and 1.00. The asymmetry arises because in a convective boundary layer updrafts have large vertical velocities but occupy less horizontal area while downdrafts occupy a large area but have smaller velocities [22]. Nonetheless, Hays *et al*. [8] did not discriminate between upward and downward movements. A direct comparison with these empirical observations (and more generally with dive patterns for other marine predators [4,5]) can be made by pooling simulation data for upward and downward movements. A truncated power-law provides a near perfect fit to this pooled simulation data *r* (electronic supplementary material, figure S2). The maximum-likelihood estimate for the power-law exponent *μ* = 1.25. A truncated power-law with *μ* = 1.25 was also found to represent accurately simulation data obtained using the model of Rotach *et al*. [23] for convective conditions with shear (horizontal flow) and persists when intermittency (fluctuations in the rate of dissipation of turbulent kinetic energy, *ɛ*) are accounted following the approaches of Pope & Chen [30] and Borgas & Sawford [31]. For neutrally stable conditions, movement patterns are effectively Brownian but in accordance with the analysis of Reynolds [32] exhibit *μ* = 4/3 Lévy scaling over a limited range.

## 4. Discussion

The complex movement patterns made by *R. octopus* were here shown to be consistent with their using fast simulated annealing to locate the highest prey concentration in the water column; a search strategy that was found to be robustly optimal. Fast simulated annealing explains: (i) why the movement patterns appear to change over time; (ii) why movement bouts are characterized by a variety of Lévy walks and normal exponential walks; (iii) why *R. octopus* bounce (occasionally return to their original depth after prospecting away), a behaviour that is also seen in some fish [33,34]; and (iv) why movement bouts resembling Lévy walks with *μ* ≈ 2 are rare. Lévy walks with *μ* ≈ 2 are expected when searching indiscriminately for patchily distributed resources [11]. This correspondence between predictions and observations suggests that the jellyfish search patterns do not, after all, sit outside of the Lévy search paradigm. It supports the view that strata of high densities of *Rhizostoma* prey are rare so that jellyfish spend long periods searching for high prey densities [8]. By focusing on these layers, the jellyfish can achieve high foraging success compared with when they forage less selectively on strata with relatively low prey densities. This is consistent with the general expectation from optimal searching theory that all foragers should search within the ‘patch’ with the highest density of food rather than within patches with lower densities of food, even if all patches contain enough food to satiate the forager [35]. This is simply because searching in the highest density patch increases the probability of finding food with each search attempt, and decreases the expected time to achieve satiation. The forager thereby maximizes its' energy intake per unit time. Fast simulated annealing could also be used by *R. octopus* to locate the strongest of many horizontally dispersing olfactory trails emanating from more distant prey, as some jellyfish can orient to olfactory cues from prey [36].

An antecedent of the hallmark of fast simulated annealing in searching behaviours can be found in the movement patterns of displaced honeybees when searching for their hive, and in the movement patterns of the Australian desert ant *Melophorus bagoti* when searching for their nest in visually unfamiliar surroundings [2,37]. Both of the insects adopt a stereotypical search strategy that begins at the location where the insect initially expects to find the target of its search, and is comprised of loops that start and end at this location, and are directed in different azimuthal directions. Loop-lengths are, to good approximation, distributed according to an inverse-square power-law, corresponding to ‘looping’ (bouncing) Lévy walks with *μ* ≈ 2. This movement pattern, which is known to be optimal [38], would be produced if organisms used fast simulated annealing to locate a single target. Honeybees use the same strategy, sometimes in combination with olfactory searching, when attempting to locate a food source after a known food source becomes depleted [39]. Data from preliminary simulations suggest that two-dimensional fast simulated annealing for the location of the best of many possible food sources is near optimal when the Lévy exponent *μ* ≈ 2. The employment of fast simulated annealing by honeybees is consistent with the suggestion of Dornhaus & Chittka [40] that the ability of these bees to recruit foragers through the waggle dance evolved in order to exploit highly rewarding but ephemeral flower patches in the tropics. This would make searching for high-quality patches a strong selective force. Food quality assessment must be reasonably precise because if bees vigorously advertised relatively poor patches then the colony will decline or even die out. Similarly, long-term survival is crucially dependent upon the ability of scout bees to eventually locate the best location for a future hive. In controlled field studies, scouts were observed to continue searching after finding a candidate location and to eventually find the best location [41]. Looping Lévy walks consistent with fast simulated annealing have also been identified in bumblebees [42]. The employment of fast simulated annealing by bumblebees is consistent with numbers of nectar foraging bumblebees following ‘ideal free distributions’ [43]. Ideal free distribution theory predicts that foragers should distribute themselves on a variable resource so that the individual food intake rate is the same at all local areas or patches [44]. Dreisig [43] reported that ‘the number of nectar foraging bumblebees on plants varying in size in 12m^{2} plots within a meadow was proportional to the number of open flowers per plant and on plants varying in rate of nectar production per flower proportional to the nectar production per plant’. These distributions result in a nectar gain per flower which is independent of plant size and rate of nectar production per flower. Such ideal free distributions would be achieved if each newly arriving bee forages at the location that currently has the highest gain and in doing so reduces the expected gain of that location. Nonetheless, the adaptive value of the searching behaviour is the avoidance of poor locations rather than the equalization of gain. It is perhaps natural to assume that the dynamics underpinning ideal free distributions simply involves foragers moving or tending to move in the direction of increasing patch quality. This local dispersal rule does, however, result in trapping in local maxima and this prevents the establishment of ideal free distributions in all but the simplest of situations [45]. Other candidate mechanisms for the establishment of ideal free distributions rely on optimal patch departure rules which require that foragers can ‘learn’ about the environmental average [46]. Looping Lévy walks and fast simulated annealing overcome these limitations. The possibility that bumblebees use fast simulated annealing to locate the best patch is also consistent with the flight patterns of naive bumblebees in the presence of multiple feeders [47]. Upon finding a feeder, bumblebees made looping flights away from that feeder, as though they were searching for other feeders. Eventually, all the feeders were located in this way.

This association of *R. octopus* movement patterns with fast simulated annealing leaves open the question of the underlying mechanism, i.e. the identification of the plastic physiological or neurological process that can tune up for the execution of advantageous Lévy walk movement patterns. Candidate mechanisms for Lévy walks in marine predators such as bony fish and turtles, honeybees and even in *E. coli* have been identified [48–51]. These studies have shown that programming for the Lévy strategy does not need to be very sophisticated or clever on the organism's part, as active Lévy walk search patterns can be derived from simple reflexive behaviours, errors in distance estimation and intrinsic noise in the chemotactic pathway. *Rhizostoma octopus* could also approximate optimized Lévy walk movements as composite random walks as in the cases of mussels and the Australian desert *M. bagoti* [37,51]. In this regard, it should be noted that the most effective composite random walk searches are those that most closely resemble optimal Lévy walk searches [51].

Fast simulated annealing in *R. octopus* may not be an advantageous search strategy when convective currents are present because its execution would require that the jellies actively resist being swept along by the convective currents and because the convective currents will make the prey distribution dynamic. Taking advantage of these currents would be an energetically favourable way to search the water column. Here, it was shown that riding convective currents in this way gives rise to movement patterns that are well represented by Lévy walk movements with *μ* = 1.25 and so is consistent with some observations of *R. octopus* within shallow (10 m), warm (summertime) waters [8]. Such behaviour may have been observed in the jellyfish, *Aurelia labiate* [52]. Purcell *et al*. [52] observed that these relatively small jellyfish (25–40 cm) all swam vertically in the same direction, either up or down, and speculated that this was because these jellyfish were following convective currents in features like Langmuir circulation cells. This speculation is supported by observations of swimming jellyfish being concentrated in convergences between Langmuir cells [53,54].

Antecedents of this new result can be found in the work of Solomon *et al*. [55], who reported on evidence for Lévy walks in numerical simulations of time-periodic Rayleigh–Bénard convection; a laminar flow that is a precursor to turbulent convective flows [56]. This system admits a Hamiltonian description of the kind that provides a generic mechanism for the production of Lévy walks [57], which was realized, for example, in the controlled laboratory-scale experiments of Solomon *et al*. [57]. Lévy walk characteristics arise when the time-periodicity is tuned so as to produce a ‘resonance’. In the simulations of Solomon *et al*. [55] such synchronous lateral oscillations cause long-range transport across many cells and this causes the diffusivity to diverge. The claim by Solomon *et al*. [55] for resonant Lévy walks rests on the fact that the step-length distribution at the first resonance looks linear over about one decade when plotted on log–log scales. This claim stands up to scrutiny as the Akaike information criterion convincingly favours a power-law distribution over an exponential distribution as a model of the simulation data (results not shown). The maximum-likelihood estimate for the power-law (Lévy) exponent is 1.9. Lévy walks were also found at other resonant frequencies. Resonant Lévy walks are, however, unlikely to arise in naturally occurring time-periodic oscillatory flows, such as shallow tidal systems [58], because these flows will not generally be tuned for resonance. Resonant frequencies appear to be most prevalent (but still scarce) with no-slip (rigid) boundary-conditions rather with free-slip (stress-free) boundary conditions [59].

The new results reported on herein suggest that the notion of optimal Lévy walk searching which underpins the Lévy flight foraging hypothesis can be extended well beyond its original conception in terms of actively searching for sparsely and random targets. It also vividly illustrates that Lévy walks are not necessarily the result of selection pressures for advantageous searching behaviour but can arise freely and naturally from simple processes. Both of these theoretical predictions could be tested by tracking jellyfish while simultaneously measuring the prey field, e.g. with plankton pumps. It would also be interesting to test whether honeybees and bumblebees use fast simulated annealing when searching at the landscape scale for resources. This could be done in harmonic radar studies which to date have been largely focused on searching for single targets [2,3]. The predicted emergence of Lévy walks in convective flows suggests that the Lévy family is vastly larger than previously thought and includes pollens, together with small seeds, mites, spiders and young larvae of certain insects that on warm days disperse passively within the atmospheric boundary layer. In this regard, it is interesting to note that Robert Brown came close to discovering Lévy walks. Brown reported that minuscule pollen particles suspended in nominally *still* water have seemingly random movements. Einstein's [60,61] subsequent mathematical description of these random ‘Brownian’ movement patterns now lies at the heart of the ‘correlated random walk paradigm’—the dominant conceptual framework for the modelling of animal movement patterns [62]. If the water used in the experiments of Brown had been heated from below or cooled at the top surface, then Brown would have observed Lévy walks rather than Brownian walks, and the history of Lévy walk research may have been advanced by many decades.

## Funding statement

This research is funded by the Biotechnology and Biological Sciences Research Council (BBSRC) and by the Australian Research Council.

## Acknowledgements

I thank Stephan Wolf for stimulating discussions about the foraging behaviours of honeybees and bumblebees.

- Received June 24, 2014.
- Accepted July 16, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.