Evolutionary foraging in grids: Intermittent search dynamics emerge in finite, depletable landscapes
Authors: Shailendra Bhandari, Alex Szorkovszky, Anis Yazidi, Pedro G. Lind
Organizations: Department of Computer Science OsloMet – Oslo Metropolitan University St. Olavs plass, N-0130 Oslo, Norway · Simula Research Laboratory Numerical Analysis and Scientific Computing Oslo, 0164, Norway · Department of Informatics University of Oslo Gaustadalléen 23B, 0373 Oslo, Norway · Department of Technology Kristiania University of Applied Sciences Oslo, Norway
How search strategies evolve in finite, depletable landscapes remains a question in foraging theory. We study this problem with an evolutionary simulation in which agents forage on a two-dimensional toroidal lattice containing non-renewable resources distributed uniformly or as Lévy dust. Each agent carries a heritable genome encoding step lengths, velocities, and turning angles, and selection acts on a fitness function combining energetic gain, movement cost, and coverage efficiency. By allowing movement traits to evolve without imposing a prescribed power-law step-length distribution, we test whether evolved trajectories are better described by intermittent-search or Lévy-walk dynamics. Our results indicate that evolved search is more consistent with intermittent dynamics than with strict scale-free Lévy motion in the finite depletion-driven landscapes considered here. We characterize the dynamics by fitting second- and fourth-order displacement moments to intermittent-search and Lévy-walk models. While a Lévy-like random walk fits the evolutionary trajectories well (mean adjusted R2 > 0.9 in most tested conditions), intermittent search achieves a closer fit (mean adjusted R2 > 0.99) for all tested resource distributions. This preference holds across the tested grid sizes and resource densities. Five independent evolutionary runs per environment on a 503 x 503 grid at nominal resource density ρ = 0.15 reproduce this preference for the uniform environment and five Lévy-dust environments. Evolution rapidly reshapes the movement genome toward short displacements while retaining a sparse tail of longer relocations, consistent with local exploitation punctuated by occasional transfer. The framework provides a controlled setting for studying how search rules emerge under resource limitation and may inform resource-constrained exploration in autonomous systems.
Figures & tables
Figure 1: Resource landscapes and best evolved trajectories for the L=503 torus. Top row : placement-step distributions used to generate the resource fields, shown for the uniform control and for Lévy-dust environments with μ∈{1.0,1.5,2.0,2.5,3.0} . Main panels : resource locations (gray points) overlaid with the trajectory of the best evolved agent in each environment (colored line). Changing μ changes the spatial correlation structure of the resource field by altering the frequency of long placement steps: smaller μ gives heavier-tailed placement steps, whereas larger μ suppresses long placement steps. The inset in the μ=3.0 panel shows a magnified local patch.
Figure 2: Evolution of fitness and the step lengths on the 503×503 grid. Left: mean fitness over 1500 generations for each resource environment, with shaded standard errors across agents within a single evolutionary run at each generation. Right: distributions of step-length genome entries at generations 0, 200, and 1500. Selection rapidly concentrates the encoded step-length pool at short values while retaining a sparse tail of longer entries.
Figure 3: Moment-based comparison of evolved trajectories on the 503×503 grid. Empirical second moments m2 and fourth moments m4 are compared with best-fit IS and LW models. Columns correspond to the six resource environments. Across environments, IS gives a closer fit to the empirical moment curves.
Grid
Environment
Γ
RˉIS2
RˉLW2
211×211
uniform
0.028±0.004
0.9985±0.0004
0.970±0.004
211×211
μ=1.0
0.019±0.002
0.9984±0.0002
0.980±0.002
211×211
μ=1.5
0.022±0.003
0.9981±0.0003
0.976±0.003
211×211
μ=2.0
0.024±0.005
0.9993±0.0002
0.975±0.005
211×211
μ=2.5
0.030±0.007
0.9987±0.0005
0.969±0.007
211×211
μ=3.0
0.08±0.02
0.9986±0.0006
0.92±0.02
Table 1: Model-comparison summary across environments and grid sizes at ρ=0.15 . For each case, we report Γ=RˉIS2−RˉLW2 , together with the adjusted coefficients of determination. For L=211 and L=1009 , values are mean ± standard deviation over the final 500 best-agent trajectories from one evolutionary run. For L=503 , values are mean ± standard deviation across five independent run-level means, each computed from the final 500 best-agent trajectories. The runs used independent resource maps and GA initializations. Positive Γ indicates preference for the IS model. All results are reported with significant figures in accordance with their respective standard deviations.
Figure 4: Movement traits and return-interval statistics of the best evolved agents across grid sizes after 1500 generations. Columns correspond to resource environments. Top row: velocity-genome entries of the best evolved agents for L=211 , 503 , and 1009 . Middle row: turn-angle increments recorded along the best evolved trajectories. Bottom row: running maximum of return intervals Mjreturn versus return-event index j .
Appendix figures & tables10 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 5: Evolutionary stabilization on the 503×503 grid. Left: best fitness as a function of generation. Middle: JS divergence JS(Pg(s),P0(s)) between the population-level step-length distribution at generation g and the initial distribution. Right: final energetic and coverage efficiencies, ηE and ηC , for the best evolved agent in each environment. Vertical markers denote the operational plateau generation estimated from a rolling-slope criterion applied to fitness and genome redistribution.
Environment
ρ=0.15
ρ=0.30
ρ=0.45
Uniform
0.150
0.299
0.449
μ=1.0
0.138
0.256
0.359
μ=1.5
0.131
0.245
0.344
μ=2.0
0.117
0.221
0.310
μ=2.5
0.097
0.184
0.263
μ=3.0
0.066
0.145
0.194
Appendix
Table 3: Realized resource occupancy ρreal , the fraction of distinct resource-bearing sites before depletion, for the saved 503×503 maps at the three original nominal densities. Each entry describes the corresponding resource map.
Figure 6: Additional-size results for the 211×211 torus. Top left : mean fitness over 1500 generations for each resource environment, with shaded standard errors across agents within the single evolutionary run at each generation. Top right : distributions of genome-encoded step-length entries at generations 0, 200, and 1500. Middle : resource maps overlaid with the trajectory of the best evolved agent in each environment. Bottom : empirical second and fourth displacement moments compared with best-fit IS and LW models.
Figure 7: Additional-size results for the 1009×1009 torus. Top left : mean fitness over 1500 generations for each resource environment, with shaded standard errors across agents within the single evolutionary run at each generation. Top right : distributions of genome-encoded step-length entries at generations 0, 200, and 1500. Middle : resource maps overlaid with the trajectory of the best evolved agent in each environment. Bottom : empirical second and fourth displacement moments compared with best-fit IS and LW models.
Figure 8: Robustness of evolved search dynamics across resource densities on the 503×503 torus. Top : best evolved trajectories overlaid on the corresponding resource maps for the highest tested density, ρ=0.45 . Bottom left : final-population distributions of genome-encoded step lengths for ρ∈{0.005,0.010,0.15,0.30,0.45} , shown separately for each resource environment. Bottom right : model-selection score Γ across densities, with error bars.
Environment
ηC , baseline F=ηEηC
ηC , ablation F=ηE
Retained ηC (%)
Γ , ablation F=ηE
Uniform
0.9278±0.0065
0.9283±0.0048
100.1%
0.0190±0.0013
μ=1.0
0.9282±0.0057
0.9232±0.0055
99.5%
0.0180±0.0021
μ=1.5
0.8997±0.0053
0.8920±0.0103
99.1%
0.0216±0.0034
μ=2.0
0.8585±0.0083
0.8489±0.0073
98.9%
0.0313±0.0049
μ=2.5
0.8129±0.0076
0.8105±0.0112
99.7%
0.0373±0.0100
μ=3.0
0.7834±0.0118
0.7422±0.0154
94.7%
0.0509±0.0097
Appendix
Table 4: Coverage-term ablation at L=503 and ρ=0.15 . Values are mean ± standard deviation over the final 20 best-agent generations. They represent correlated late-generation best-agent evaluations from one evolutionary run per objective, rather than independent evolutionary replicates.
Environment
Γcueon
Γcueoff
RˉLW,cueoff2
Uniform
0.0175±0.0014
0.0010±0.0001
0.9988
μ=1.0
0.0185±0.0017
0.0029±0.0008
0.9969
μ=1.5
0.0204±0.0024
0.0059±0.0024
0.9933
μ=2.0
0.0326±0.0046
0.0145±0.0062
0.9818
μ=2.5
0.0251±0.0060
0.0267±0.0175
0.9720
μ=3.0
0.0424±0.0085
0.0797±0.0483
0.9181
Appendix
Table 5: Fixed-genome cue-off evaluation at L=503 and ρ=0.15 . Genomes evolved under the standard pcue=0.5 condition were re-evaluated with pcue=0 , without further evolution. The cue-on and cue-off Γ values are the mean ± standard deviation over the final 500 matched best-agent generations. The final column reports the mean adjusted LW coefficient of determination in the cue-off condition. These results are separate from the five-run aggregates in Table 1 .
Figure 9: Synthetic validation of the moment-based classifier and intermittent-search parameter estimation. Left: distributions of the model-comparison score Γ=Rˉ2IS−Rˉ2LW for 1000 synthetic L’evy-walk trajectories and 1000 synthetic intermittent-search trajectories. Positive values favor IS, whereas negative values favor LW. The confusion matrix shows that 999 of 1000 trajectories from each generating model were classified correctly. Right: estimated versus generating values of λBD for synthetic IS trajectories. The dashed line marks exact recovery. The inset summarizes the regression slopes and coefficients of determination for D , vB , λBD , and λDB .
Env.
F
ηE
ηC
Return
Hit
New-site
q0.95/mdn.
q0.95/mdn.
q0.95/mdn.
μ=1.0
0.40 ± 0.02
0.43 ± 0.01
0.92 ± 0.02
25.4
3.0
2.0
μ=1.5
0.39 ± 0.01
0.43 ± 0.01
0.90 ± 0.01
18.6
3.5
2.0
μ=2.0
0.40 ± 0.01
0.46 ± 0.01
0.86 ± 0.01
81.5
6.0
2.0
μ=2.5
0.37 ± 0.02
0.48 ± 0.02
0.81 ± 0.01
55.6
6.0
2.1
μ=3.0
0.37 ± 0.05
0.49 ± 0.04
0.76 ± 0.06
42.9
5.0
3.0
Appendix
Table 6: Waiting-time statistics for the best evolved agents. Fitness and efficiency values - cf. Eqs. ( 2 ) and ( 3 ) - are reported as mean ± standard deviation across grid sizes . Upper-tail spread is summarized by the median across grid sizes of q0.95/median . Return intervals quantify recurrence to previously visited sites, resource-hit intervals quantify the temporal spacing between resource encounters, and new-site discovery intervals quantify the timing of first visits to previously unvisited sites.
Figure 10: Empirical survival distributions of waiting-time observables in evolved search trajectories. Left : return intervals Δtreturn , measuring the time between successive visits to the same previously visited site. Middle : resource-hit intervals Δthit , measuring the time between successive resource encounters. Right : new-site discovery intervals τn , measuring the waiting time between first visits to successive previously unvisited sites.