Stochastic Optimal Foraging: Tuning Intensive and Extensive Dynamics in Random Searches

Recent theoretical developments had laid down the proper mathematical means to understand how the structural complexity of search patterns may improve foraging efficiency. Under information-deprived scenarios and specific landscape configurations, Lévy walks and flights are known to lead to high search efficiencies. Based on a one-dimensional comparative analysis we show a mechanism by which, at random, a searcher can optimize the encounter with close and distant targets. The mechanism consists of combining an optimal diffusivity (optimally enhanced diffusion) with a minimal diffusion constant. In such a way the search dynamics adequately balances the tension between finding close and distant targets, while, at the same time, shifts the optimal balance towards relatively larger close-to-distant target encounter ratios. We find that introducing a multiscale set of reorientations ensures both a thorough local space exploration without oversampling and a fast spreading dynamics at the large scale. Lévy reorientation patterns account for these properties but other reorientation strategies providing similar statistical signatures can mimic or achieve comparable efficiencies. Hence, the present work unveils general mechanisms underlying efficient random search, beyond the Lévy model. Our results suggest that animals could tune key statistical movement properties (e.g. enhanced diffusivity, minimal diffusion constant) to cope with the very general problem of balancing out intensive and extensive random searching. We believe that theoretical developments to mechanistically understand stochastic search strategies, such as the one here proposed, are crucial to develop an empirically verifiable and comprehensive animal foraging theory.

In the following, we calculate every quantity necessary to determine the search efficiency η and the average distance traversed between consecutive encounters L , generally expressed by Eqs. (4) and (6) in the main manuscript. We review the developments due to Buldyrev et al. [1,2], which have been also explored in Ref. [3]. The results will be given in a general form, as function of the distribution of step lengths p(ℓ). In Appendix S2 we apply them to some relevant cases of reorientation strategies, such as Lévy, stretched exponential, log-normal and gamma distributions.

