PRL 119, 250601 (2017)

week ending
22 DECEMBER 2017

PHYSICAL REVIEW LETTERS

Emergence of Lévy Walks from Second-Order Stochastic OptimizationŁukasz Kuśmierz* and Taro Toyoizumi†
RIKEN Brain Science Institute, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan(Received 8 October 2017; published 18 December 2017)
In natural foraging, many organisms seem to perform two different types of motile search: directed
search (taxis) and random search. The former is observed when the environment provides cues to guidemotion towards a target. The latter involves no apparent memory or information processing and can bemathematically modeled by random walks. We show that both types of search can be generated by acommon mechanism in which Lévy flights or Lévy walks emerge from a second-order gradient-based
search with noisy observations. No explicit switching mechanism is required—instead, continuoustransitions between the directed and random motions emerge depending on the Hessian matrix of the costfunction. For a wide range of scenarios, the Lévy tail index is α ¼ 1, consistent with previous observations
in foraging organisms. These results suggest that adopting a second-order optimization method can be auseful strategy to combine efficient features of directed and random search.
DOI: 10.1103/PhysRevLett.119.250601

Many organisms must actively search for resources inorder to survive and produce offspring. The foraging theory
examines the various search strategies implemented byorganisms depending on their abilities and the environments
in which they live. In a directed search, the greater
involvement of sensory and information processing abilities
enable more complicated strategies. In contrast, in theboundary case of a memoryless and senseless forager, theonly option is to wander randomly in the environment(random search). Even in this case, however, different
strategies exist, depending on the character of the random
motion. A natural candidate model for the random strategy is
Brownian motion that describes a wide range of natural
phenomena, including the movement of inanimate particles
under thermal noise. A prominent feature of Brownian
motion is the linear growth of the variance of the position
with time. However, empirical data indicate that for organisms the observed growth is often faster. Lévy walks (LWs)[1–3] and similar Lévy flights (LFs) [2,4,5] have beensuccessfully applied to fit experimental data obtained fromthe movement patterns of many organisms and their cells,
including T cells [6], microglia [7], starved slime mold
(Dictyostelium discoideum) [8,9], swarming bacteria [10],fruit flies [11], honey bees [12,13], wandering albatrosses[14,15], marine predators [16], and humans [17–19] (also ina human’s gaze [20] and word association [21] trajectories).
In many different random search scenarios, LWs and LFs
have been shown to be advantageous over normal diffusion[22–32] and alternative superdiffusive models [33]. Theseobservations have led to the so-called Lévy flight optimal
foraging hypothesis, which states that LFs (or LWs) represent evolutionary adaptations due to their distinct advantages over other random search strategies [22,34,35].
Recently, this view has been disputed, because none of thementioned organisms is senseless and all of them are able to0031-9007=17=119(25)=250601(6)perform some forms of directed search (taxis); for example,
T cells and isolated bacteria perform chemotaxis [36–40],
whereas fruit flies perform phototaxis [41], geotaxis [42],and chemotaxis [43,44]. Indeed, a number of studies haveshown that characteristics of LFs and LWs may emerge
naturally on large scales from more realistic case specificmodels of movement [45], including simple deterministic
and semideterministic walks in complex environments[35,46–50], self-avoiding random walks [51–53], diffusion
with a time-varying diffusion constant [54–56], and amultiplicative, self-accelerating process [9,57,58]. It has
also been suggested that in some species a power-law
distribution of lengths of straight line segments of their
movement patterns, a hallmark of LWs and LFs, is aconsequence of either the Weber-Fechner law in odometry[59], a power-law distribution of switching times between
competing activities [60–65], or a so-called aerial lottery[66–68]. Moreover, in some cases the apparent superdiffusive character of the population dynamics may be an artifactof averaging over an ensemble of the diffusive motions ofindividuals with diverse characteristics [69].
These studies suggest that LWs and LFs naturally arise inmany realistic biological settings, but they do not argue
why an apparent common behavior is observed acrossspecies and environments. Recently, a generalization of the
LF optimal foraging hypothesis was proposed that explicitly combines directed and random search strategies.
Specifically, an ad hoc combination of taxis for choosing
a direction and random, heavy-tailed distributed step
lengths was shown to be efficient under some search
conditions [27,70]. In contrast, here we propose a novel
mechanism by which LWs and LFs can emerge from ageneric, locally optimal, directed search strategy. In ourmodel, the directed search is realized as a taxis driven by
local observations of a cost function (e.g., a repellent

250601-1© 2017 American Physical Society

week ending
22 DECEMBER 2017

PHYSICAL REVIEW LETTERS

PRL 119, 250601 (2017)

concentration minus an attractant concentration) whose
minima correspond to targets. Inspired by the second-order
gradient-based optimization techniques known from computer science, we assume that the search is based on noisy
gradient and Hessian estimates. As we show below, thisgenerically leads to heavy tails of the step distribution. In
contrast to previous models, our model predicts a continuous crossover between random Lévy searches and directed,
deterministic taxis depending on the amount of informationon the target location provided by observations.
This Letter is organized as follows. First, we fix thenotation and introduce a one-dimensional version of ourmodel. Next, we list different scenarios in which we are able
to prove the existence of the heavy tails. We then discuss how
the tails are affected by the landscape and the observation
methods. Finally, we discuss a multidimensional generalization of the model followed by concluding remarks.
Let ðxn Þ∞
n¼0 be a sequence generated by the Newton
optimization rule

