Exploration, movement, and search for resources are ubiquitous among organisms in nature [1–3]. Classical theories of search [4], such as optimal foraging theory [5, 6], have mostly focused on long time limits and typically assume that natural selection favors search strategies that maximize long-term encounters with nutrients. However, many phenomena in ecology [7] and other fields of biology operate in transient regimes [8, 9], extending over time scales that never reach the asymptotic stationary state [10]. Another typical assumption is to consider random walks [11–13] or diffusion processes [14] to describe the movement of searchers navigating the landscape based on local information [15–17]. Yet, in many instances, searchers can obtain and store [18] non-local information gathered through sensory cues [19, 20] or through anticipation of environmental changes [21–23] (Figure 1). The question then arises as to how such finite perceptual range can influence both the dynamics of movement and the search efficiency over the finite time scales relevant in biology.

Here we study the role of finite time scales associated with ecological movement and search; specifically, the effect of limited spatial perception when the search time is itself finite. To formalize these aspects, we propose an optimal navigation (ON) model, which allows us to extend the description of search as a biased random walk [11, 14, 24], and reinterpret it in the framework of optimal control theory. The ON model includes a time horizon that quantifies the perceptual range of the searchers along their trajectory and fixes a non-local optimization target for the agents. In the limit of vanishing time horizon (i.e., as the spatial perception shrinks and the information becomes local), the ON model recovers the classic Keller-Segel (KS) drift-diffusion model [25] of local search strategies (i.e., with instantaneous sensing and alignment to the point-wise gradient).

Using simulations and analytical results, we find that a population of non-local searchers moving toward a nutrient patch exhibits distinct transient behavior, clustering faster at the hotspot than local searchers, thereby increasing their search efficiency. Our results show that the maximum efficiency gain occurs when the perceptual range of the searchers matches the environmental length scale over which the nutrient concentration changes significantly. As the search time becomes asymptotically large or small, the efficiency gain from the non-local strategy diminishes, and the searchers behave effectively as local responders. If the environmental length scale changes, we show that the efficiency gain can be maintained as long as searchers can adjust their sensitivity. This means that finite perception remains advantageous to searchers that can rescale their response dynamically, or to populations that contain a diversity of responses. Finally, we consider the effect of interaction between searchers, and show that non-local information consistently reinforces the dominant strategy in the population and leads to improved search efficiency overall, even though multimodality (subpopulations) can appear during the transients. Our framework provides and optimization perspective on a range of collective phenomena in population biology and, more generally, on biologically-inspired search and exploration algorithms, thus shedding light on the role of spatial perception on finite-time search.

A classic model for the dynamics of a population of searchers using local gradient alignment is given by the Keller-Segel equation, which we briefly recap here.

Consider a population of searchers moving in a closed, bounded, d-dimensional domain Ω ⊂ ℝ^d. The searchers move by responding to two concentration fields: to a primary stimulant S₁(x) (e.g., nutrient) with sensitivity χ₁, and to a secondary stimulant S₂(x, t) (e.g., pheromone), released by the agents themselves, with sensitivity χ₂.

On long space and time scales relative to the microscopic motion, one can describe the biased random walk x(t) ∈ Ω of searchers by a Langevin equation [24]:

where ξ(t) is a white noise process, D is the coefficient of diffusion, and

is the velocity of the searcher. The parameters D, χ₁ and χ₂ are typically inferred experimentally from trajectories of the agents [14, 26], and can sometimes be expressed in terms of microscopic parameters [27]. Note that Equations (1) – (2) describe a searcher that uses local information, since it aligns its velocity instantaneously to the gradients of S₁ and S₂.

The secondary stimulant S₂ introduces interaction between the searchers. If S₂ is assumed to diffuse faster than the searchers, its evolution is given by

where

and Γ is the Gamma function (see Appendix A).

