## Abstract

‘Optimal’ behaviour in a biological system is not simply that which maximizes a mean, or temporally and spatially averaged, fitness function. Rather, population dynamics and demographic and environmental stochasticity are fundamental evolutionary ingredients. Here, we revisit the problem of optimal foraging, where some recent studies claim that organisms should forage according to Lévy walks. We show that, in an ecological scenario dominated by uncertainty and high mortality, Lévy walks can indeed be evolutionarily favourable. However, this conclusion is dependent on the definition of efficiency and the details of the simulations. We analyse measures of efficiency that incorporate population-level characteristics, such as variance, superdiffusivity and heavy tails, and compare the results with those generated by simple maximizing of the average encounter rate. These results have implications on stochastic search problems in general, and also on computational models of evolutionary optima.

## 1. Introduction

The question of how an organism should forage to find its food is fundamental and apparently accessible, yet remains a contentious and unresolved issue both in theoretical ecology (e.g. Codling *et al.* 2008; Plank & James 2008; Sims *et al.* 2008; James *et al.* in press) and in the interpretation of movement data (e.g. Edwards *et al.* 2007; Sims *et al.* 2007; Edwards 2008). Recent empirical studies claim to observe animal movements obeying the so-called Lévy walks defined as random walks, where the distance *l* between successive changes in direction obeys a probabilistic power law, *P*(*l*) ∼ *l*^{−μ}, 1 < *μ* ≤ 3. These studies are supported by analytical results which, for certain idealized foraging scenarios, predict optimality for Lévy walks where the exponent *μ* = 2 (Viswanathan *et al.* 1999, 2000; Bartumeus *et al.* 2002; Sims *et al.* 2008); but see also Benhamou (2007) and James *et al.* (2008) for conflicting conclusions.

Central to all the above arguments of optimality has been the notion of search (or foraging) efficiency, defined as the expected number of encounters per unit time, or the expected time to first encounter (James *et al.* in press). In such studies it is assumed, often implicitly, that the best foraging strategy is one which maximizes this foraging efficiency. In terms of evolution by natural selection, this maximization principle may be an oversimplification (Darwin 1859; Real 1980). Evolution favours strategies that maximize the probability of being a winner in Darwin's ‘struggle for existence’, and not those that simply make the organism ‘good on average’ (Yoshimura & Shields 1987; Yoshimura & Clark 1991). Currey *et al.* (2007) illustrate this concept in the context of bone mineral content, showing that, in an evolutionary trade-off between bone stiffness (deterministic) and toughness (stochastic), evolution favours over-investment in the deterministic property; this over-investment can be quantified by taking into account selection pressure, stochasticity in the map from genotype to phenotype, and the nature of the trade-off.

The foraging of planktonic fish larvae for patchily distributed copepod prey in a turbulent (and therefore stochastic) environment forms the motivation for this study. This process underpins the parameterization of stock–recruitment relationships used to inform and develop sustainable fisheries management. Ecologically, the system involves large numbers and massive mortality: a mature female cod will produce millions of eggs per year, only a few of which ever survive to adulthood because of high levels of background mortality (of the order of 5% per day (Chambers & Trippel 1997)). To be successful, a larva therefore needs to feed and grow rapidly to escape this mortality (for example, to reach metamorphosis to the next life-history stage). In this sense, the mean fish is a dead fish; only those larvae that are able to grow the fastest are likely to contribute to recruitment in the next generation (Pitchford & Brindley 2001; Pitchford *et al.* 2005). We use this ecological system as a basis for the models and simulations in this study. We take the term ‘recruitment’ to mean an individual successfully finding a threshold amount of prey items without being eaten by a higher predator. The faster an individual reaches this threshold, the less time it spends exposed to the background mortality, i.e. the larger its recruitment probability.

Here we show how, for foragers with only very limited information, a Lévy foraging strategy can be evolutionarily desirable when compared with a simple random walk. This is not because Lévy motion maximizes the expected encounter rates, but rather because its superdiffusive character results in a large spread in the evolutionarily relevant metric; in this case, a stochastic hitting time relating to recruitment probability. We show that mean-field and diffusive growth approximations are insufficient to describe the behaviour at the population level. Considering the details of our evolutionary model, we argue that, although the Lévy strategy may be optimal, a homogeneous population of Lévy foragers does not necessarily emerge; dimorphism is possible. Although based on the evolutionary foraging problem, this research has wide implications, both on the general problem of random searching with incomplete information and heterogeneous environments and on evolutionary models and their computational implementation.