where the tail index 0 < α < 2. In the following, we showthat for a wide range of scenarios the random walk definedby (1) and (2) is equivalent to a (possibly inhomogeneous)
LW or LF (depending on how it is mapped into a
continuous time process [71]) with α ¼ 1. We shall first
analyze the case when β ¼ 0 and fðxÞ ¼ const, so that only
noise is sampled. Assuming that the noise is Gaussian andthat both f 0 and f 00 are measured independently and withoutbias, we can writexnþ1 ¼ xn þ Δn ;where γ ¼ ðσ G =σ H Þ. Comparing (6) with (4), we see that in
our case α ¼ 1. More generally, let us assume that thenumerator XG and denominator XH in (2) are independent
random variables with PDFs ρXG and ρXH , respectively.
This is the case if ξG and ξH are conditionally independent
given the current position of the walker. The asymptotic
form of the PDF of Δ ¼ −XG =X H is given by


Z ∞
Z ∞
z1
ρΔ ðzÞ ¼
dz1 ρXG ðz1 Þ
dz2 ρXH ðz2 Þδ z þ
z2
−∞
−∞


Z ∞
1
z
¼ 2
dz1 ρXG ðz1 ÞρXH − 1 jz1 jzz −∞hjXG jiρXH ð0Þ¼
þ oðz−2 Þ;ð7Þz2R
where the last equality holds if hjXG ji≡ ρXG ðzÞjzjdz < ∞and 0 < ρXH ð0Þ ¼ limz→0 ρXH ðzÞ < ∞. The condition
hjXG ji < ∞ is equivalent to the statement that the tailsof ρXG ðzÞ decay faster than z−2 . If this condition is not
fulfilled, the appearance of heavy tails in the distribution of
Δ is trivial. In our case, however, the heavy tails of ρΔ
appear due to a nonzero probability of XH being arbitrarilyclose to zero. The described mechanism is very general, as
it does not assume that the noise distribution has heavy
tails. Intuitively, the division in (2) takes the role of a noise
amplifier. Clearly, first-order methods, such as the steepest
descent (3), do not involve a division by a random variable
and therefore do not generically lead to heavy tails.
In the case of correlated XG and XH , the presence ofheavy tails cannot be ensured, in general. For instance, if
XG ¼ −YXH for some random variable Y, then the resulting Δ has the same distribution as Y. However, as we will
now show, the heavy tails are still present in the genericcase of the bivariate normal distribution of XG and XH :

ð1Þ

with
ðnÞ

Δn ¼ −

f 0 ðxn Þ þ ξG

ðnÞf 00 ðxn Þ þ β þ ξH;

ð2Þ