Taking all together, the time evolution of the population density ρ(x, t) of searchers obeying (1)– (3) can be described with a Fokker-Planck equation [28] known as the Keller-Segel (KS) model. In dimensionless variables x → x/L and t → t/T_D, where TD=L2/D is the diffusion time of the searchers, the KS equation reads

where the parameters Pe₁ = χ₁S_1,av/D and Pe₂ = χ₂S_2,av/D are Péclet numbers quantifying the ratio of diffusive to advective forces on the searchers, and S_1,av, S_2,av are the average stimulant concentrations. Given an initial distribution ρ₀(x): = ρ(x, 0), a stimulant profile S₁(x), and parameters Pe₁ and Pe₂, Equation (5) can be solved numerically using standard techniques. Here, we use a first-order in time, second-order in space forward Euler scheme [29] with upwind discretization and Δx = 0.01 and Δt = 10⁻⁶. We denote the solution of the KS model by ρ_KS(x, t). See Appendix A for a fuller derivation and details of the non-dimensionalization.

The variational rewriting (11) and its approximation scheme (12) leads us to formulate the optimal navigation (ON) search model, as follows. Consider a population of searchers that move by performing the optimization (12) over a finite time horizon τ ≥ Δt > 0, which reflects the perceptual range of the agents. Then the time evolution of the population corresponds to a sequence of constrained optimization problems [33], i.e., a succession of JKO solutions, each over time τ.

Starting from the initial density ρ₀(x), we construct the evolution of ρ(x, t), such that each iteration k ≥ 0 finds m(x, s): = ρ(x, kτ + s) for s ∈ [0, τ] by solving the minimization problem:

Note that the constraint is the continuity equation ensuring the conservation of ρ as in (6), whereas the cost function J contains a transportation cost, which constrains the average motion to geodesics between optimal states, and an end-point term involving the evaluation of the free energy at τ (10). Although we use it here for a particular form of the free energy functional, the formulation is generic: through suitable choice of F, the ON model (13) converges to a broad class of conservation laws as long as they can be recast as continuity equations and possess a variational structure [30–32, 34].

We denote the solution of the ON model (13) by ρ_ON(x, t; τ), and compute it using Algorithm 1, a gradient descent algorithm inspired by Burger et al. [35] and presented in detail in Appendix B.

Physically, the ON model (13) describes the motion of searchers that optimize their displacement over paths bounded by the time horizon τ. From the proof of the JKO scheme [32], it follows directly that the ON model recovers the local KS model as τ → 0:

For finite horizon τ > 0, the time evolution of the ON model departs from the KS solution due to the effect of non-local information on the movement of the searchers, as explored below.

We use the ON model (13) to study how the finite perception of the agents (encapsulated in the time horizon τ > 0) affects the search at the population level. We first consider non-interacting searchers insensitive to the secreted stimulant, i.e., Pe₂ = 0. The case of interacting searchers is presented in section 6.

Our numerics start with a uniform initial condition ρ₀(x) = 1 and we compute ρ_ON(x, t; τ), the time evolution (13) of the ON population of non-local searchers with time horizon τ > 0. We also compare it to the time evolution (5) of a KS population of local searchers, or equivalently the ON model with τ = 0.

Both the KS and ON models converge to the same stationary solution ρ_∞(x) asymptotically as t → ∞:

for all τ and S₁, given by the Gibbs-Boltzmann distribution (15), where Z is a normalization constant. This result is well known for the KS equation [28]. To check that ρ_∞(x) is also the stationary solution of the ON model, note that at stationarity dW(ρ∞,ρ′)=0 in (12), and the result (16) follows from solving for the minimizer.

The approach to stationarity, on the other hand, reveals differences between the KS and ON models. To develop our intuition, let us first consider a linear attractant gradient S₁(x) = αx. Such a landscape is fully characterized by the local gradient ∇S₁(x), and hence non-local searchers have no advantage since τ > 0 provides no further information than what is known by local searchers. In line with this expectation, when we solve the ON model (13) using the method of characteristics (section S2 in the SI), we find that the drift velocity predicted by the KS and ON models coincide for all τ, i.e.,

