key: cord-0938694-fb10pxon
authors: Vilk, Ohad; Campos, Daniel; M'endez, Vicencc; Lourie, Emmanuel; Nathan, Ran; Assaf, Michael
title: Phase transition in non-Markovian animal exploration model with preferential returns
date: 2021-11-30
journal: Physical review letters
DOI: 10.1103/physrevlett.128.148301
sha: a7376d7ef5011af68a51e6b7549a0369ed0197f7
doc_id: 938694
cord_uid: fb10pxon

We study a non-Markovian and nonstationary model of animal mobility incorporating both exploration and memory in the form of preferential returns. We derive exact results for the probability of visiting a given number of sites and develop a practical WKB approximation to treat the nonstationary problem. We further show that this model adequately describes empirical movement data of Egyptian fruit bats (Rousettus aegyptiacus) when accounting for inter-individual variation in the population. Finally, we study the probability of visiting any site a given number of times and derive the corresponding mean-field equation. Here, we find a remarkable phase transition occurring at preferential returns which scale linearly with past visits. Following empirical evidence, we suggest that this phase transition reflects a trade-off between extensive and intensive foraging modes.

Introduction. Movement is a vital part of life and is key in a wide range of physical, biological and ecological systems. Theoretical and empirical frameworks are thus amply used to study the mechanisms underlying movement patterns in all organisms [1] . In particular, individualbased modeling of movement has played a crucial role in studying dynamic systems across multiple spatiotemporal scales [2] [3] [4] . These models can be applied to infer behaviors and draw causal links between observed phenomena and their underlying mechanisms beyond phenomenological description of the observed patterns [5] .