where the cost function f∶R → R is to be minimized. TheðnÞðnÞ
rule with β ¼ ξG ¼ ξH ¼ 0 performs a gradient descentor a gradient ascent on f, depending on its curvature.
A positive constant β (damping) is added to the denominator in order to turn this algorithm into a minimizer. Notethat the steepest descent method0ðnÞ~ n ¼ − f ðxn Þ þ ξG
Δ
β

ð3Þ

is recovered from (2) in the limit of β → ∞ if f has aðnÞðnÞbounded second derivative. Terms ξG and ξH account fornoise: If the optimization is to be performed in the physicalworld, derivatives of f are based on noisy measurements.
Similarly, in many optimization problems solved on acomputer, especially in machine learning, a function to
be optimized is estimated with finite precision. With thesedefinitions, the sequence ðxn Þ∞
n¼0 denotes a one-dimensional
discrete-time random walk.
Trajectories of LWs consist of linear segments (orinstantaneous jumps in the case of LFs) Δn , which are
independent and identically distributed random variables(hence we omit the time index n in the discussion of
distributions and write simply Δ). The probability density
function (PDF) of Δ is characterized by heavy tails, i.e., forlarge jzjρΔ ðzÞ ∼ jzj−1−α ;ð4ÞξΔ¼− G;
ξH

ð5Þ

where ξG and ξH are independent Gaussian variables withzero mean and standard deviations σ G and σ H . The reader
can easily verify that Δ is in this case characterized by the
Cauchy distributionρΔ ðzÞ ¼

250601-2

1 γ;
π γ 2 þ z2

ð6Þ

PRL 119, 250601 (2017)

PHYSICAL REVIEW LETTERS
jPj1=2
1⊺
ρX ðxÞ ¼
exp − ðx − μÞ Pðx − μÞ ;
2
2π

ð8Þ

where X ¼ ðXXHG Þ is a two-dimensional Gaussian random
vector, μ ¼ ðμμHG Þ is a vector of its expected values, P ¼
ðPP1121 PP1222 Þ is a symmetric, positive-definite precision matrix,
and jPj is itsRdeterminant. The PDF of Δ can be calculated
as ρΔ ðzÞ ¼ d2 xρX ðxÞδ½z þ ðx1 =x2 Þ and in the special
case of μG ¼ μH ¼ 0 simplifies to the shifted Cauchy
distributionρΔ ðzÞ ¼

1
jPj1=2:π P22 − 2P12 z þ P11 z2

ð9Þ

In general, ρΔ takes the formZh 1i IðzÞ
jPj1=2 ∞ρΔ ðzÞ ¼
dxjxj exp − ð~x − μÞ⊺ Pð~x − μÞ ¼ 2 ;22z
2πz −∞
ð10Þxwhere x~ ¼ ð−x=z
Þ. Since 0 < limz→∞ IðzÞ < ∞, we see that
yet again ρΔ ðzÞ ∼ 1=z2 for large z.
It may now seem like the second-order methods should
always lead to LFs or LWs given noisy observations, whichmight prevent them from being an efficient search strategy.
However, this is not the case if the regularization factor β
and curvature f 00 ðxÞ in the denominator strongly temper
heavy tails in (2). For example, in the case of independent
Gaussian ξG and ξH , the large z limit of the cumulative
distribution of the step size in (7) is given by

PðjΔj > zÞ ≈2hjXG jiρξH ½−β − f 00 ðxÞ
2hjX ji
2¼ pﬃﬃﬃﬃﬃﬃ G e−ðc =2Þz
2π σ H z
ð11Þ

with c ¼ (β þ f 00 ðxÞ)=σ H . Namely, heavy tails are stillpresent, but they are suppressed by the exponential factor2
e−c =2 . Thus, for jβ þ f 00 ðxÞj ≫ σ H the probability of large
random displacements is extremely low. Equation (11)provides a hint as to why the noisy second-order search
may be efficient: If β and σ H are chosen such that c ≈ 0 at
the maxima, where f 00 is negative, and c ≫ 1 at the minima,
where f 00 is positive, heavy tails are present in the vicinityof the maxima facilitating fast escapes, whereas around theminima heavy tails are strongly suppressed, allowing for aneffective local exploration.
We now address the question of how the method of
estimating f 0 and f 00 from noisy measurements of f can
influence our results. The simplest possible model in 1D
consists of three observations. Let us assume that theobservations are performed at x0 − δx, x0 , and x0 þ δx,
resulting in the following noisy measurements: y− ¼fðx0 −δxÞþξ− , y0 ¼ fðx0 Þþξ0 , and yþ ¼ fðx0 þδxÞþξþ ,where ξ• represent multivariate Gaussian noise. If weassume that δx is small enough, we can write the following