This observation implies that non-local search is only advantageous in landscapes with non-zero curvature. To illustrate this point in more detail, let us consider a static Gaussian concentration of stimulant with characteristic length scale σ ≪ 1 over the domain Ω = [−1/2, 1/2]^d (Figure 2A):

where det(·) is the matrix determinant and x^T is the transpose of the vector x. For simplicity of the computations that follow, we take Σ=σ2Id, where I_d is the d-dimensional identity matrix. Such a Gaussian landscape serves as a simple model of a stimulant patch emanating from a point source, and its characteristic length scale σ indicates regions of steep attractant gradients near the source (||x|| ≪ σ) and regions of shallow gradients in the tails of the distribution (||x|| ≫ σ).

In Figure 2 we show that, as the time horizon τ is increased, the population of ON searchers exhibits a faster approach to stationarity, as measured by the normalized L²-distance between ρ_ON(x, t; τ) and ρ_∞(x) as a function of time for different values of τ:

Figures 2A– C presents the one-dimensional case (d = 1). The effect of the horizon in accelerating the convergence to stationarity is shown in Figure 2B. Note also that for small values of τ, an intermediate, quasi-steady distribution develops during the transient (e.g., τ = 10⁻⁵ in Figures 2B,C). This long-lived intermediate behavior is the result of the population evolving on two timescales [14]: searchers near the maximum of S₁(x) (||x|| ≪ σ) are driven by advection due to the steep gradient, whereas those far from the maximum (||x|| ≫ σ) are driven by diffusion in shallow gradients, and hence move more slowly toward the maximum. Due to the slow diffusive searchers, the stationary state is only reached at a longer time scale t ~ 1. As the horizon τ is increased, this dual behavior (diffusion- or advection-dominated) is lost: the searchers escape quickly the diffusion-dominated part of the domain and, as a result, the distribution approaches stationarity increasingly faster with no appreciable quasi-steady transient distribution (e.g., τ = 10⁻³ in Figures 2B,C). The same behavior is observed also in the two-dimensional case (d = 2) in Figures 2D–F.

Such transient states can be important in biological systems, which typically operate on time scales far from the asymptotic long-time regime [7–9]. In our setting, this situation arises when the search time T (which is analogous to the foraging effort in ecology) is smaller than the diffusion-dominated convergence time, i.e., when T ≪ 1. In such a situation, non-local (ON) searchers have an advantage over local (KS) searchers since they converge faster to areas with high concentration of attractant. To quantify this effect, we consider the amount of stimulant S₁ encountered over the search time T

and define the relative search efficiency as

where U^(0,T) is the uptake of the population of KS searchers. Therefore, U(τ, T) > 1 indicates a gain in search efficiency, that is, increased stimulant encountered by the population due to the perceptual horizon τ > 0.

Our numerics in Figure 3A show that, given a finite search time T, the search efficiency (19) reaches a maximum U∗(τσ∗,T) for searchers operating with an optimal horizon τσ∗, which depends on the length scale σ for a given dimension d (Figure 3A). The presence of a maximum follows from the asymptotic behavior U(τ, T) → 1 for τ → 0 and τ → ∞. The latter limit follows from the invariance of ρ_∞(x) under τ, and the fact that the integral (19) is asymptotically dominated by the steady state. The presence of a maximum is also observed in the two-dimensional case, but U∗(τσ∗,T) decreases in higher dimensions.

The dependence of τσ∗ with the length scale of the landscape σ obtained numerically from our simulations is shown in Figure 3B with solid (1D) and open (2D) circles. To understand this dependence, consider a searcher at x(t) obeying the Langevin equation (1) under the ON model. The reachable set until t + τ is within a ball of radius x_ms(τ) defined by