Most theoretical models rely on Markovian assumptions to capture the properties of animal trajectories. However, memory and similar cognitive mechanisms are key to understanding patterns observed in animal foraging [6, 7] . A range of taxa, from insects to primates, have been shown to exhibit spatial memory, and many of them are known to repeatedly return to previously visited sites as a part of their regular foraging strategies, mammals being a paradigmatic example [4, [8] [9] [10] [11] . Notably, memory patterns must be properly balanced by the organisms with some level of behavioral plasticity to enhance flexibility and exploration (see, e.g., [12] ). For all these reasons, correctly incorporating memory within stochastic models is an important research line for improving both predictive and descriptive tools of movement [6, 13-15]. Indeed, it has been shown that heuristic models of memory can be derived from microscopic consideration for limited cases [16, 17] . At the same time, new experimental methods are allowing to disentangle, even under field conditions, memory effects on movement from other cognitive/perception mechanisms [9] . As long as such experimental and theoretical advances are able to nourish * To whom correspondence should be addressed: michael.assaf@mail.huji.ac.il each other, new levels of detail in our understanding of living organisms can be potentially reached.

Dealing with memory and similar cognitive mechanisms still represents a significant theoretical challenge. Stochastic models that allow the individual to return to its original position (resettings) have attracted much attention recently [18] [19] [20] , but these only incorporate memory in an elementary way. More complex self-avoiding random walks or preferential returns (PR), where the individual returns to any previous location with a probability proportional to the number of previous visits have also been studied [21, 22] . These models are non-Markovian (and typically also non-stationary), requiring that the individual identifies and keeps record of its entire trajectory. While the propagator of these models [23] , and some properties of relocation times [24] , have been computed, characterizing the revisits complete statistics to each particular location remains an open problem.

Here we study a non-Markovian and non-stationary mechanistic model of animal mobility, explicitly incorporating both the tendency of an individual to return to previously visited locations (PR) and the tendency to explore new sites. Versions of this model have been used to model the mobility of humans [21] and monkeys [25] , the latter suggesting that monkey movements are non-random due to the use of memory and visitation patterns driven by resource availability. We generalize the model by accounting for stochasticity, incorporating inter-specific variability in the population, and allowing for nonlinear PR [26] . We provide analytical solutions to this non-Markovian, non-stationary model that go well beyond previous mean-field results. In particular, we present several approaches to analytically find the (nonstationary) probability of having visited n sites at time t and study the statistics of how revisits are distributed through the available locations. Our approach, based on the WKB (Wentzel-Kramers-Brillouin) approximation, is thus useful to deal with explicitly time-dependent and arXiv:2111.15231v1 [cond-mat.stat-mech] 30 Nov 2021 non-stationary problems. Remarkably, by allowing for nonlinear PR we find a phase transition as a function of the strength of the PR, where above some threshold the most visited site dominates the dynamics, receiving practically all new visits. We suggest that this phase transition reflects a balance between the tendency to return to known sites and the will to explore new ones [19] . We further verify our predictions using simulations, and show that our theoretical results adequately describe the space use patterns and the revisitation dynamics of Egyptian fruit bats (Rousettus aegyptiacus) to fruit trees.

Model. Our model consists of two elements [21]: (i) exploration -with probability P new the animal visits a new site, and (ii) PR -with probability 1 − P new the animal visits a previously visited site i with probability Π i (m i ), where m i is the number of previous visits to site i. Following empirical evidence in humans and animals [21, 25] we assume that

(1)

Here n is the number of previously visited sites, and β > 0 and 0 < q < 1 control the animal's tendency to visit new sites indicating a power-law decay controlled by conformity exponent β. On the other hand, the PR exponent α > 0, governs the tendency to return to a previously visited location. Furthermore, and without loss of generality, we order the sites by rank such that i = 1 is the most visited site with m 1 visits. Notably, we assume that the number of available sites is always larger than the number of visited sites, and that no significant resource depletion within a site occurs.

Cumulative number of sites. The probability P (n, t) of having visited n sites in t 1 time steps follows ∂P (n, t)/∂t = q(n − 1) −β P (n − 1, t) − qn −β P (n, t).

(2)

Although this master equation is interpreted here in the context of movement between spatially distributed sites, it can equivalently describe a birth-death process of population of size n, where the growth rate is proportional to n −β [27]. In particular, for β = 0 the birth-death process is ∅ q − → A and for β = −1 the birth-death process is A q − → 2A. While these special cases have known exact solutions, in this manuscript we are primarily interested in the regime β > 0, which describes a growth which decreases [or saturates, see Eq. (1)] with the number of sites (or with the population size). To the best of our knowledge this regime has not been analytically studied.

An equation for the first moment n = n nP (n, t) can be obtained from Eq. (2) by multiplying the latter by n, summing over all n's, and using the definition for n , resulting in d n /dt = q n −β which under the meanfield approximation n −β n −β is solved by [21]

predicting a power-law dependence on the time of measurement. A similar derivation for the second moment yields n 2 = n 2 + n such that the variance follows σ 2 n ≡ n 2 − n 2 = n , i.e., the variance of the number of sites is equal to the mean as in a Poisson process. This result, however, turns out to be inaccurate as it involves various uncontrolled assumptions, and is not consistent with simulations, as elaborated below.

An exact solution to Eq.

(2) can be found by Laplace transforming the equation and solving the resulting recurrence equation [28] . The exact solution for β = 0 has the form [see Supplementary Information (SI) Sec. S1.A] P (n, t) = (−1) n−1 n β n k=1 k −β e −qt/k β n j=1,j =k

(4)

For special values of β this result simplifies to

(5) Although Eq. (4) is an exact solution, it is given in form of a summation of large terms of alternating sign, which converges due to a precise balance between the terms. Thus, in practice this result is very slow to converge for n 1 in any finite-size system (e.g., python, matlab, and mathematica) and may result in a significant lack of accuracy. Moreover, virtually any approximation made in calculating the individual terms may cause large errors for n 1 [29]. To circumvent these issues, we develop a time-dependent WKB approximation.

Time dependent WKB. In the limit of a large number of sites n 1 [30-32], we substitute the timedependent ansatz P (n, t) ∼ e −S(n,t) into Eq. (2). Neglecting terms of order O(n −1 ) we obtain a classical Hamilton-Jacobi equation for the action function S(n, t): ∂ t S = H(n, ∂ n S) ≡ H(n, p) where we have defined the Hamiltonian H(n, p) = q (1 − e −p ) n −β , and denoted p = −∂ n S as the conjugate momentum. Instead of directly solving the Hamilton-Jacobi equations, we use the Hamilton approach for the classical equations of motioṅ

We write the action on a classical trajectory as [31]:

where the energy E ≡ H[n(t), p(t)] is constant along a dynamical trajectory given by p(n) = log q/(q − En β ) .

To find the energy we solve the equation of motion (6) on a given dynamical trajectory on which the energy is constant. After some algebra (SI, Sec. S1.C) this yields

is the incomplete beta function. This calculation of the probability of having visited n 1 sites at time t, is one of our main results. In the low energy limit, E 1, S(x) becomes (SI)

which can be shown to solve the Hamilton-Jacobi equation in the limit |x − 1| 1. In Fig. 1a we find good agreement between the exact result for the PDF [Eq. (4)], time-dependent WKB approximation [Eqs. (8, 9) ], and simulations (see also Fig. S1 in the SI). Here the exact result and WKB approximation are practically indistinguishable, whereas the low energy approximation predicts the PDF well only in its Gaussian vicinity. Notably, for n 1 the accuracy of the exact result rapidly deteriorates due to summation of (alternating) very large terms and accumulation of errors, making the time-dependent WKB approach highly advantageous in this case [29] . Equation (8) predicts a different variance than that predicted by the mean-field approach above. The variance can be found by approximating S(x) in the Gaussian vicinity of n = n . Doing so, we find S(x) (β + 1/2) (x − 1) 2 , which yields a variance of σ 2 n = n /|S (x)| x=1 = n /(1 + 2β). This entails that the distribution is narrower by a factor of (1 + 2β) than that predicted by the moment equations. In the inset of To account for between-individual variation, we generalize our model by allowing different β values acoss individuals. Assuming β is sampled from a normal distribution N (β 0 , σ 2 ) with mean β 0 and variance σ 2 1 (indicating the inter-individual variability in β around β 0 , see empirical results analysis below), P n (t) satisfies where P β (n, t) is the probability for a given β, see e.g., Eq. (8). Although analytical progress is possible in the limit of σ 1 (SI, Sec. S1.D), Eq. (10) can generally be solved numerically (Fig. 1b) . Notably, we checked that while variability in a population (σ > 0) will not significantly affect the mean number of sites, it dramatically affects the variability across a population (SI, Fig. S2 ).

Statistics of number of visits to a site. Having computed the statistics of number of sites, we now turn to study the probability W i (m i , t) of having m i visits at time t to an already visited site i, which follows

where P new and Π i are given by Eq.

(1). In general, Π i depends on the number of visits to other sites, such that Eq. (11) couples between different sites. Below, we focus on the limit t 1, where P new → 0 (SI, Sec. S2). The case of α = 1. In the limit t 1 we have i m i = t, namely, the total number of visits to all sites equals the total number of time steps t. Thus, Π i (m i ) = m i /t and P new = 0 [33] . Equation (11) is then

The case of α = 1. Here, we a priori assume that n j=1 m j α ∼ t ξ , where ξ is a priori unknown and satisfies α < ξ < 1 (see SI, Sec. S2.A for proof). The average number of visits to any site i can then be shown to asymptotically follow

Plugging this back into the sum over m j α we find that ξ = ξ 0 (1 + ), where ξ 0 = (1 + αβ)/(1 + β) and 1 is an unknown function of α, β (Fig. S3 ). For 1 − α [β = O(1)], we further find m i ∼ t β/(1+β) , which is independent of α; yet, the condition 1−α breaks down as α approaches 1 (SI, Sec. S2.B and Fig. S4 ). Importantly, in the limit t 1, for any α < 1 we find that all sites scale similarly with time. In contrast, for α > 1 not 

These results reveal a phase transition at α = 1 (see also SI, Sec. S2.D), where for weak PR (α < 1) the frequency of visits to the most visited site f 1 ≡ m 1 / n j=1 m j is only a small fraction of the total number of visits, while for strong PR (α > 1) f 1 approaches 1 as t is increased and site 1 dominates (Fig. 2a) . Importantly, in addition to the phase transition for f 1 , the next most visited sites (f k , for k = 2, 3, ...) peak around α = 1 (Figs. 2b and S5). Here, as α approaches 0 the number of visits to all sites becomes similar, while for α > 1 these visitation frequencies tend to zero.

Movement of fruit bats. To study the relevance of our model for real-life systems and to obtain insights into the phase transition, we compare our predictions to resource use patterns and visitation dynamics of wild fruit bats tracked by ATLAS during winter and summer [34] . In Fig. 3a -b we fit the mean and variance of the number of visited sites (fruit trees) as a function of the number of movement steps (defined here as distinguishable movement between trees, see SI) to our theory [35] . We find that during the summer β 0 and σ are higher than dur-ing the winter, entailing lower rate of visits to new sites (higher levels of conformity) and larger inter-individual variability, respectively. Notably, many empirical studies have shown that individual preferences and decisions can affect movement and behavior, such that individuals do not have identical movement patterns [36, 37] . This may explain the inter-individual variability observed in both summer and winter. In Fig. 3c-f we show that during summer m 1 ∼ t 0.97 and m 2 ∼ t 0.99 (as well as σ 2 m1 ∼ t 1.94 and σ 2 m2 ∼ t 2.16 ), which matches the theory for α = 1. In contrast, during winter m 1 ∼ t 0.89 and m 2 ∼ t 0.87 (σ 2 m1 ∼ t 1.96 , σ 2 m2 ∼ t 1.80 ), which is consistent with α values slightly below 1. These seasonal differences may be attributed to the fact that bats during the summer feed of highly abundant and palatable fruits -mulberries or figs with high levels of sugar contentand hence do not need to explore for feeding alternatives (high β 0 ) and can strongly rely on a limited number of sites (α = 1). Similarly, when food is abundant it makes sense that most bats will forage on the same bountiful locations, as they are not constrained by intense resource competition, thus reducing inter-individual variability in behaviour. In contrast, during winter there is less motivation to return to less favorable fruits (chinaberries) and a higher motivation to explore alternative trees, such as non-seasonal (unpredictable) fruit from the Ficus family.

In light of the phase transition at α = 1, and in agreement with experimental results, we hypothesize that in animal movement the value of α will tend towards 1. This maximizes the frequencies of visits to preferred sites (Fig. 2b ), yet avoids an exclusive choice of a preferred site which occurs at α > 1 (see Fig. 2a ). In this manner the animal combines intensive search patterns (committing to few site) with extensive searches (returning to all sites with some probability), a balance which is essential to optimize between energy expenditure and risk management [19, 38, 39] . Indeed, for fruit bats we find α ≈ 1, and similar results were obtained for humans [21] and monkeys [25] . Furthermore, the strategy of avoiding an exclusive site resembles bet-hedging strategies, e.g., bacterial persistence [40] . In addition, the value of α may be correlated to the total number of known sites: for large β 0 (few sites, summer) the bats will show stronger PR, while for smaller β 0 (more sites, winter) the strategy may tend towards a more uniform visitation rate to all trees.

More generally, we expect our results to also provide key insights into revisit dynamics in other areas like human mobility [41] [42] [43] , COVID-19 spread [44] , human migration [45] and languages dynamics [46] . [27] Similar growth rates appear in self-inhibitory gene regulatory networks, where a protein inhibits its own growth, see e.g., Refs. [47, 48] .

[28] In the SI, Sec. S1.B we alternatively develop a a generating function approach, and Eq.

(2) becomes a fractional integro-differential equation, which is solvable in specific cases. 

Here we provide additional details and results to support the derivations presented in the main text. Below, the notations and acronyms are the same as in the main text and the equations and figures refer to those therein.

In this section we detail the exact solution and the WKB approximation to Eq.

(2) of the main text dP (n, t) dt = q(n − 1) −β P (n − 1, t) − qn −β P (n, t).

(S1)

Here we derive Eqs. (4) and (5) of the main text by Laplace transforming Eq. (S1) and solving the resulting recurrence equation. First we define J(n, t) = qP (n, t)/n β , which turns Eq. (S1) into: In some special cases it is more convenient to solve Eq. (S1) using the generating function approach. Particularly, in the limit of a large number of sites, the equation in the generating function domain becomes a fractional integrodifferential equation for 0 < β ≤ 1, see below.

We define the generating function G(z, t) = n z n P (n, t), such that 2. Solution for β = −1 For β = −1 a similar derivation to Eq. (S11) yields a partial differential equation for the generating function: ∂ t G(z, t) = q(z − 1)z∂ z [G(z, t)], whose solution is G(z, t) = z/[z + e qt (1 − z)] [51] . Using Eq. (S10) we find P (n, t) = e −nqt e qt − 1 n−1 .

(S13)

The average number of sites here is n = e qt , i.e., we find exponential growth, as expected for a growth rate that is linear in n. Here, the variance is σ 2 n = e qt (e qt − 1) e 2qt = n 2 , which is significantly broader than that of the Poisson distribution. This result also agrees with Eq. (5) of the main text.

For β = 1, Eq. (S11) is rewritten in explicit integro-differential form:

Here, we Laplace transform Eq. (S14) in time

to obtain an integral equation. This equation can be solved iteratively by the Neumann series method [52] to give:

where γ(a, b) is the lower gamma function. Using Eq. (S10) we can inverse Laplace transform Eq. (S16) to obtain P (n, t) = 1 (n − 1)! n k=1 (−1) n−k k n−1 n k e − qt k , (S17) in agreement with Eq. (5) of the main text.

Here we derive Eqs. (8) and (9) of the main text. We employ the time-dependent WKB approximation in the limit of a large number of sites n 1 [31, 32]. Substituting the time-dependent ansatz P (n, t) ∼ e −S(n,t) into Eq. (S1) and neglecting terms of order O(n −1 ) we obtain a classical Hamilton-Jacobi equation for the action function S(n, t):

where H is the Hamiltonian and p = −∂ n S is the conjugate momentum. Instead of directly solving the Hamilton-Jacobi equations, we use the Hamilton approach for the classical equations of motion [Eq. (6)]

We write the action on a classical trajectory as [31]:

where the energy E ≡ H[n(t), p(t)], is constant along a dynamical trajectory given by p(n) = log q/(q − En β ) . To find the energy we solve the equation of motion (S19) on this given dynamical trajectory, which yieldṡ For n 1 and β > 0 the right hand side of Eq. (S21) varies very slowly with time [O(n −β )] (as shown below, the energy E also scales as n −β ), such that the solution for Eq. (S21) can be approximated as n = (qn −β − E)t + C. Here, C is a slowly-varying function of time, and includes constants such that the energy corresponding to the mean-field solution n = n obeys E(n = n ) = 0 [31]. Having shown that n ∼ t 1/(1+β) , we find C = n β/(1 + β), which indeed varies with time slower than t. Substituting this back into the equation for n and solving for the energy yields

where we have expressed t in terms of n and defined x ≡ n/ n . Substituting the energy (S22) into Eq. (S20) and solving the integral yields S(n, t) = n S(x) ,

where B(z; a, b) is the incomplete beta function, defined as B(z; a, b) = z 0 u a−1 (1 − u) b−1 du, and we define f (x) = 1 − x β (β(x − 1) + x). This result coincides with Eq. (8) of the main text and is valid in the limit of n 1.

To get better insight for the Gaussian vicinity of P n (t), we solve Eq. (S19) in the low energy limit E 1. This yields an approximated solution for the energy in the form

where we have again expressed t in terms of n . Equation (S24) is indeed small in the limit |x − 1| 1, which is the Gaussian vicinity of P n (t). By further approximating the dynamical trajectory as p(n) En β /q + E 2 n 2β /(2q 2 ), we substitute this back into Eq. (S20). Performing the integral and substituting Eq. (S24) yields the following action

which coincides with Eq. (9) of the main text. Equation (9) can be shown to solve Eq. (S18) in the limit |x − 1| 1. In Fig. 1a and S1a we compare the WKB solutions to simulations. Notably, while in Fig. 1a (main text) we are able to plot the exact results, in Fig. S1 , due to the larger values of n, the exact result cannot be plotted with standard computational tools. 

Here we analytically solve Eq. (10) of the main text in the limit of small variance σ. In the main text we write the probability of having visited n sites at time t as [Eq. (10)]

which can be numerically solved (Figs. 1b, S1b and S2). Analytical progress can only be made in limit of small variance, σ 1/ n 0 , where n 0 = [(1 + β 0 )qt] 1/(1+β0) is the mean number of sites given β = β 0 . For simplicity we focus on the small energy regime, yet similar calculations can be made with the full expression for the action [Eq. (8)]. Substituting Eq. (9) into Eq. (S26) yields

where S β (x) is given by Eq. (9). This integral can be solved for σ 1/ n 0 , using the saddle point approximation, which yields P (n, t) ∼ e − n 0 S β 0 (n/ n 0 )+ n 2 0 σ 2 S1(n/ n 0 ) ,

with

Here, the mean number of sites obeys n = n 0 1 + O(σ 2 ) , whereas the variance obeys

Thus, while inter-individual variability will almost not affect the mean number of sites, it does significantly affect the variance of the number of sites (by a factor of n 0 compared to that of the mean), see Fig. S2 .

Here we provide details on the mean-field equation for the mean number of sites. In particular we explicitly derive and solve this equation for all α values in the limit of t 1, and provide evidence of a phase transition at α = 1. Our  FIG. S3 . Comparison between the value of ξ as obtained from 100 simulation of length t = 10 5 (symbols) to the theoretical prediction (dashed lines). Plotted for β = 0.5 (blue crosses) and β = 1 (red triangles). In the inset we plot the relative error between the predicted value and the one obtained in simulations.

starting point is Eq. (11) of the main text

where P new and Π i are given by Eq.

(1) in the main text.

A. The case of α = 1

In the main text we assumed that P new → 0 and solved Eq. (S31). Here we provide a solution to the mean-field equation without neglecting P new . In mean field, we obtain an equation for the first moment m i by multiplying Eq. (S31) by m i and summing over all m i . For α = 1 this yields

where we a priori (to be checked a posteriori ) assume that j m j m i for any site i such that the denominator can be taken out of the sum over m i , and that n j=1 m j n j=1 m j . To find the value of n j=1 m j ≡ Q we sum over both sides of Eq. (S32) to obtain a differential equation for Q: ∂ t Q = (1 − q n −β ), an equation which is solved by Q = t − n . Substituting this back into Eq. (S32) gives

which is solved, assuming site i is first visited at time t i [i.e., with an initial condition m i (t i ) = 1], by [21]

Here n ti is the average number of sites at time t i , and on the right we approximated the solution for t 1 and discarded terms of order O(t 1/(1+β) ). This final result agrees with the one found in the main text. Importantly, as all sites have a linear dependence on t, we verify a posteriori that j m j m i for any site i, as contribution from all visited sites will not diminish at long times. For α < 1 we solve the mean-field equation at t 1 such that P new → 0,

where, similarly to the case of α = 1, we a priori assume that j m j α m i α for any site i, and we further assume n j=1 m j α = At ξ with α < ξ < 1 (to be proved a posteriori, see below). Notably, such a scaling was found to hold in numerical simulations. For initial condition m i (t = t 0 ) = 1, Eq. (S35) is solved by

Note that the asymptotic scaling of this result at t t i depends on the value of ξ, where for ξ < 1, Eq. (S36) predicts an asymptotic scaling of m i ∼ t (1−ξ)/(1−α) [1 + O(t ξ−1 )], for all sites. Now, as all sites scale similarly with t, it is evident that j m j α m i α , thus verifying our initial assumption. Using this solution for m i , we find that up to some unknown factor n j=1 m j α ∼ t 1/(1+β) t α(1−ξ)/(1−α) , entailing that ξ = α(1−ξ)/(1−α)+1/(1+β) = (1+αβ)/(1+β).

In Fig. S3 we show that this prediction agrees with simulations for two different values of β, up to a maximum of 3% relative error. However, we note that this relative error becomes crucial when ξ is substituted back into the scaling m i ∼ t (1−ξ)/(1−α) found above. Let us denote ξ 0 = (1 + αβ)/(1 + β) such that ξ = ξ 0 (1 − ), where 1 is a small correction that depends on α and β (see Fig. S3 ). Substituting ξ into m i ∼ t (1−ξ)/(1−α) one readily obtains (1 − ξ)/(1 − α) = β/(1 + β) + ξ 0 /(1 − α). Here, the approximation is valid only as long as β/(1 + β) /(1 − α), or alternatively 1 − α [assuming β = O(1)]. As we numerically find that = O(10 −1 ) this condition is hard to satisfy as α approaches 1, and in this regime it is preferable to find ξ directly from simulations.

In Fig. S4 we plot the number of visits to the five most visited site for different values of α. We show that for α < 1 all sites converge to the same number of visits, while for α > 1 the most visited site diverges from all other sites (see below). All α values show good agreement with the theory presented above. Similarly, in the bottom panels of S4 we show evidence for our a priori assumption that n j=1 m j α = At ξ , again with good agreement to the theory. 

The Egyptian fruit bat (Rousettus aegyptiacus, EFB) is a long-lived, widely distributed Old World fruit bat [53] . Like other fruit bats, individual EFBs tend to feed on a small subset of available trees and repeatedly revisit them for weeks and even months [11, 54] affirms that EFBs rely heavily on individual memory. Additionally, it has recently been shown that EFBs obtain a "cognitive map," which encompasses information about a large number of tree locations, suggesting that memory expands beyond the trees used at a given time [11] . Bats were tracked at a 0.125Hz sampling rate for an average tracking period of 23.7 nights and up to 131 nights. The data also includes nearly all fruit trees in the study area (14,314 trees and 18,111 orchard trees), which enabled us to identify specific tree visits.

To segment the data into movements and tree visits, we first filtered raw EFB tracks for localization errors based on the covariance matrices attributed to each ATLAS fix [55] . Localization that exceeded the highest realistic speed threshold for this species (20 m s ) were removed. Visits to trees were defined based on track segmentation using the firstpassage algorithm to determine the center of a "cloud of fixes" where the animal has spent a specified number of observations within a certain radius (for source code and details see https://github.com/ATLAS-HUJI/R/tree/master/ AdpFixedPoint). Finally, the median coordinates of each cloud were related to the closest tree in the dataset.

To make the seasonal classification most relevant for bats' foraging, we defined winter and summer based on the known peak of fruiting periods of the main seasonal tree species the bats frequently visit in the study area. These are the mulberry (Morus negra) and common fig (Ficus carica) species during May-September (summer) and Chinaberry (Melia azedarach) during November-February (winter). During each fruiting period we used 10 day periods for each bat to ensure sufficient statistics (many of the bats did not have longer tracks during a single season). Based on the field work it is reasonable to assume that during such 10 day periods no significant resource depletion occurs. To fit between the data and the theory we used a standard least square fit procedure from python scipy package.

Wkb theory of large deviations in stochastic populations

11) is similar to the master equation in Ref

However, they only consider the case of β = 0

To have proper statistics and ensure no significant resource depletion, we analyse 10-day periods for each bat. For details on the species, seasonality and the segmentation and fitting procedures, see SI Sec

Note that the power-law dependence for the mean has also been experimentally observed in humans [21] and monkeys [25]; however, these studies did not analyze fluctuations around n

Repeatability and heritability of exploratory behaviour in great tits from the wild

Quantifying individual variation in behaviour: mixed-effect modelling approaches

Stochastic optimal foraging: Tuning intensive and extensive dynamics in random searches

Composite random search strategies based on non-directional sensory cues

Bacterial persistence as a phenotypic switch

Finding disease outbreak locations from human mobility data

On estimating the predictability of human mobility: the role of routine

The universal visitation law of human mobility

Human mobility in response to covid-19 in France, Italy and UK

A general approach to detecting migration events in digital trace data

Empirical observations of ultraslow diffusion driven by the fractional dynamics in languages

Noise in transcription negative feedback loops: simulation and experimental analysis

Assaf, Dynamics of simple gene-network motifs subject to extrinsic fluctuations

Characterizing the accuracy of a self-synchronized reverse-gps wildlife localization system

Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables

Stochastic methods

Mathematical methods for physicists

Rousettus egyptiacus

Resource ephemerality drives social foraging in bats

A guide to pre-processing high-throughput animal tracking data

OV and MA were supported by the Israel Science Foundation Grant No. 531/20. MA was also supported by the Humboldt Research Fellowship for Experienced Researchers of the Alexander von Humboldt Foundation. Fieldwork with ATLAS was funded by the Minerva

with initial condition P (n, 0) = δ n,1 ⇒ G(n, 0) = z. Substituting Eq. (S1) into the definition of G(n, t) yields ∂G ∂t = q n z n (n − 1) −β P (n − 1, t) − z n n −β P (n, t) = q(z − 1) n z n n −β P (n, t) q(z − 1) n D −β z z n−β P (n, t) = q(z − 1)D −β z z −β G(z, t) , (S11)z a (z − ξ) β−1 f (ξ)dξ, and we used the following relationsHere, the equality on the left is valid for β ≤ 1 while the approximation on the right holds for n 1 and is exact in the special cases of β = 0, 1. Notably, for β < 0, i.e., when the growth rate increases in n, Eq. (S11) is a fractional differential equation for G, while for β > 0, i.e., when the growth rate decreases in n, it is a fractional integrodifferential equation. Although Eq. (S11) is hard to solve analytically for general β, and the direct method presented above is more suited, it can be solved for β = −1, 0, 1.

For β = 0, Eq. (S11) simplifies to ∂ t G(z, t) = q(z −1)G(z, t) and is accordingly solved by G(z, t) = ze −qt(1−z) . Using Eq. (S10) we find that P (n, t) follows a Poisson distribution (S9). In particular, in this case we have n = σ 2. Importantly, for α > 1 and β > 0 it follows that α > 1/(1 + β), such that m 1 α ∼ t α t 1/(1+β) ∼ n j=2 m j α , thus verifying our initial assumption. As discussed in the main text, these results suggest a phase transition at α = 1, see Fig. 3 and S5, and the next subsection.

Here we prove that there is a phase transition at α = 1, with no a priori assumptions on the solution (see previous sections). To this end, we define Q ≡ n i=1 m α i , so the probability that the most visited site will be visited again in the next time step t will be p 1 (t) = m α 1 /Q, and for any site in general we have p i (t) = m α i /Q. Next, we write an expression for the expected value of p 1 (t) in the next time step:where we have defined f α (m i ) ≡ (m i + 1) α − m α i . The first term on the right hand side represents the case for which the most visited site is visited in the next time step t + 1, and the second term corresponds to the case where any other site is chosen instead.In the case α = 1 we have f 1 (m i ) = 1 for any m i . Taking into account that p 1 (t) = m α 1 /Q, Eq. (S39) leads after some algebra to p 1 (t + 1) = p 1 (t) independently of the specific set of values {m i } we have. Thus, the probability of revisiting the most visited site will be kept constant through time (and the same can be proved for any other site).For α > 1 we note that f α (m i ) increases monotonically with m i . This, together with the fact that p 1 (t) = 1 − n i=2 p i (t) allow us to write the inequalityFinally, introducing p 1 (t) = m α 1 /Q into the previous inequality, after some algebra we obtainWe thus conclude that p 1 (t + 1) > p 1 (t) regardless of the specific set {m i } we have. The probability of revisiting the most visited site thus always increases with time on average, leading eventually to its dominance over the others. For α < 1 we proceed in a similar manner. Here, f α (m i ) will decrease monotonically with m i , so we can write p 1 (t + 1) < 1 − (f α (m 2 ) − f α (m 1 )(Q − m α 1 ) (Q + f α (m 1 ))(Q + f α (m 2 )) p 1 (t).This leads to p 1 (t + 1) < p 1 (t), such that for α < 1, on average the probability of revisiting the most visited site will decrease with time, thus leading to a much more homogeneous distribution of revisits among all sites available.