week ending
22 DECEMBER 2017

formulas for the maximum likelihood estimates of the firsttwo derivatives:
yþ − y−≈ f 0 ðx0 Þ þ ξG ;2δxy þ yþ − 2y0≈ f 00 ðx0 Þ þ ξH ;f̂ 00 ðx0 Þ ¼ −δx2
f̂ 0 ðx0 Þ ¼

ð12Þ

where ξG ¼ ðξþ −ξ− Þ=2δx and ξH ¼ ðξþ − 2ξ0 þ ξ− Þ=δx2 .
Hence, ξG and ξH are two jointly Gaussian randomvariables. As shown in (9), this generally yields the LW
or LF with α ¼ 1. This reasoning is still valid in scenarioswith more measurements, whenever the desired estimatesare based on linear combinations of noisy observations.
Note that the more measurements used in the estimators,
the better the Gaussian model of noise.
Finally, we turn our attention to the multidimensionalcase. For simplicity, we assume that the search space is RD
with D ∈ N. The jump vector, in analogy to (2), takes thefollowing form (as before, since we focus on a single step,
we omit the step index):
Δ ¼ −A−1 ∇f̂ðxÞ ¼ −½HfðxÞ þ βI þ ξH −1 ½∇fðxÞ þ ξG ;
ð13Þ
where HfðxÞ denotes the Hessian of f, ξG is a noise vector,and ξH is a symmetric noise matrix. Under mild conditions,in the limit of D → ∞ the noise eigenvalues λðξH Þ follow
the Wigner semicircle distribution [72–75]. If the curvatureand damping are much weaker than the noise, they can
influence the distribution only insignificantly, so that
0 < ρλðAÞ ð0Þ < ∞ still holds. Let Q be an orthogonal
matrix diagonalizing A. The kth component of QT Δ ¼−QT A−1 QQT ∇f̂ðxÞ is proportional to 1=λðAÞk , and thus,
according to (7), its distribution has the heavy tail 1=z2 . We
can thus conclude [76,77] that the distribution of jjΔjj also
has the heavy tail 1=z2 . In the continuous time limit, thisleads to a superdiffusive, multidimensional LW or LF[76,78–82]. Note that the components of Δ are notindependent and the spectral measure [76,81] takes anontrivial form, which will be the subject of future studies.
In contrast, if the shift of the eigenvalue distribution related
to the curvature and damping factor is strong enough, thetails can be cut off completely, due to the bounded supportof the Wigner semicircle distribution. In this case, thecontinuous time limit process corresponds to diffusive
search. Although for any D < ∞ the cutoff formallydisappears, this shows that the heavy tails can be strongly
tempered by the damping factor and the curvature.
Importantly, the heavy tails may be tempered in the
directions of large curvatures while being preserved in
the other directions, thus providing a flexible combinationof random and directed search mechanisms.
We finally test our results using computer simulations
in a simple example of a search in an unbounded

250601-3

PRL 119, 250601 (2017)

PHYSICAL REVIEW LETTERS

FIG. 1. An example of search processes in an unbounded twodimensional space. The optimized function f is a sum of 20
Gaussians uniformly distributed within a 100 × 100 cell, which isperiodically repeated across an infinite space. The targets aresparse, so that in most places f is flat and does not provide any
information about the position of the targets due to measurement
noise. The top two plots represent a cell of f (yellow corresponds
to high values, blue to low values) and 20 exemplary trajectories
of searches starting from randomly chosen positions within thecell (blue dots) and finishing at a target (green dots) or at somerandom position without finding the target (red dots) due to thetime limitation (104 steps). The bottom left plot depicts the
distributions of jump lengths averaged over time and an ensembleof 104 trajectories. The second-order search produces a powerlaw tail with an exponent α ≈ 1. The bottom right plot shows thescaling of a displacement with time. As expected, for short times
the first-order search leads to a diffusive behavior which scales as
n1=2 , whereas the second-order search with strong noise σ H leads
to a superdiffusive behavior which scales as n. For longer times,
the median displacement saturates due to trapping at the targets.