where the two terms represent the displacement due to diffusion and to an effective drift velocity v̄ over time τ, respectively. To estimate v̄(τ), we assume that the search time T is large enough such that individual searchers have explored the whole domain (e.g., for T = 0.1 this is fulfilled when τ > 10⁻³ as seen in Figure 2B). The effective drift velocity can then be approximated by the velocity of an average searcher (over the domain) that maximizes its gain up to time τ,

where the integral is evaluated along a one-dimensional cross-section of the Gaussian landscape

From (20) and (21), we obtain an estimate of the horizon τ necessary to search over a distance x_ms with the ON model.

The relevance of this estimate is shown in Figure 3B, which shows that for small σ, the maximum search efficiency is attained when the mean-squared displacement of the searchers equals the length scale of the environment:

as obtained with our approximation. The ON model thus predicts that the most efficient searchers are those that tune their horizon such that they traverse the characteristic length scale of the environment within one optimization step. Shorter or longer optimizations lead to a decreased search efficiency.

The dependence of this behavior on the dimension d is also captured by (20), which tells us that the ballistic and diffusive terms balance when v̄(τ)2τ/(2d). Thus for a given S₁, τ and Pe₁ the motion of searchers becomes gradually diffusion dominated as the dimension d increases. As a result, the relative search efficiency decreases in higher dimensions as shown in Figure 3.

The search efficiency of the ON model depends on the length scale σ of the Gaussian landscape: the ON gain diminishes as σ increases (Figure 3). However, as we now show, an ON searcher can retain the same search efficiency under a Gaussian landscape with varying length scale by adjusting Pe₁, its sensitivity to the stimulant.

To see this, consider an ON searcher starting at x₀ exposed to its nutrient micro-environment until time T. The effective gradient for this searcher depends on the starting position x₀ and is given by:

For a fixed exploration time T, an increase in the stimulant length scale σ leads to shallower effective gradients (Figure 4A). Using asymptotic techniques, we show in section S1 in the SI that the effective gradient (23) for symmetric Gaussian profiles has a well-defined behavior in the two limiting regimes:

Together with the form of the dynamics (1), this suggests the following scaling for the Péclet number:

with d ≤ α(T, d) ≤ 2 + d.

We have tested this scaling by obtaining the ON solution (13) over a given T for Gaussian profiles with different σ using the renormalized Péclet number (25). We then compute the relative search efficiency (19) for this solution, U(τ,T;Pe~1(σ,T)). Our numerics in Figure 4B show that the search efficiency curves for the renormalized parameter (25) for different σ collapse on a single curve:

The exponent α(T, 1), i.e., for d = 1, is obtained numerically (Figure 4C) is consistent with the expected asymptotic limits (24). Note that, the effective gradient (23) is a function of x_ms(τ), which changes with dimension. Therefore, although for any d a scaling relationship exists, it will be different from the curve in Figure 4C depending on d.

Hence the search efficiency can be made invariant for different environmental length scales σ by rescaling the Péclet number (25). Alternatively, adjusting α(T, d) can be viewed as responding to a “renormalized landscape” [Pe~1S1(x)] in order to maintain the ON search efficiency. This is intuitive in limiting cases: when the search time is small (T → 0), the efficiency remains unchanged on landscapes with similar local gradients near the center (x₀ ≪ 1); when the search time is large (T → ∞), the efficiency is invariant for landscapes with similar effective gradient over the whole domain (see inset of Figure 4C).

This result suggests that searchers can optimize their search efficiency by adjusting their response sensitivity [as in the scaling (25)] so as to balance the relative effect of the advection and diffusion velocities or, in other words, the relative importance of gradient optimization vs. noisy exploration. Since the diffusion coefficient D is typically independent of S₁(x) [27], the adjustment of Pe₁ could be achieved by varying the sensitivity as a function of the stimulant, i.e., χ₁(S₁(x)) (see Appendix A). In the Discussion, we explain possible biological mechanisms to achieve this effect.