Calculation of L and |ℓ|
In order to calculate the average distance L traversed by a searcher starting from the position x 0 = a until reaching one of the borders located at x = 0 and x = λ, one must deal with the the survival probability that the walker has not yet encountered any of the boundary targets after n moves: where ρ n (x n ) represents the pdf for the walker to lie between x n and x n + dx n after n moves. Actually, an interesting connection with the first-passage-time problem [4] can be established since L can be also interpreted as the average distance traversed until the searcher reaches for the first time any of the boundaries. The complementary probability of finding any of the targets in some move n ′ ≥ n+1 is thus P n ′ ≥n+1 = 1 − P not n , so that the probability of finding a target precisely after n moves reads P n = |P n ′ ≥n+1 − P n ′ ≥n | = |P not n − P not n−1 |.
This result leads to so that ρ n−1 (x) − ρ n (x) can be interpreted as the pdf to encounter a target precisely after n moves.
To calculate L , we must weight with P n the average distance traversed up to an encounter after n steps. By denoting |ℓ| (x) as the average length of a single step starting at position x, which is a quantity independent of n, some algebra [1,2] Notice that a recursion formula for ρ n (x), can be obtained from the general result leading to The above result can be put in a much simpler form by defining the integral operator [1,2] [ that is, It is interesting to note that [Lρ n ](x) ≤ |L max | λ−rv rv ρ n (x ′ )dx ′ , where |L max | means [L 1] calculated at the middle point x ′ = λ/2, which gives rise to its maximum value. Therefore, we find that |L max | = λ−rv rv p(λ/2 − x)dx < 1. This property allows the series expansion where I denotes the unitary operator: [Iρ](x) = ρ(x), so that, by inserting the initial condition ρ 0 ( This closed analytical expression is actually essential to determine the efficiency of the search, according to Eq. (4) in the main manuscript. In Eq. (11), the average single step length is given by in the case r v + ℓ 0 ≤ a ≤ λ − r v − ℓ 0 (i.e. for the largest part of the search space). The first two integrals represent moves to the left and to the right which are not truncated by the encounter of a target. The third and fourth represent moves truncated by the encounter of the targets, respectively, at x = 0 and x = λ. In fact, due to the perceptual range or radius of vision of the searcher, these sites are detected as soon as the walker reaches the respective positions x = r v and x = λ − r v . In addition, since p(ℓ) = 0 if |ℓ| < ℓ 0 , then in the complementary cases |ℓ| (a) is given only by the second, third and fourth (first, third and fourth) integrals for The exact formal expression (11) can be numerically evaluated through a spatial discretization [1,2] of the continuous range 0 ≤ x ≤ λ. In order to accomplish it, we consider positions x which are multiple of some discretization length ∆x, i.e. x = j∆x, with j = 0, 1, ..., M and ∆x much smaller than any relevant scale of the problem (ℓ 0 , r v , λ). In this case, the targets at x = 0 and x = λ are respectively associated with the indexes j = 0 and j = M = λ/∆x (M is the integer number of intervals of length ∆x in which the range 0 ≤ x ≤ λ is subdivided). Similarly, we define ℓ 0 = m 0 ∆x and r v = m r ∆x, with m 0 and m r integers. The continuous limit is recovered by taking ∆x → 0 and M → ∞, with λ = M ∆x fixed.
With these considerations in mind, Eq. (5) can be discretized to in which the structure of a sequence of matrix products appears by regarding the quantities a k,j as the matrix elements [A] k,j of a symmetric matrix A, with null diagonal elements and dimension (M − 2m r − 1) × (M − 2m r − 1). As a consequence, the discrete equivalent of Eq. (11) becomes where I is the (M − 2m r − 1) × (M − 2m r − 1) unity matrix and (I − A) −1 is the inverse of the matrix (I − A). Finally, by noticing that [A] k,j is the discrete representation of the probability p(x − x ′ )dx ′ of performing a move of length between |x − x ′ | = |k − j|∆x and |x − x ′ | + dx ′ = (|k − j| + 1)∆x, we obtain [A] j,j = 0 and [A] k,j = 0 if |k − j| < m 0 due to the minimum step length ℓ 0 .
Calculation of p 0 and p λ Two other relevant quantities are the probabilities for a searcher that starts at x 0 = a to find, respectively, the target at x = 0 or x = λ, denoted by p 0 (a) and p λ (a). We start by noticing [1,2] that, when the target at x = λ is found after n steps, then where ρ n−1 (x n−1 )dx n−1 represents the probability for the walker to be located in the interval [x n−1 , x n−1 + dx n−1 ) after n − 1 moves, and we multiply the probability that the next (n-th) move will reach the target at x = λ and terminate the walk. When summing over all n, a procedure similar to that described in the last subsection takes place, so that we obtain from which we also readily determine p 0 (a) = 1 − p λ (a). The numerical implementation of Eq. (17) also follows the above discussed procedure.

Calculation of L 0 and L λ
We now obtain the average distances L 0 and L λ traversed by a searcher starting at x 0 = a to find the target at x = 0 or x = λ, respectively. We first notice that Eq. (12) can be rewritten as since p(ℓ) = 0 for |ℓ| < ℓ 0 . Our aim is to separate in Eq. (18) the search walks that terminate at x = 0 from those that terminate at x = λ. This is already achieved in the second and third integrals above, which represent truncation by target encounter at x = 0 and x = λ, respectively. However, the first integral in Eq. (18) actually represents moves which are not truncated, i.e. the search walk does not end at these moves. In fact, in this case the process stops only at some subsequent step. By introducing the identity p λ (x) + p 0 (x) = 1 in the the first integral in Eq. (18) and substituting the result into Eq. (11), we obtain Above, we actually separate the events associated with the finding of the target at x = λ and at x = 0. By comparison with Eq. (6) in the main manuscript, we identify [1,2] L λ (a) = 1 and with p λ (x) and p 0 (x) = 1 − p λ (x) given by Eq. (17). We shall now write, say L λ , in the discrete form, proper for the numerics. From Eq. (20) we generally have where the column vector F ι0 is given by the discrete limit of the continuous function We proceed by separating F (a) into two parts [1,2]: The simplest one is As previously explained, after integration the setting of discrete parameters should be considered. We next focus on By transforming this integral into a discrete sum, taking the discrete limit of |x − a|, which is |i − ι 0 |∆x, and substituting the probability p(x − a)dx of giving a move of length |x − a| by the corresponding matrix element A i,ι0 (= A ι0,i ), the (M − 2m r − 1) × 1 column vector F 1,ι0 is obtained: From Eq. (24) the sum of F 1,ι0 and F 2,ι0 results in the column vector F ι0 , which is inserted into Eq. (22) to give L λ ι0 . A similar approach leads to the discrete version of L 0 . However, we can always make use of the discrete form of Eq. (6) in the main manuscript, so that Appendix S2: Application to specific reorientation strategies We apply the general calculations of Appendix S1 to some relevant cases, where the reorientation strategy (characterized by p(ℓ)) is given by a power-law Lévy (discussed in [1][2][3]), and the new cases of upper truncated Lévy, stretched exponential, log-normal and gamma distributions. With their explicit calculations a comparison between efficiencies, as well as the study of the role of nearby/distant visits and the search dynamics in each case, become possible.

Truncated Lévy power-law distributions
We start by considering the following pdf of move lengths where the theta function is such that, e.g. Θ(|ℓ| − ℓ 0 ) = 0 if |ℓ| < ℓ 0 and Θ(|ℓ| − ℓ 0 ) = 1. The presence of Θ(|ℓ| − ℓ 0 ) above assures the minimum move length ℓ = ℓ 0 . The introduction of the other theta function Θ(|ℓ| − τ ) imposes a maximum move length ℓ = τ , as justified below. The normalization constant A is calculated from Eq. (1) in the main manuscript, to give Clearly, in the limit τ → ∞ the pdf (29) presents no upper cutoff, with possible move lengths extending up to infinity. In such a case, the power-law dependence of Eq. (29) represents the long-range asymptotic limit of the complete family of α-stable Lévy distributions of index α = µ − 1 [5,6]. As the second moment of pdf (29) with τ → ∞ diverges for 1 < µ ≤ 3, the central limit theorem does not hold, and anomalous superduffisive dynamics takes place, governed by the generalized central limit theorem.
Indeed, Lévy walks and flights in free space are related to a Hurst exponent [5,6] H > 1/2, whereas Brownian behavior (diffusive walks with H = 1/2) emerges for µ > 3. In particular, for Lévy walks one finds H = 1 for 1 < µ ≤ 2, with ballistic dynamics emerging as µ → 1 (the case µ = 2 corresponds to the Cauchy distribution). Therefore, by varying the single parameter µ in Eq. (29), the whole range of search dynamics can be accessed (from Brownian to superdiffusive and ballistic).
We also remark that, as one considers the search path as a whole, the truncation of steps by the finding of a statistically large number of targets generates an effective truncated Lévy distribution [7], with finite moments and emergence of a crossover towards normal dynamics in the long term, as justified by the central limit theorem.
Below we illustrate the calculations for some intervals and µ = 2. The extension to the other cases is straightforward. By substituting Eqs. (29) and (30) into Eq. (12), the results read as follows: i. for ii. for a > r v + τ and a > λ − r v − τ : iii. for a ≤ r v + τ and a ≤ λ − r v − τ : Discrete space expressions associated with the above equations for |ℓ| (a) are readily found by following the prescription described for L and |ℓ| in Appendix S1. For example, from Eq. (30) and the definition τ = m τ ∆x the discrete limit of Eq. (32) reads The matrix A can be also determined by substituting Eqs. (29) and (30) into Eq. (15), so that with A ij = 0 if i = j, |i − j| < m 0 or |i − j| > m τ . At this point we recall that the substitution of |ℓ| ι0 into Eq. (14), along with Eq. (36), leads to L ι0 and, therefore, also to the efficiency η, Eq. (4) in the main manuscript, in the case of truncated Lévy searches. In order to compute the quantity p 0 (and p λ = 1 − p 0 ) we first need to calculate the probability of moves of length ℓ ≥ λ − r v − a, P (ℓ ≥ λ − r v − a). Specifically, for the case of truncated Lévy walks we obtain The discrete version of these expressions can be taken straightforwardly.

Stretched exponential distributions
For stretched-exponential search walks with minimum move length ℓ 0 , p(ℓ) reads with φ > 0, 0 < β ≤ 1, and, from Eq. (1) in the main manuscript, Note also that the limit value β = 1 corresponds to the exponential distribution. The exponential factor in Eq. (38) introduces an effective attenuation of the long-tailed power-law decay. Indeed, we observe that the β → 0 + and φ → 0 + cases tend to the τ → ∞ Lévy distributions (29) with µ → 1 + and µ = 1 − β, respectively. By substituting Eqs. (38) and (39) into Eq. (12) we find for the interval whereφ = φℓ β 0 and γ(k, x) represents the lower incomplete gamma function, The discretization of Eq. (40) is taken analogously to the Lévy case. The matrix A for the stretched exponential distribution can be determined by substituting Eqs. (38) and (39) into Eq. (15), to obtain with A ii = 0. In order to compute p 0 and p λ , we find for the stretched exponential pdf and P (ℓ ≥ λ − r v − a) = 1/2 otherwise.