two-dimensional space; see Fig. 1. First, the first-ordersearch in (3) does not produce power-law tails, and its
searching trajectory corresponds to that of normal diffusion
except near a target.
pﬃﬃﬃThe median distance from the initial
position scales as n with time step n. Next, as predictedby our theory, the second-order search with a sufficiently
large Hessian noise (σ H ≈ f 00 þ β) produces heavy tails of
the jump length distribution. In this case, the median
distance from the initial position scales linearly with time,
which is a characteristic of the LF process with α ¼ 1. If the
Hessian noise σ H is a few times smaller than f 00 þ β, thebehavior of the second-order method is in the middle
between the above two extreme cases. In this case, theheavy tail of the step size distribution is present but
somewhat tempered, in line with (11). As a consequence,
the growth ofpthe
ﬃﬃﬃ median distance from the initial position
is initially ∼ n, similarly to normal diffusion, but slowly

week ending
22 DECEMBER 2017

builds up in time as ∼n due to infrequent LF-like jumps.
Note that, in the current simulation setup, the median
distance saturates before the generalized central limit
theorem predicts its linear growth, because the trajectories
are trapped by the targets sooner. This indicates that thesecond-order method with an appropriate noise level (or thedamping constant) can find a target faster than the firstorder search or the ad hoc combination of directed searchwith power-law step sizes. The detailed analysis of the
optimal choice of the damping factor will be given elsewhere, but, intuitively, it is beneficial to perform the firstorder search in the direction of a convex surface andperform the LF search along the direction of a concave orplateau surface.
In summary, we analyze a stochastic version of Newton’s
optimization method. We argue that noise in the estimatesof the Hessian leads to a heavy-tailed distribution of jumps,
an indicator of LWs or LFs. We present conditions in
D ¼ 1 and D → ∞ under which the appearance of theheavy tails is guaranteed, and we corroborate these findings
with computer simulations in the biologically relevant caseof D ¼ 2.
Our model explains how a seemingly common behavior(LW with the fixed tail index α ¼ 1) can emerge from ageneric and locally optimal search strategy in the presenceof noise. This proposal is consistent with the role ofevolution and adaptation under selection pressures inacquiring an advantageous search strategy. However, unlike
some earlier proposals [22,35,83] that explore over theentire range of plastic α, our model gives only two possible
rigid values for α: α ¼ 1 generally and α ¼ 2 in the limitof ðf 00 þ βÞ=σ H → ∞. Moreover, our approach is distinct
from the Lévy flight foraging hypothesis [22], because thepossible evolutionary optimization is not carried within thefamily of random search strategies (e.g., over a rangeof α) but explicitly involves directed searches. Indeed,
multiple specific mechanisms shaping directed search havebeen previously shown to produce LWs with an α ¼ 1, e.g.,
movements in narrow, confined environments [45], in bulkmediated effective surface diffusion [45,84], and in patchy
environments if foragers use information about patch
quality [49]. Our proposal has two advantages over suchfindings: First, the second-order gradient-based optimization method is a well-established generic search strategy
that works efficiently in many different environmental
conditions. Second, by including taxis, our model suggests
a specific continuous crossover between the random anddirected search strategies as we further discuss below.
Our results suggest that some organisms may performtaxis according to Newton’s (or some other second-order)
optimization method, which should be possible to verify
experimentally. The resulting random walks are inhomogeneous and anisotropic, with less jerky motion along
directions with larger curvatures or weaker measurement
noise. These are distinct features of our model that can be

250601-4

PRL 119, 250601 (2017)

PHYSICAL REVIEW LETTERS