## 2. Material and Methods

The methods are based on the synthesis of rigorous individual-based simulations with analytical results based on solving hitting time problems for simple diffusion-with-drift stochastic differential equations.

We simulate a cruise predator foraging destructively for identical discrete prey items in a three-dimensional environment. The key output of the simulation is the time, *t*_{hit}, taken by an individual forager to successfully find a predetermined threshold number *N* of prey items (mathematically, the random variable *t*_{hit} is the ‘hitting time’). Each simulation starts by generating a fixed prey field, according to a specified statistical distribution, and placing a forager within it. The forager chooses a direction to move in, with all directions equally likely, and then moves in a straight line of length *l* at a fixed speed, where *l* is drawn from a specified probability distribution. If the predator moves for *l* without encountering a prey item, it then chooses a new direction, and new distance *l* to move in that new direction. If the forager perceives a prey, i.e. if the predator passes within a perceptive radius *r*_{p} of a prey item, then the predator moves to the prey item and (destructively) consumes it. It then chooses a new random direction and new step length *l* to move in that new direction. Once *N* prey are consumed, the simulation stops and the individual's hitting time *t*_{hit} (the total time elapsed) is recorded. This type of forager movement has strong connections with the zero-information (i.e. each new search direction is chosen independently of any previous experience) Lévy foraging models of Viswanathan *et al.* (2000) translated into three dimensions.

We then use these individual-based data to assess recruitment probability, defined as the probability that the larva escapes external mortality for the complete duration of *t*_{hit}. Assuming mortality to occur as a Poisson process with rate *m* then an individual larva's probability of recruitment is exp(−*mt*_{hit}). It is this nonlinear filtering of individual-based hitting times that results in the most rapidly growing larvae being strongly favoured; only the left-hand tail of the *t*_{hit} distribution contributes significantly to recruitment (Pitchford *et al.* 2005). Previous work assumed that the *t*_{hit} distribution arose from an underlying Gaussian process, but the results presented here show that, for a Lévy forager, this assumption is misleading, with important evolutionary consequences.

The three elements of the computer simulations are the environment, the prey field and the forager itself. The environment is a three-dimensional volume of continuous space, measuring 100 arbitrary spatial units along each side. Three types of boundary condition are considered: absorbing, periodic and reflecting, all of which have been previously used in the literature. A predator reaching an absorbing boundary stops, then chooses a new random movement direction (conditional on this being within the environment) and a new movement of length *l*. A predator reaching a periodic boundary is continued in the same direction from the opposite boundary with its total movement length *l* unchanged. A predator reaching a reflecting boundary is perfectly reflected back into the environment, with its total movement length *l* unchanged.

The prey field consists of 1000 prey randomly distributed either uniformly randomly within the environment or distributed patchily. To distribute the prey uniformly, we place a prey item at coordinates (*x,y,z*) where *x*, *y* and *z* are independently chosen from a uniform distribution over [−50,50]. To distribute the prey patchily, a prey item is placed at the end of each step of a Lévy walk, with parameters 1.5, 2.0 and 2.5 for light, medium and heavily patchy prey fields, respectively, and with a random starting position. Examples of uniform and heavily patchy prey fields are depicted in figure 1. See appendix A.1 for the method employed to generate random numbers from a Lévy distribution, and also from a bounded Lévy-like distribution (henceforth called bounded Lévy).

The forager has a fixed speed of 1 spatial unit per time unit, a perception radius *r*_{p} of 1 spatial unit, chooses its new direction uniformly and is required to eat *N* = 50 prey. The distributions from which it can select its time between reorientations are: a delta distribution (centred at *δ* = 2^{n}, *n* = 0, … ,7) simulating a fixed-step random walk; a bounded Lévy distribution (bounded above by 100 spatial units); and a standard Lévy distribution (both with parameters *μ* in the range 1.1, 1.2, … , 3.0). Again, see appendix A.1 for methodological details on random number generation.