Is non-local search advantageous over search times relevant for ecological systems? The invariant search efficiency U(τ,T) characterizes the performance of a searcher that is tuned to the intrinsic length scale of the stimulant landscape during its search time T. In Figure 4D, we show the dependence of the maximum renormalized search efficiency (26) with the search time T. As expected, for short search time T → 0 and long search time T → ∞, the efficiency U∗(T) is equivalent to the local search strategy, i.e., U∗(T)→1. However, between both extremes, searchers benefit from finite perception. Our numerics show that the optimal search time is T∗=argmaxU∗(T,τ∗)~0.1<1=TD. Hence finite perception is maximally advantageous for search times smaller than the diffusion time, a fact that is typical in ecological systems (section 7).

Until now, we have considered a population of non-interacting searchers reacting only to an external stimulant S₁(x). Now we consider searchers which also interact with each other through an attractive secondary stimulant S₂(x, t) (i.e., Pe₂ > 0). The effect of such interaction can be potentially contradictory: agent interaction may increase aggregation; however, aggregation might not increase search efficiency for the population if a large proportion of the agents clumps away from the stimulant source if the sensitivity to stimulant S₂ becomes larger.

To explore these effects numerically, we use the ON model (13) with a weak interaction in the free energy (10)

and compute the time evolution of the population for various β and τ. We restrict ourselves to the weak interaction case (β ≪ 1) to facilitate our numerics. Specifically, to prevent infinite density concentration at finite time, a well-known artifact of the KS model, by applying a small regularizing factor ωρ² with 0 < ω ≪ 1 to the free energy F[ρ] (10). This is a volume exclusion term that models the fact that agents occupy a finite volume in space [36, 37] [see (B5) in Appendix B]. The weak interaction assumption (β ≪ 1) is important for numerical performance so as to avoid high velocities near the concentration point that would necessitate a finer time discretization to maintain the CFL condition (see Appendix B).

A summary of the numerics of the two-dimensional ON case with different horizons τ and levels of interaction β over different search times T is presented in Figure 5. In all cases, the presence of interaction reduces the tails of the population density and increases aggregation near the maximum of S₁ at the center of the domain. This becomes more noticeable as the search time grows and we approach stationarity (Figure 5, right column). During transients, however, agent interaction induces multimodality in the population (Figure 5, left and middle column). This implies that some searchers move away from the maximum of the stimulant S₁(x) and aggregate into transient subpopulations. This behavior arises due to the non-linear response of the searchers to the gradient of S₁(x): for steep gradients of S₁(x), agents are driven by attraction to S₁, whereas for shallow gradients of S₁(x) the agents are driven by the interaction with the rest of the population through the secreted stimulant S₂. In contrast, for non-interacting agents (β = 0) the distribution is always uni-modal (Figure 5).

Figure 6 shows the dependence of multimodality on T and τ at three different interaction strengths β = 5 × 10⁻⁵, 10⁻⁵, 5 × 10⁻⁶. We find that interaction introduces multimodality, as part of the population clumps away from the source of the attractant. For higher interaction strength (β = 5 × 10⁻⁵) multimodality is only present during early transients, and is independent of the time horizon τ, since the strongly interacting population rapidly converges to the peak of the primary stimulant. As the interaction strength β decreases, multimodality lasts longer, but it decreases as the time horizon τ increases. This is in line with the expectation that increased spatial perception τ or search time T leads to overlapping information of the searchers about the environment, whereas, in contrast, when τ and/or T are small, the searchers remain isolated within their local environment, thus leading to multimodality by local aggregation. As the interaction decreases further multimodality disappears.