taken advantage of in experiments aiming to assess whether
foraging organisms employ second-order derivatives.
The candidate organisms that use taxis and, in someconditions, perform LWs with α ≈ 1 include microglia[7], Dictyostelium discoideum [8,9], and Drosophila[11]. The strategy we introduce combines the characteristics of two algorithms that are known to be efficient in
directed (second-order optimization) and random (LFs or
LWs) search scenarios. This method should therefore
perform well in a broad range of scenarios of stochastic
optimization, which may be of interest for the machinelearning community.
This work was supported by Brain/MINDS from Japan
Agency for Medical Research and Development and
RIKEN Brain Science Institute.

*nalewkoz@gmail.com
taro.toyoizumi@brain.riken.jp[1] M. F. Shlesinger, J. Klafter, and Y. M. Wong, J. Stat. Phys.
27, 499 (1982).
[2] M. F. Shlesinger and J. Klafter, in On Growth and Form
(Springer, New York, 1986), pp. 279–283.
[3] M. F. Shlesinger, B. J. West, and J. Klafter, Phys. Rev. Lett.
58, 1100 (1987).
[4] B. B. Mandelbrot, The Fractal Geometry of Nature(Freeman, New York, 1983).
[5] R. Metzler and J. Klafter, Phys. Rep. 339, 1 (2000).
[6] T. H. Harris et al., Nature (London) 486, 545 (2012).
[7] Y. Y. Grinberg, J. G. Milton, and R. P. Kraig, PLoS One 6,e19294 (2011).
[8] H. Takagi, M. J. Sato, T. Yanagida, and M. Ueda, PLoS One
3, e2648 (2008).
[9] A. M. Reynolds, Physica (Amsterdam) A389, 273 (2010).
[10] G. Ariel, A. Rabani, S. Benisty, J. D. Partridge, R. M.
Harshey, and A. Be’Er, Nat. Commun. 6, 8396 (2015).
[11] A. M. Reynolds and M. A. Frye, PLoS One 2, e354 (2007).
[12] A. M. Reynolds, A. D. Smith, R. Menzel, U. Greggers,
D. R. Reynolds, and J. R. Riley, Ecology 88, 1955 (2007).
[13] S. Wolf, E. Nicholls, A. M. Reynolds, P. Wells, K. S. Lim,R. J. Paxton, and J. L. Osborne, Sci. Rep. 6, 32612 (2016).
[14] G. M. Viswanathan, V. Afanasyev, S. V. Buldyrev, E. J.
Murphy, P. A. Prince, and H. E. Stanley, Nature (London)
381, 413 (1996).
[15] N. E. Humphries, H. Weimerskirch, and D. W. Sims,
Methods Ecol. Evol. 4, 930 (2013).
[16] D. W. Sims et al., Nature (London) 451, 1098 (2008).
[17] D. Brockmann, L. Hufnagel, and T. Geisel, Nature (London)
439, 462 (2006).
[18] M. C. Gonzalez, C. A. Hidalgo, and A.-L. Barabasi, Nature(London) 453, 779 (2008).
[19] I. Rhee, M. Shin, S. Hong, K. Lee, S. J. Kim, and S. Chong,IEEE/ACM Trans. Netw. 19, 630 (2011).
[20] D. Brockmann and T. Geisel, Neurocomputing;Variable Star
Bulletin 32–33, 643 (2000).
[21] M. E. Costa, F. Bonomo, and M. Sigman, Front. Integr.
Neurosci. 3 (2009).
†