The combination of model parameters presented above allows us to generate a large number of independent hitting times over the foraging strategies in a computationally manageable time. Specifically, these options provide 576 different scenarios, each of which was run for 10 000 individuals (generating *M* = 10 000 independent hitting times *t*_{hit}), to arrive at the numerical results summarized in figure 2 and §3.1. These hitting times are evaluated in an evolutionary context in figure 3 and §3.2 using mean-field methods, diffusion approximations and direct calculation. Using recruitment probability and selection pressure as metrics of fitness or evolutionary success, we thereby assess which foraging strategy, involving essentially identical predators moving at identical speeds, should be regarded as optimal in a given environment.

## 3. Results

### 3.1. Hitting time analysis

The foraging process as modelled is highly stochastic; for our default parameter values calculated hitting times *t*_{hit} ranged from 600 to 160 000 time units. The choice of environmental boundary (absorbing, periodic or reflecting) has no qualitative effect on the data, and very little quantitative effect (data not shown). Therefore, results only for the absorbing boundary are shown here. Similarly, results for fixed-step and unbounded Lévy walks only are shown here; the bounded Lévy walk results are qualitatively identical to those of unbounded Lévy walks, but show minor quantitative differences that serve to strengthen the argument in favour of Lévy foraging (see below). The choice of unbounded Lévy walks and absorbing boundary remove any artefacts caused by extremely large move lengths *l* (James *et al.* in press), and in effect *l* will always be truncated by the finite absorbing boundaries in this case. For other boundaries, the bounded Lévy walk can be used to avoid these complications. Further, we present only those results from the uniform and heavily patchy prey field because, as anticipated, the light and medium patchy prey fields generate results that lie between these two extremes.

Figure 2 summarizes these hitting time data. We plot the mean (black), standard deviation (red, with crosses) and 10.0, 1.0, 0.1 per cent quantiles (i.e. the most successful 1000, 100 and 10 individuals, with a dashed blue line, where success is defined as having a small hitting time, *t*_{hit}) for simulations using Lévy foraging (figure 2*c,d*), and also for fixed-step random walks (figure 2*a,b*), in uniform (figure 2*a,c*) and heavily patchy (figure 2*b,d*) prey fields. Figure 2 shows that if, as in previous studies, one were to equate optimal foraging with minimum mean hitting time then choosing Lévy parameter *μ* = 1.1 in figure 2*c* and *μ* = 1.4 in figure 2*d*, or a fixed step length of *δ* = 64 in figure 2*a* and *δ* = 128 in figure 2*b*, would be optimal. The minimum mean hitting times are 16 870, 16 785, 17 092 and 16 150 for figure 2*a*–*d*, respectively. A ballistic strategy emerges as *δ* → ∞ and as *μ* → 1, and these plots emphasize this link between the two foraging strategies under consideration. An optimal foraging argument based on mean efficiency would therefore conclude that: in a uniform prey environment, the choice of a movement strategy is relatively unimportant but that large fixed-step random walks may be marginally favoured; in a patchy environment, a Lévy strategy with *μ* = 1.4 is favoured, but again with only a weak selective advantage attached to the putative optimal exponent as fixed-step and Lévy strategies with parameters *δ* = 64 128 and *μ* = 1.1–1.8 all lie within 5 per cent of this optimum.

### 3.2. Probability of recruitment

The central argument of this research is that, from an evolutionary point of view, mean foraging efficiency can be misleading or irrelevant. Figure 3 illustrates this fact, by calculating the probability of recruitment (our first metric for the success of a foraging strategy) on the vertical axis for various foraging scenarios on the horizontal axis. As in figure 2, the left-hand panels show results for a fixed-step random walk and the right-hand panels for a Lévy forager. The data summarized in red apply to a heavily patchy environment, those in blue (marked with crosses) correspond to a uniform environment; the results for the Lévy forager (figure 3*d*–*f*) are discussed first.

We assume in all cases that, while it is foraging, the forager is exposed to a constant rate of external mortality *m*, specifically that individual mortality occurs as a Poisson process with rate *m* = 0.25 × 10^{−3} d^{−1}. This figure is broadly representative of recruitment probabilities of *O* (1%) observed in marine fish species (Chambers & Trippel 1997). Given this mortality rate, the probability of recruitment of an individual could be estimated in many ways, and three possible estimates are explained and quantified below.