Importantly, despite the presence of such mildly multimodal transients with clumped subpopulations, agent interaction β leads to an overall improvement in the search efficiency U^*(τ^*, T) at a shorter time horizon τ^* (Figure 7). This behavior follows from the increased concentration of searchers around the center (on average), and suggest that when resources are sparsely distributed, agent interaction in conjunction with a finite perceptual range, could play a role in improving the collective sensing of the environment toward improved population efficiency.

In this work, we studied how finite perception influences the dynamics and efficiency of collective random search in a population of agents. Using concepts from optimal control and random walks, we proposed a model that encapsulates the spatial information the searchers possess as a time horizon for an optimization problem. Simulations of the dynamics of population search show that non-local information affects the movement strategy, as compared with the standard Keller-Segel model based on local optimization. Although non-local search does not change the stationary state, it leads to qualitatively different transient responses of possible relevance in biological systems [7]. For example, marine bacteria have been observed to aggregate at food patches much faster (~ 10¹− 10² s) [38, 39] than the timescale to reach the steady state distribution (TD~104-105 s, based on D ~ 10²−10³μm²/s [40] and typical inter-patch distance L ~ 10³μm [41]). Similarly, in several rodent species diffusion coefficients of D ~ 200 m²/day and home range of L ~ 70 m have been reported [42]. Therefore, the time it takes to reach home by diffusion (~ 25 days) is much longer than their typical response time (~ 1 day).

We find that non-local information can account for the increase in the search efficiency under transient search times: the maximum efficiency is reached when the mean-squared displacement of the searchers matches the environmental length scale of the stimulant. When the time horizon vanishes or when the search time is infinite, our model recovers the response of local searchers. This is in accordance with the fact that when long-range cues are unreliable, local response leads to highest efficiency [15].

We also showed that the search efficiency can be made invariant to changes in environmental length scales by suitably scaling the response sensitivity. As a consequence, a searcher with a given perceptual range may always achieve its maximum efficiency by dynamically adjusting its sensitivity to the environmental stimuli. This can be achieved either by dynamically rescaling the responses via adaptation at the agent level [43], or by the presence of a distribution of sensitivities among the agents at the population level. For example, it has been shown that phenotypic heterogeneity (or plasticity) across a population can be used to achieve maximum search efficiency in patchy environments [44].

Finally, we considered the effect of interaction between searchers with finite perception and showed that interaction can lead to unimodal or multimodal population distributions on transient timescales. Multimodality appears even in the presence of unimodal stimulant landscapes due to a trade-off between following the environmental gradient or the rest of the population. In our numerics, interaction always improved the overall search efficiency of the population.

Recent theoretical and experimental studies [45, 46] also suggest that rodents (and other higher animals) store spatial information of their environment, and in doing so commonly prefer longer (non-trivial) future paths as opposed to paths leading to immediate rewards. These results are in line with a predictive optimization based on spatial knowledge, as in the ON model, where the searcher weighs up local cues with those at a distance to inform the planning of the future trajectory as opposed to an immediate (gradient) optimization. Clearly, the perceptual horizon in the ON model is a very coarse abstraction of the spatial information the organism perceives (for example in its visual field) or is encoded in its neural or cellular memory. Going beyond our framework would require explicit assumptions about the sensory or cognitive mechanisms involved in search and navigation, an area of interesting future research which is out of the scope of this work.

Our work opens up several directions of research. Beyond our simple setup, it would be of interest to study search on temporally-fluctuating or patchy nutrient landscapes [47] using non-local strategies. Random search theory based on local response predicts that when a searcher is positioned equally far from two nutrient patches it is equally likely to explore either patch. However, on transient timescales, nonlocal searchers are expected to explore the patch with denser resources with higher probability. This direction will be the object of future work.

Movies accompanying Figure 2F and a Matlab implementation of the ON model in 1D and 2D is available at: https://doi.org/10.6084/m9.figshare.7477568.v1.

AG and MB conceived and developed the study. AG performed the research under the supervision of MB and JC. AG and MB wrote the manuscript, with inputs from JC.