week ending
22 DECEMBER 2017[22] G. M. Viswanathan, S. V. Buldyrev, S. Havlin, M. G. E.
Da Luz, E. P. Raposo, and H. E. Stanley, Nature (London)
401, 911 (1999).
[23] S. V. Buldyrev, S. Havlin, A. Y. Kazakov, M. G. E. da Luz,
E. P. Raposo, H. E. Stanley, and G. M. Viswanathan, Phys.
Rev. E 64, 041108 (2001).
[24] F. Bartumeus, J. Catalan, U. L. Fulco, M. L. Lyra, and G. M.
Viswanathan, Phys. Rev. Lett. 88, 097901 (2002).
[25] E. P. Raposo, S. V. Buldyrev, M. G. E. da Luz, M. C. Santos,
H. E. Stanley, and G. M. Viswanathan, Phys. Rev. Lett. 91,
240601 (2003).
[26] F. Bartumeus, M. G. E. da Luz, G. M. Viswanathan, and J.
Catalan, Ecology 86, 3078 (2005).
[27] Z. Pasternak, F. Bartumeus, and F. W. Grasso, J. Phys. A 42,
434010 (2009).
[28] E. P. Raposo, F. Bartumeus, M. G. E. Da Luz, P. J. RibeiroNeto, T. A. Souza, and G. M. Viswanathan, PLoS Comput.
Biol. 7, e1002233 (2011).
[29] V. V. Palyulin, A. V. Chechkin, and R. Metzler, Proc. Natl.
Acad. Sci. U.S.A. 111, 2931 (2014).
[30] N. E. Humphries and D. W. Sims, J. Theor. Biol. 358, 179(2014).
[31] L. Kusmierz, S. N. Majumdar, S. Sabhapandit, and G.
Schehr, Phys. Rev. Lett. 113, 220602 (2014).
[32] Ł. Kuśmierz and E. Gudowska-Nowak, Phys. Rev. E 92,
052127 (2015).
[33] A. M. Reynolds, J. Phys. A 42, 434006 (2009).
[34] F. Bartumeus, Fractals 15, 151 (2007).
[35] G. M. Viswanathan, E. P. Raposo, and M. G. E. Da Luz,
Phys. Life Rev. 5, 133 (2008).
[36] D. D. Taub, K. Conlon, A. R. Lloyd, J. J. Oppenheim, and
D. J. Kelvin, Science 260, 355 (1993).
[37] M. D. Gunn, K. Tangemann, C. Tam, J. G. Cyster, S. D.
Rosen, and L. T. Williams, Proc. Natl. Acad. Sci. U.S.A. 95,
258 (1998).
[38] J. Adler, Science 153, 708 (1966).
[39] H. C. Berg and D. A. Brown, in Chemotaxis: Its Biology andBiochemistry (Karger, Basel, 1974), pp. 55–78.
[40] H. Szurmant and G. W. Ordal, Microbiol. Mol. Biol. Rev.
68, 301 (2004).
[41] Y. Zhu, A. Nern, S. L. Zipursky, and M. A. Frye, Curr. Biol.
19, 613 (2009).
[42] M. S. Miller and T. S. Keller, J. Gravit. Physiol. 6, P99(1999).
[43] B. J. Duistermars, D. M. Chow, and M. A. Frye, Curr. Biol.
19, 1301 (2009).
[44] S. Wasserman, A. Salomon, and M. A. Frye, Curr. Biol. 23,
301 (2013).
[45] A. M. Reynolds, Phys. Life Rev. 14, 59 (2015).
[46] D. Boyer, O. Miramontes, G. Ramos-Fernandez, J.
Mateos, and G. Cocho, Physica (Amsterdam) A 342,
329 (2004).
[47] M. C. Santos, D. Boyer, O. Miramontes, G. M. Viswanathan,E. P. Raposo, J. L. Mateos, and M. G. E. da Luz, Phys. Rev. E
75, 061114 (2007).
[48] A. M. Reynolds, Phys. Rev. E 78, 011906 (2008).
[49] A. M. Reynolds, J. R. Soc. Interface 9, 1568 (2012).
[50] A. M. Reynolds, Proc. R. Soc. A 470, 20140408 (2014).
[51] M. F. Shlesinger, J. Chem. Phys. 78, 416 (1983).
[52] A. M. Reynolds, Sci. Rep. 4, 4409 (2014).

250601-5

PRL 119, 250601 (2017)