Figure 3*d* uses a simplistic mean-field approach. Defining *T* to be the mean hitting time for each scenario (*μ* = 1.1, … , 3.0) (i.e. *T* = *E*(*t*_{hit}) = ∑*t*_{hit}/*M*), one could estimate the probability of recruitment by
3.1
As it is based entirely on mean efficiency, this calculation yields an optimum strategy, in the heavily patchy environment, of *μ* = 1.4 with a probability of recruitment of 0.0176. The value *μ* = 1.4 is to be expected as this measure is a simple transform of the mean-field analysis found in §3.1. In a uniform environment, the values of *μ* close to 1.1 are favoured, again as anticipated from the analysis above. However, we argue that these estimates can be misleading; they undervalue the contribution of the smaller *t*_{hit} observations that have higher probabilities of recruitment, and similarly overvalues the contribution of the largest *t*_{hit} observations (Pitchford *et al.* 2005).

It is possible to improve the accuracy of this calculation by taking into account the variability in the hitting times. Using the method of (Pitchford *et al.* 2005), which assumes that each individual's growth arises as a solution to a diffusion-with-drift Gaussian process (i.e. that growth is diffusive in character), we define the metric used in figure 3*e* as
3.2
where the drift, *r*, and diffusion, *σ*, parameters are given by
with *V*(*t*_{hit}) = ∑(*t*_{hit} − *T*)^{2}/*M*. Recruitment probabilities in the heavily patchy environment, where there is larger variability in hitting times, are noticeably increased. For a Lévy forager, recruitment probability increases by around 80 per cent, with an optimal Lévy exponent of *μ* = 1.8 and probability of recruitment of 0.0318. In the uniform environments, values of *μ* close to 1 are favoured with *μ* = 1.1 having a probability of 0.0165.

The assumption underlying equation (3.2) is that individual growth trajectories can be characterized by diffusive processes. This assumption does not necessarily hold, especially when a forager moves with Lévy-like movements in a patchy environment (Burrow *et al.* 2008); the distribution of sizes at a given time is characteristically heavy-tailed. Figure 3*f* uses all of the simulated hitting time results, rather than making an implicit diffusive assumption, to estimate the probability of recruitment as
3.3
where *f* is the underlying probability distribution for each foraging strategy. In this definition we have, by necessity, had to apply a numerical estimation from the simulated data. In the patchy environments, the maximal recruitment probability is increased twofold over the mean-field estimate (3.1) with the maximum probability of recruitment, at *μ* = 1.8, being 0.035 (figure 3*f*). The values *μ* = 1.5–1.9 lie within 5 per cent of this maximum and those between *μ* = 1.4 and 2.3 lie within 10 per cent. As discussed below, such Lévy foragers comfortably outperform fixed-step strategies in patchy environments, but not in uniform environments. Estimates of evolutionary fitness based on mean-field arguments such as equation (3.1) can be both qualitatively and quantitatively misleading. This maxim is applicable in a wider range of biological situations, for example Codling (2008) and Pitchford *et al.* (2007) make the same argument in the analysis of a marine reserve-fishery system.

We apply the same definitions of *p*_{mean}, *p*_{stoch} and *p*_{dist} to the fixed-step random walk scenarios in figure 3*a*–*c*. The probabilities of recruitment for fixed-step walks in uniform prey fields (blue lines in figure 3*a*–*c*) maintain the same shape and show a slight increase in probability of recruitment from *p*_{mean} to *p*_{stoch} to *p*_{dist}. In a uniform prey environment, fixed-step and Lévy foragers (figure 3*d*–*f*) both increase their recruitment probability as they tend towards ballistic motion (having larger fixed steps or *μ* → 1). The *p*_{mean} and *p*_{stoch} metrics for fixed-step random walk in the patchy environment indicate that close-to-ballistic motion is optimal, whereas *p*_{dist} shows that a locally diffusive (short step) strategy is optimal. However, all of these fixed-step strategies in patchy environments are out-performed by Lévy strategies.

Overall, the results presented in figure 3 lead to the conclusion that Lévy foraging strategies (with *μ* ≈ 1.8) are optimal in patchy prey environments, whereas ballistic motions (large fixed-step random walks) are favoured in uniform environments. The quantitative advantage of Lévy foraging in a patchy environment is potentially large, but this advantage becomes apparent only when the full distributional data (and not averages) are accounted for.

### 3.3. Selection pressure

In this section we analyse an alternative mortality model that operates at the population level and not on each phenotype or foraging strategy individually. We use selection pressure *p* to represent mortality in the model, defined to be the fraction or percentage of a population to survive based on a fitness criterion. This definition has a stronger biological foundation, through the ‘survival of the fittest’ paradigm, and allows for more complicated conclusions to be drawn; for example, population-level optimality may be defined by a distribution of strategies rather than a single strategy.

We create a heterogeneous population of 280 000 individuals over 28 phenotypes by amalgamating our simulated hitting time results for fixed-step and unbounded Lévy foragers. Each of the 28 foraging strategies (Lévy with *μ* = 1.1, … , 3.0 or fixed step with *δ* = 2^{n}, *n* = 0, … , 7) is represented by 10 000 independent members of the population. Note that each individual strategy contributes approximately 3.6 per cent of the initial population and all of the Lévy strategies together contribute 71.4 per cent. Below, we compare the performances of this population in both the uniform and heavily patchy environments with absorbing boundaries. Figure 4 describes the surviving population structure under selection pressures of *p* = 10 per cent (figure 4*a*,*d*), *p* = 1 per cent (figure 4*b*,*e*) and *p* = 0.1 per cent (figure 4*c*,*f*). The largest hitting times in each surviving population (i.e. the worst performing survivor) are 10 473, 5791, 2940, 15 199, 12 729 and 11 226 for figure 4*a*–*f*), respectively. We first consider the population in the heavily patchy environment (figure 4*a*–*c*).

At a 10 per cent level, the Lévy walks with *μ* between 1.4 and 2.3 each contribute over 4 per cent surviving population, with the peak at *μ* = 1.8 contributing 1306 out of 28 000 (4.7%) and all of the Lévy walks together constituting 80 per cent of the surviving population. As the selection pressure increases from 10 per cent through 1 per cent to 0.1 per cent, the surviving Lévy walk members tend towards diffusive/Brownian motion (*μ* → 3.0) and, at the extreme level, 0.1 per cent, the Lévy walks form a smaller proportion of the surviving population (70%). In contrast, the fixed step-length walks show a tendency towards short step lengths as selection pressure increases. Under extreme selection pressure, 0.1 per cent, fixed step-length walks contribute 30 per cent of the surviving population with *δ* = 1 and 2 dominating with 65 out of 280 (23%).

In other words, the level of selection pressure may dictate the eventual evolutionarily favourable outcome, as in Currey *et al.* (2007). If the selection pressure is extreme, then the optimal strategy, in a patchy environment, may be to use a fixed-step random walk with very small step lengths (the best strategy is to hope to initially be in a dense prey patch, and then not to leave it). Considering the data presented in figure 4, such fixed-step strategies would be very unlikely to be eradicated from the population, resulting in a mixed population of Lévy and fixed-step foragers over evolutionary time. In the example of fish recruitment, the observed highly variable stock-recruitment data (e.g. Minto *et al.* 2008) suggest that selection pressure on larvae and juveniles varies stochastically from year to year, which could make the evolution of single strategy foraging still less likely.

In the uniform prey field environment, figure 4*d*–*f* shows that ballistic motion (*δ* → ∞ and *μ* → 1) forms the largest proportion of the surviving populations for each selection pressure. Considering the quarter of the strategies that are closest to ballistic motion, i.e. *δ* = 64 128 and *μ* = 1.1, … , 1.5, we notice that the proportion of the surviving population with these strategies increases as selection pressure increases, from 14 236 out of 28 000 (51%) through 1590 out of 2800 (57%) to 181 out of 280 (65%) when *p* = 0.1 per cent. The conclusion of ballistic motion being consistently optimal irrespective of the selection pressure is in contrast with the patchy prey field, where different strategies are optimal at different pressures. Again, we note that, for all of these experiments (figure 4*a*–*f*), the numerical evidence supports the theoretical expectation that large *δ* and *μ* close to 1 will exhibit similar properties as both approximate ballistic motion.

This method of choosing a surviving population strongly echoes the selection step in a genetic algorithm. A genetic algorithm is a biologically inspired evolutionary search algorithm that, at each generation, selects a proportion of a population to breed based on fitness criteria and then these survivors combine and mutate to form the successive generation. A key feature of genetic algorithms is the stochasticity that can be applied in the selection, combination and/or mutation functions. Our model captures this stochasticity in the fitness criteria through the generation of hitting times by the simulator. We show only the survivors from the first generation (i.e. without any combinations or mutations applied) in figure 4 and can only postulate about the long-term behaviour of any particular genetic algorithm.

## 4. Discussion

The results presented above show that Lévy searching can be evolutionarily favoured for a stochastic forager in a high mortality, patchy environment where the key fitness criterion is the amount of time taken to encounter a fixed number of objects. They also show that fixed-step foragers can be evolutionarily favoured by including a different performance metric in the definition of fitness. These alternate conclusions of ‘optimality’ both arise not because such foragers are best on average (typically they are not), but because of the emergent heavy-tailed nature of the hitting time distributions, specifically the distribution of *t*_{hit} as *t*_{hit} → 0. This raises important questions about the notion of evolutionary optimality, and models thereof.

The existence of an external mortality effect is central to the results in §§3.2 and 3.3; otherwise, there is no evolutionary incentive for any foraging strategy. Our choice of the rate parameter *m* in §3.2 yields recruitment probabilities of the correct order of magnitude for fish, but the results are qualitatively robust. Larger rates of external mortality amplify the importance of the tails of the hitting time distributions, increasing the difference between average foraging efficiency and recruitment probability. The choice of the selection pressure at the population level in §3.3 represents the mortality by a more direct method. In this model, the effect of setting *p* = 1 per cent, the same order of magnitude of survival as in §3.2, leads to different conclusions about the optimal foraging strategy at a population level.

Our model concentrates on the question of when a fixed speed (equivalent to a fixed energy cost per unit time) forager should choose to reorientate itself. In reality, the forager can control its speed (energy cost) as well as its direction. This yields an interesting trade-off reminiscent of Currey *et al.* (2007)—faster swimming incurs an increased deterministic energy cost to overcome drag, but results in a stochastic benefit (increased stochastic encounter rate)—and the Currey *et al.* (2007) model predicts an evolutionary over-investment in minimizing the deterministic cost. It is worth mentioning that, although the optimal foraging predictions of Pitchford *et al.* (2003) are qualitatively robust, they suggest optimal swimming speeds at the upper end of observations. A synthesis of the evolutionary models of Currey *et al.* (2007) with the three-rate probabilistic model of Pitchford *et al.* (2003) will result in reduced optimal swimming speeds, and forms a focus of ongoing work.

### 4.1. Distributional data

Throughout this paper we advocate using the entire distributional data from individual strategies or populations to draw conclusions about optimality. In figure 5, we show a selection of the hitting time distributions that give rise to the analysis above. Figure 5*a* shows the probability distribution functions (PDFs) for three key foraging strategies in the heavily patchy environment, namely *δ* = 1, *μ* = 2 and *μ* = 3. These three distributions clearly provide the full picture for the appropriate mean, quantile and standard deviation information summarized in figure 2*c*,*d*. Figure 5*b* shows the comparison of cumulative distribution functions (CDFs) between the same three distributions for low hitting times (i.e. the best-performing foragers), and it is the distributions of these best-performing foragers that provides greatest insight. It is also worth noting that foraging using the Lévy step-length distribution, which itself has a heavy tail, produces the hitting time distribution with the smaller, less heavy tails when compared with the hitting time distributions generated by the fixed-step distributions (i.e. by comparing *μ* = 2 with *δ* = 1 in figure 5*a*).

The changes in probability of recruitment from *p*_{mean} (equation (3.1)) to *p*_{dist} (equation (3.3)), i.e. from figure 3*a*,*d* to 3*c*,*f*, by including an exponential reward for performance has two effects: (i) an increase in absolute probability from *p*_{mean} to *p*_{dist} and (ii) an increase in the relative performance of the fixed-step walks in comparison with Lévy walks. It is the significant increase in relative performance of the fixed-step foragers that first indicates that distributional metrics may confer an advantage over mean-field analysis. The relative performance of *δ* = 1 compared with *μ* = 2 from *p*_{mean} to *p*_{dist} rises from 0.007 to 0.51. By disproportionately (via the exponential mortality rate) rewarding the fastest foragers, we see that fixed-step foragers, with short step lengths, become viable possibilities for optimal strategies and that is purely because of the heavier tail of the best-performing foragers (figure 5*b*). This advantage for fixed-step foragers, of heavier low hitting time tails, in the probability of recruitment analysis cannot overcome the disadvantage of heavier high hitting time tails and, in §3.2, leads to the conclusion that Lévy foragers outperform fixed-step strategies.

The advantage of heavy tails of low hitting times is emphasized when using the concept of selection pressure to determine relative performance of foraging strategies. The two vertical dotted lines in figure 5*b* show the 0.1 per cent and 1 per cent threshold hitting times for the populations in figure 4*a*–*c*. By examining the cumulative distributional data, we see the patterns from the selection pressure analysis in more detail. We see the dominance of *δ* = 1 and *μ* = 3 at the highest selection pressure (figure 4*c*), the relatively stronger performance of *μ* = 2 when *p* = 1 per cent and *μ* = 2, having higher cumulative density than *δ* = 1 between the 1 per cent and 10 per cent levels. The selection pressures at which each strategy outperforms each of the others (if they ever do) are clearly available from the cumulative density data and so more detailed analysis may be undertaken, such as for heterogeneous population demographics or for evolutionarily stable strategies (as in Currey *et al.* (2007)). In figure 5*b* *δ* = 1 outperforms *μ* = 3 when *p* ≤ 0.7 per cent (threshold at 5239) and outperforms *μ* = 2 when *p* ≤ 2.3 per cent (7174). Although not shown in figure 5*b*, *μ* = 2 eventually outperforms *μ* = 3, when *p* ≥ 5 per cent (8766).

Hence, a key differential in determining optimal foraging strategies is the high performance tail of the distribution; having a heavy tail of best-performing foragers can provide a significant advantage at both the individual and populations levels. The question that now arises is how to incorporate the ‘heaviness’ of the tails of each distribution into a relevant metric, i.e. is *δ* = 1 heavier than *μ* = 3 because it has an advantage when *p* ≤ 0.7 per cent or is *μ* = 3 heavier because it has an advantage for the significantly larger range *p* ≥ 0.7 per cent? We have used the distributional data in the performance metrics above, through the probability of recruitment (equation (3.3)) or selection pressure, to show that optimality in biological situations is not necessarily captured by mean-field analysis alone.

### 4.2. Model limitations

The foragers in our simulations employed perfect, spherical perception enabling them to sense every prey item within their perception range of 1 spatial unit. In the ocean there is a plethora of environmental cues, as well as imperfect and directional perception, informing each movement decision a forager makes. Factors such as light levels, water density, salinity and turbulence, effect perception and efficiency of movement are not considered here.

The destructive nature of the foraging changes the underlying prey field; at every foraging event, the prey distribution changes on both a local and global level. Parallels may be drawn with chemical reaction kinematics where reaction rates fall with time due, in part, to the heterogeneity of the mixing of the reactants (Kopelman 1988). The relatively low proportion of prey consumed (*N* = 50, 5%) alleviates this problem in our simulations; the original patchy prey distribution changes slightly (but is not destroyed) and relevant hitting time data are generated for multiple predator–prey encounters. Other complex interactions are also beyond the scope of this basic model, but are amenable to more specific mathematical models. Accounting for turbulence and prey movement as well as consumption by multiple foragers leads to ever more complex and stochastic time-dependent prey fields (Pitchford & Brindley 2001), while environmental factors can lead to emergent prey field topologies, for example (Durham *et al.* 2009). Considering coevolutionary aspects such as group- versus individual-based foraging (Wood & Ackland 2007) and infochemical signalling across trophic levels (Steinke *et al.* 2006) may also be important. In general, there is no universality in the prey fields and ecological interactions. This weakens the argument that Lévy foraging should be universally optimal purely because of its alleged high mean efficiency in ideal environments. It is still possible for highly stochastic (Lévy-like) movement strategies to be evolutionarily favoured in complex and dynamic environments, but this will rely on them having a heavy tail at the high performance end of the distribition.

As each forager lives in a combination of sensory gradients and emergent prey fields, its perception of the direction to, and proximity of, prey exists at multiple scales, possibly leading to composite foraging behaviour. This reliance on the underlying environmental factors in determining multiple behavioural patterns informs the current debate of pattern versus process and the relative efficiency of composite random walks (Benhamou 2007; Plank & James 2008). A key extension to the analysis above is the possibility to investigate emergent patterns arising from the processes that will inform the practical foraging research; such as how do fractal prey fields arise and what are the emergent distributions of interforaging event times?

It should be noted that we have based our models on larval fish foraging but have deliberately simplified them to emphasize the new methods of data analysis for stochastic search problems that we advocate above. We have minimized the number of model parameters, as much as possible, to keep the effects of the model manageable (bearing in mind the pattern versus process debate) and to keep the content of this study focused on the analysis of the hitting times. Our qualitative conclusions are also robust to changes to the constants in the model, i.e. *N*, *r*_{p} and the size of the environment. The results are also qualitatively unchanged when fixed-step random walks are replaced by exponentially distributed step lengths with the same mean. By selecting Lévy walks and fixed step-length movement strategies, we provide two contrasting foraging techniques that keep the analysis relevant to current optimality research (in the Lévy tradition of Viswanathan *et al.* (2000)) and give an easily accessible alternative through a classic diffusive search pattern. One key simplification is the selection of a new search direction for each movement step that is independent of any previous movement direction, and this leads to the models being classed as zero-information searches. Although there is experimental evidence to support larval foragers exhibiting rudimentary cognitive ability (Dytham & Simpson 2007), the zero-information assumption is commonly applied in foraging simulations (Zollner & Lima 1999; Viswanathan *et al.* 2000; Benhamou 2007) and we continue with that assumption here.

The model as defined above, with its limitations, will affect our conclusions of optimality under our two metrics: probability of recruitment and selection pressure. Indeed, as mentioned above, the superior performance of the fixed-step random walk with *δ* = 1 under extreme selection pressure (figure 4*c*) could be regarded as an artefact of experiment, essentially corresponding to a forager being lucky enough to start in, or close to, a patch and then staying in, or around, it to gain enough prey from it.

## 5. Conclusions

We have shown that the mean efficiency is a misleading metric to use in evolutionary optimal foraging problems, principally because it ignores population dynamics and demographic and environmental stochasticity. This notion is not new (Real 1980), but its consequences have not been accounted for in recent models of search behaviour. Evolutionary conclusions should be drawn from the distributional behaviour of the significant section or the entire population, for example using equation (3.3) or selection pressure (§3.3) techniques. These methods reward the strategies that have the best-performing individuals, i.e. those achieving low hitting times in the extreme left of the population distribution. Mean (3.1) or diffusive probabilistic (3.2) approximations can ignore or underestimate important effects generated by the tails of the distributions. Using selection pressure in §3.3 to compare strategies in direct competition with each other in heterogeneous populations shows that either search strategy type (Lévy and fixed step) could dominate, depending on the intensity of the pressure, and this method allows for the possibility of strategy dimorphism.

One conclusion of our work is that Lévy foragers can be evolutionarily optimal in patchy environments with high mortality; the heavier high performance tails (low hitting times) of their distributions can confer an evolutionary advantage by allowing larger numbers of lucky individuals to succeed. Equally important, we have also shown that the optimality of Lévy foragers is not a universal conclusion; other strategies that produce heavier tails may have an evolutionary advantage, e.g. the case of short fixed-step random walkers that perform better for only a small range of extreme selection pressures but have significantly worse performance at every other scale. The problem of how to exploit a stochastic and uncertain environment, leading to probabilistic distributions of realized fitness which translate to reproductive success in future generations, is common in biological systems. We have focused on a foraging problem inspired by planktonic fish larvae but the underlying message is more generic: evolution rewards those with a chance to operate at the extremes of performance. How extreme performance is measured, modelled and analysed is fundamental.

## Acknowledgements

We wish to thank all seven referees for their time and invaluable feedback during the review process. M.P. would also like to express his thanks for the Pump Priming Funds from the University of York for providing financial support during this project. A.J.W. research is supported by an RCUK fellowship. Special thanks go to Alex James and Richard Brown at the University of Canterbury, New Zealand, for verifying our numerical results.

## Appendix A

#### A.1. Generating random numbers

Let *U* be a uniform random variable on the range (0, 1). Then we define an unbounded Lévy random variable (*UL*) and a bounded Lévy random variable (*BL*), with Lévy exponents *μ*, as
and
where *l*_{min} and *l*_{max} are the minimum and maximum of the range of the Lévy random variables. We have used the values *l*_{min} = 1 and *l*_{max} = 100 in our simulations.

These algorithms have been used to generate the patchy (fractal-like) prey fields and to simulate predator movements for the Lévy forager.

#### A.2. Supplementary material

Supplementary material, available on the project website (http://www.yorkevolution.org/publications), includes the raw hitting times data for all 5 760 000 simulations and Matlab code to manipulate data and create graphics equivalent to figures 2–5 for other model parameters.

## Footnotes

- Received February 17, 2010.
- Accepted March 4, 2010.

- © 2010 The Royal Society