PHYSICAL REVIEW LETTERS[53] D. W. Sims, A. M. Reynolds, N. E. Humphries, E. J.
Southall, V. J. Wearmouth, B. Metcalfe, and R. J. Twitchett,Proc. Natl. Acad. Sci. U.S.A. 111, 11073 (2014).
[54] A. Ott, J. P. Bouchaud, D. Langevin, and W. Urbach, Phys.
Rev. Lett. 65, 2201 (1990).
[55] N. Srivastava, D. A. Clark, and A. D. T. Samuel, J. Neurophysiol. 102, 1172 (2009).
[56] L. C. M. Salvador, F. Bartumeus, S. A. Levin, and W. S.
Ryu, J. R. Soc. Interface 11, 20131092 (2014).
[57] T. S. Biró and A. Jakovác, Phys. Rev. Lett. 94, 132302(2005).
[58] I. Lubashevsky, R. Friedrich, and A. Heuer, Phys. Rev. E 79,
011110 (2009).
[59] A. M. Reynolds, P. Schultheiss, and K. Cheng, Behav. Ecol.
Sociobiol. 67, 1219 (2013).
[60] E. Korobkova, T. Emonet, J. M. G. Vilar, T. S. Shimizu, and
P. Cluzel, Nature (London) 428, 574 (2004).
[61] Y. Tu and G. Grinstein, Phys. Rev. Lett. 94, 208101(2005).
[62] A. L. Barabási, Nature (London) 435, 207 (2005).
[63] F. Matthäus, M. Jagodič, and J. Dobnikar, Biophys. J. 97,
946 (2009).
[64] F. Matthäus, M. S. Mommer, T. Curk, and J. Dobnikar,
PLoS One 6, e18623 (2011).
[65] A. M. Reynolds, Physica A (Amsterdam) 390, 245 (2011).
[66] G. G. Katul, A. Porporato, R. Nathan, M. Siqueira, M. B.
Soons, D. Poggi, H. S. Horn, and S. A. Levin, Am. Nat. 166,368 (2005).
[67] M. W. Shaw, T. D. Harwood, M. J. Wilkinson, and L. Elliott,
Proc. R. Soc. B 273, 1705 (2006).
[68] A. M. Reynolds, Am. Nat. 181, 555 (2013).
[69] S. Petrovskii, A. Mashanova, and V. A. Jansen, Proc. Natl.
Acad. Sci. U.S.A. 108, 8704 (2011).

week ending
22 DECEMBER 2017[70] S. G. Nurzaman, Y. Matsumoto, Y. Nakamura, S. Koizumi,
and H. Ishiguro, in Proceedings of the IEEE International
Conference on Robotics and Biomimetics, 2008, ROBIO
2008 (IEEE, New York, 2009), pp. 806–811.
[71] V. Zaburdaev, S. Denisov, and J. Klafter, Rev. Mod. Phys.
87, 483 (2015).
[72] E. P. Wigner, Ann. Math. 67, 325 (1958).
[73] D. Voiculescu, Inventiones Mathematicae 104, 201 (1991).
[74] J. H. Schenker and H. Schulz-Baldes, arXiv:math-ph/0505003.
[75] L. Erdős, Russ. Math. Surv. 66, 507 (2011).
[76] G. Samorodnitsky and M. S. Taqqu, Stable Non-Gaussian
Random Processes: Stochastic Models with Infinite Variance(CRC Press, Boca Raton, FL, 1994), Vol. 1.
[77] J. P. Nolan, Modeling financial data with stable distributions, in Handbook of Heavy Tailed Distributions in
Finance, edited by S. T. Rachev (Elsevier-North Holland,
Amsterdam, 2005).
[78] M. Teuerle and A. Jurlewicz, Acta Phys. Pol. B 40, 1333
(2009).
[79] M. Teuerle, P. Żebrowski, and M. Magdziarz, J. Phys. A 45,
385002 (2012).
[80] K. Szczepaniec and B. Dybiec, Phys. Rev. E 90, 032128(2014).
[81] K. Szczepaniec and B. Dybiec, J. Stat. Mech. (2015)
P06031.
[82] B. Dybiec and K. Szczepaniec, Eur. Phys. J. B 88, 184(2015).
[83] M. E. Wosniack, M. C. Santos, E. P. Raposo, G. M.
Viswanathan, and M. G. da Luz, PLoS Comput. Biol. 13,
e1005774 (2017).
[84] O. V. Bychuk and B. OShaughnessy, J. Chem. Phys. 101,
772 (1994).

250601-6