Paper deep dive
Tuning Agent-Based Predator-Prey Models Toward Lotka-Volterra Dynamics
Corinna Mandl, Siddharth Chaturvedi, Marcel van Gerven
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 7/18/2026, 4:05:55 PM
Summary
This paper presents a method for tuning agent-based predator-prey models to exhibit classical Lotka-Volterra population dynamics. Using the ABMax framework, the authors simulate sheep and wolves as active agents with continuous-time recurrent neural network (CTRNN) controllers. They optimize ecological parameters (e.g., metabolic rates, birth/death thresholds) using a feature-based loss function that rewards sustained oscillations, phase lag, and bounded populations, effectively bridging microscopic agent rules with macroscopic Lotka-Volterra cycles.
Entities (8)
Relation Signals (6)
ABMax → implements → Agent-Based Model
confidence 95% · The model is implemented in ABMax, a JAX-based agent-based modelling framework
Lotka-Volterra → istargetfor → Agent-Based Model
confidence 94% · tune the microscopic parameters of a continuous predator-prey ABM so that its population trajectories satisfy LV-like constraints
CTRNN → controls → Sheep
confidence 93% · sheep and wolves are active agents with ... recurrent neural network-based controllers
CTRNN → controls → Wolf
confidence 93% · sheep and wolves are active agents with ... recurrent neural network-based controllers
CMA-ES → optimizes → ABM-level ecological parameters
confidence 91% · used covariance matrix adaptation evolutionary strategy (CMA-ES) to minimise Equation (14) over ABM-level parameters
JAX → powers → ABMax
confidence 90% · ABMax, a JAX-based agent-based modelling framework
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Recent growth in compute power has made it increasingly feasible to use large-scale agent-based models to simulate complex adaptive systems. A central difficulty is that such models contain many local rules and parameters, where small changes can lead to runaway behaviour, population collapse, or saturation at artificial bounds. We study this problem in a continuous predator-prey system where sheep and wolves are active agents with local sensing, internal energy, and recurrent neural network-based controllers. We ask whether environmental and demographic parameters can be tuned so that the resulting population dynamics resemble classical Lotka-Volterra cycles. We optimise these parameters with a feature-based loss that rewards sustained oscillations, phase lag, bounded populations, and long-term persistence, first for random controllers and then for evolved controllers in a more naturalistic setting. The model is implemented in ABMax, a JAX-based agent-based modelling framework that enables efficient batched simulation on hardware accelerators.
Tags
Links
- Source: https://arxiv.org/abs/2606.13639v1
- Canonical: https://arxiv.org/abs/2606.13639v1
PDF not stored locally. Use the link above to view on the source site.
Full Text
35,513 characters extracted from source content.
Expand or collapse full text
11institutetext: Department of Machine Learning and Neural Computing Donders Institute for Brain, Cognition and Behaviour, Radboud University Nijmegen, the Netherlands 11email: siddharth.chaturvedi@donders.ru.nl Tuning Agent-Based Predator-Prey Models Toward Lotka-Volterra Dynamics Corinna Mandl These authors contributed equally. Siddharth Chaturvedi0 Marcel van Gerven Abstract Recent growth in compute power has made it increasingly feasible to use large-scale agent-based models to simulate complex adaptive systems. A central difficulty is that such models contain many local rules and parameters, where small changes can lead to runaway behaviour, population collapse, or saturation at artificial bounds. We study this problem in a continuous predator-prey system where sheep and wolves are active agents with local sensing, internal energy, and recurrent neural network-based controllers. We ask whether environmental and demographic parameters can be tuned so that the resulting population dynamics resemble classical Lotka-Volterra cycles. We optimise these parameters with a feature-based loss that rewards sustained oscillations, phase lag, bounded populations, and long-term persistence, first for random controllers and then for evolved controllers in a more naturalistic setting. The model is implemented in ABMax, a JAX-based agent-based modelling framework that enables efficient batched simulation on hardware accelerators. 1 Introduction Agent-based models (ABMs) connect local interactions among individual agents with population-level dynamics, making them useful for studying complex adaptive systems in ecology, artificial life, neuroscience, and economics [11, 6]. However, their flexibility also makes them difficult to control. Many ABMs contain nonlinear interactions whose aggregate outcomes can be multi-modal, skewed, or fat-tailed. In such settings, rare large events are more likely than under thin-tailed statistics, making model dynamics sensitive to local rules, initial conditions, and parameters [9, 5]. A central challenge is therefore to constrain ABM dynamics without removing the interactions that make them interesting. Predator-prey systems provide a useful testbed for this problem. They are simple enough to analyse, but rich enough to produce sustained population-level structure through ecological feedback. Sheep support wolves, wolves suppress sheep, and the decline of one population changes the growth conditions of the other. This feedback loop is captured by the classical Lotka-Volterra (LV) equations [21, 12], which provide a principled macroscopic target for stabilising an agent-based predator-prey model. Recent work has increasingly treated the relation between microscopic rules and macroscopic behaviour as a modelling problem in its own right. Macroscopic models can be derived from microscopic swarm descriptions and used to guide individual rules [18], while ABMs can be calibrated from data by fitting parameters or latent micro variables to reproduce aggregate patterns [15]. In predator-prey systems, ABMs have been reduced to coarse-grained stochastic dynamics [17], and LV models and ABMs have been used as complementary descriptions of the same process [10]. Our work follows this micro-macro tradition, but uses the desired macroscopic regime as an explicit top-down optimisation target. Rather than reducing an ABM to a coarse LV-like model, we tune the microscopic parameters of a continuous predator-prey ABM so that its population trajectories satisfy LV-like constraints. We model sheep and wolves as self-propelled active agents moving in a continuous two-dimensional arena [24]. The interaction rules are inspired by the canonical Wolf Sheep predation model [23]. Agents also have a continuous-time recurrent neural network (CTRNN) controller to control their motion. The optimisation is carried out in two separate stages. First, we optimise agent controllers that determine sensory-motor behaviour. Second, we freeze either evolved or random controllers and optimise ABM-level ecological parameters controlling energy intake, predation, metabolism, birth, and death. The rest of the paper describes the model, optimisation procedure, results, and limitations. 2 Model The model consists of circular active agents moving in a two-dimensional Cartesian plane (Fig. 1a). It is inspired by related ABMax models of active particles and ecological agent-based systems [3, 4]. Agents do not undergo hard-body collisions and can pass through each other. The system is integrated with a forward Euler scheme of step size Δt t, with time evaluated on t∈kΔt:k=0,1,…,Tt∈\k t k=0,1,…,T\. Figure 1: (a) Snapshot of model with agent interactions. (b) Different sensing scenarios. (c) Energy during the birth-death process of agents. 2.1 Population turnover The agent population =∪A=W\,∪\,S is divided into sheep S and wolves W, representing prey and predators. To preserve static array shapes during JAX compilation, each species is represented by a fixed-size buffer. The model accommodates at most NmaxN^S_ sheep and NmaxN^W_ wolves, of which N(t)N^S(t) and N(t)N^W(t) are active at time t. Population turnover follows threshold-based birth and death. Inspired by the canonical Wolf Sheep model [23], an agent i∈i is eliminated when its internal energy ei(t)e_i(t) remains below a death threshold ede_d for an uninterrupted interval tdt_d. In contrast to stochastic reproduction in the canonical model, an active agent becomes viable for reproduction when ei(t)e_i(t) remains above a birth threshold ebe_b for an uninterrupted interval tbt_b. Reproduction occurs only if an inactive slot of the same species is available. The offspring inherits the parent’s controller parameters, receives half of the parent’s energy, and is initialised near the parent by perturbing its Cartesian coordinates within [−ds,ds][-d_s,d_s] along each axis. The parent energy is also halved, introducing a natural delay before the next birth (Fig. 1c). 2.2 Position model The position of the centre of an active agent is i(t)∈ℝ2q_i(t) ^2, and its orientation is θi(t)∈[−π,π) _i(t)∈[-π,π). Positions are clipped component-wise to the rectangular arena. Dropping the agent index for readability, the agent positions are updated as ˙(t) q(t) =(t),(t)=s(t)[cosθ(t)sinθ(t)], =v(t), (t)=s(t) bmatrix θ(t)\\ θ(t) bmatrix, (1) θ˙(t) θ(t) =ω(t),ω(t)=u(t). =ω(t), ω(t)=u(t). (2) Here, (t)∈ℝ2v(t) ^2 is the Cartesian velocity and ω(t)∈ℝω(t) is the angular velocity. The translational and rotational velocity commands si(t)s_i(t) and ui(t)u_i(t) for an agent i are obtained from the controller readouts as si(t) s_i(t) =k(s)tanh(ai(1)(t))(1+ϵξi(s)(t)), =k^(s) \! (a_i^(1)(t) ) (1+ε _i^(s)(t) ), (3) ui(t) u_i(t) =k(u)tanh(ai(2)(t))(1+ϵξi(u)(t)), =k^(u) \! (a_i^(2)(t) ) (1+ε _i^(u)(t) ), (4) where ϵε is a positive scalar by which Gaussian noise terms ξ(s)(t),ξ(u)(t)∼(0,1)ξ^(s)(t),ξ^(u)(t) (0,1) are multiplied. The translational and angular speeds are scaled by positive constants k(s)k^(s) and k(u)k^(u), respectively. The quantities ai(1)(t),ai(2)(t)∈ℝa_i^(1)(t),a_i^(2)(t) are the two raw readouts of the continuous-time recurrent neural network controller. 2.3 Energy model Each active agent maintains an internal energy ei(t)e_i(t). Sheep gain energy from a circular grass patch, wolves gain energy by consuming sheep, and both species lose a constant metabolic cost. For sheep i∈(t)i (t), grass intake is ei(G)(t)=k(G)δ(G)(i(t)),e^(G)_i(t)=k^(G)δ^(G)(q_i(t)), (5) where k(G)k^(G) is the grass intake rate and δ(G)(⋅)δ^(G)(·) is 11 when the sheep lies within the grass patch of radius dGd_G centred at the origin, and 0 otherwise. Predation is a local energy transfer from sheep to wolves. Each wolf consumes at most one sheep at a time, namely the closest active sheep within its capture distance d_A. We define the capture indicator as a Boolean-step function given by δji(P)(t)δ^(P)_ji(t) as δji(P)(t)=1,if i∈(t) is nearest to j∈(t) and ‖i(t)−j(t)‖2≤d,0,otherwise.δ^(P)_ji(t)= cases1,&if i (t) is nearest to j (t) and \|q_i(t)-q_j(t)\|_2≤ d_A,\\ 0,&otherwise. cases (6) When sheep i is caught by wolf j, a fraction k(P)k^(P) of the sheep’s current energy is transferred. Thus, the predation loss of sheep i is ei(P−)(t)=k(P)ei(t)∑j∈(t)δji(P)(t),e^(P-)_i(t)=k^(P)e_i(t) _j (t)δ^(P)_ji(t), (7) and the predation gain of wolf j is ej(P+)(t)=k(P)∑i∈(t)ei(t)δji(P)(t).e^(P+)_j(t)=k^(P) _i (t)e_i(t)δ^(P)_ji(t). (8) This allows multiple wolves to target the same sheep, in which case the sheep loses the corresponding energy fraction once for each wolf that catches it. The energy updates for sheep and wolves become ei(t+Δt) e_i(t+ t) =ei(t)+ei(G)(t)−ei(P−)(t)−μ(),i∈(t), =e_i(t)+e^(G)_i(t)-e^(P-)_i(t)-μ^(S), i (t), (9) ei(t+Δt) e_i(t+ t) =ei(t)+ei(P+)(t)−μ(),i∈(t). =e_i(t)+e^(P+)_i(t)-μ^(W), i (t). (10) The internal energy of agents is further clipped between emine_ and emaxe_ . If an agent’s energy remains below the death threshold for the required duration, it is removed from the active population. If it remains above the reproduction threshold for a certain duration, and a free slot is available, it reproduces, as described in Section 2.1. 2.4 Lotka-Volterra objective In this subsection, we explain the empirically engineered loss function that is used to tune the discrete population dynamics of the predation ABM according to smooth Lotka-Volterra dynamics. The chief ingredient of the loss function is the Pearson correlation function [14]. For two equal-length time series a(t)a(t) and b(t)b(t), the Pearson correlation is defined as C(a,b)=⟨(a(t)−⟨a⟩)(b(t)−⟨b⟩)⟩(σ(a)+ϵ(C))(σ(b)+ϵ(C))C(a,b)= (a(t)- a )(b(t)- b ) (σ^(a)+ε^(C) ) (σ^(b)+ε^(C) ) (11) where ⟨⋅⟩ · denotes the temporal mean, σ(a)σ^(a) and σ(b)σ^(b) are the temporal standard deviations of a(t)a(t) and b(t)b(t), and ϵ(C)ε^(C) is a small positive constant used for numerical stability. The correlation ranges from −1-1 (anti-aligned) to 11 (perfectly aligned). In our case, it is used as ℒ(corr) ^(corr) =(1−C(x¯(t),y¯(t+ℓ)))+(1−C(Δx¯(t),Δy¯(t+ℓ))) = (1-C( x(t), y(t+ )) )+ (1-C( x(t), y(t+ )) ) +12(1−C(x¯(t)−⟨x¯⟩,Δy¯(t)))+12(1−C(y¯(t)−⟨y¯⟩,−Δx¯(t))) + 12 (1-C( x(t)- x , y(t)) )+ 12 (1-C( y(t)- y ,- x(t)) ) (12) where x¯(t) x(t) and y¯(t) y(t) are the smoothed versions of the normalised sheep and wolf population time series, respectively, with x(t)=N()(t)Nmax(),y(t)=N()(t)Nmax().x(t)= N^(S)(t)N^(S)_ , y(t)= N^(W)(t)N^(W)_ . (13) The smoothing is performed using a uniform moving-average kernel of width m. Further, Δx¯(t)=x¯(t+Δt)−x¯(t) x(t)= x(t+ t)- x(t), and similarly for Δy¯(t) y(t). We also compare the sheep population at time t with the wolf population shifted by a delay ℓ . The rationale is that in a classical Lotka-Volterra curve, the sheep population leads the wolf population, because an increase in prey availability supports predator growth only after a delay, while the subsequent increase in predators suppresses the prey population. Thus, the first two terms in ℒ(corr)L^(corr) reward the lagged alignment of prey and predator population levels and their slopes. The last two terms encode the local Lotka-Volterra pressure, namely that above-average sheep abundance should be associated with wolf growth, while above-average wolf abundance should be associated with sheep decline. Other components of the main loss function discourage degenerate solutions that satisfy the correlation loss ℒ(corr)L^(corr) without showing sustained Lotka-Volterra-like oscillations. These include the extinction loss ℒ(e)L^(e) and the ceiling loss ℒ(c)L^(c), which penalise population trajectories going below and above prescribed thresholds, respectively. The former prevents population collapse, while the latter prevents saturation at the artificial population capacity of the simulation. Next, the amplitude loss ℒ(a)L^(a) requires both populations to show sufficiently large oscillations in both halves of the rollout. The turning loss ℒ(t)L^(t) requires both upward and downward movement in each half of the trajectory, thereby rejecting monotonic growth or collapse. The crossing loss ℒ(cross)L^(cross) penalises trajectories that cross their temporal mean fewer than n(cross)n^(cross) times. Finally, the drift loss ℒ(d)L^(d) penalises long-term changes in the mean population level by comparing the beginning and end of the rollout. Thus, the final loss function for tuning the ABM is a weighted sum of all the loss components described above, ℒ=∑r∈λ(r)ℒ(r),=corr,e,c,a,t,cross,dL= _r λ^(r)L^(r), =\corr,e,c,a,t,cross,d\ (14) where each λ(⋅)λ^(·) is a positive weighting constant. With the agent controllers fixed, we used covariance matrix adaptation evolutionary strategy (CMA-ES) [8] to minimise Equation (14) over ABM-level parameters only. The ABM parameter vector subjected to optimisation is given by ϕ(ABM)=[k(G),k(P),μ(),μ(),eb(),ed(),tb(),td(),eb(),ed(),tb(),td()],φ^(ABM)= [k^(G),k^(P),μ^(S),μ^(W),e_b^(S),e_d^(S),t_b^(S),t_d^(S),e_b^(W),e_d^(W),t_b^(W),t_d^(W) ], (15) where k(G)k^(G) is the grass intake rate, k(P)k^(P) is the predation transfer rate, μ()μ^(S) and μ()μ^(W) are the metabolic costs of sheep and wolves, and eb,ed,tb,tde_b,e_d,t_b,t_d denote the energy and time thresholds for birth and death for each species. In practice, CMA-ES samples M raw candidate vectors at each generation. Each raw vector is transformed with a sigmoid function and rescaled to a predefined viable range before being evaluated in the ABM. For every candidate vector, one ABM rollout is simulated for T(ABM)T^(ABM) steps, after which the resulting sheep and wolf population trajectories are assigned a loss according to Equation (14). CMA-ES then updates its search distribution using these losses and samples a new set of candidate vectors. This process is repeated for G(ABM)G^(ABM) generations. 2.5 Sensor model Agents emit r equiangular rays from i(t)q_i(t), covering the full angular range with maximum length dℛd_R. Ray-agent intersections are computed using the standard ray-circle quadratic test [7]. If multiple active agents are intersected, only the closest intersection is retained, making sensing occlusion-aware and producing blind spots (Fig. 1b). Each ray p∈1,…,rp∈\1,…,r\ returns two channels of information. The first channel is the distance from the ray origin to the closest intercepted agent surface. The second channel is a type value indicating whether the intercepted agent is a sheep or a wolf. Let τj∈1,−1 _j∈\1,-1\ denote the type value of an agent j, with sheep encoded as 11 and wolves encoded as −1-1. The observation of ray p emitted by agent i is given by i,p(t)=[di,p(t),τj]⊤,if ray p first intersects active agent j,[dℛ,0]⊤,if no active agent is intersected.o_i,p(t)= cases [d_i,p(t), _j ] ,&if ray p first intersects active agent j,\\ [d_R,0 ] ,&if no active agent is intersected. cases (16) Here, di,p(t)d_i,p(t) is the distance along the ray to the first valid intersection, and τj(t) _j(t) is the type of agent j intersected by the ray. The external observation vector is obtained by concatenating all ray observations, i(ex)(t)=[di,1(t),τi,1(t),di,2(t),τi,2(t),…,di,r(t),τi,r(t)]⊤.o^(ex)_i(t)= [d_i,1(t), _i,1(t),d_i,2(t), _i,2(t),…,d_i,r(t), _i,r(t) ] . (17) The agent controllers also receive a subset of their internal states as input to make them aware of the proprioceptive variables. Hence, the ray observations are concatenated with some internal variables. Let ηi(t) _i(t) denote the pre-metabolic energy gain or loss of agent i, and let fi(t)∈0,1f_i(t)∈\0,1\ denote whether agent i is overlapping with at least one other active agent. The internal observation vector is i(in)(t)=[ei(t),ηi(t),fi(t),‖i(t)‖2,i⊤(t),ωi(t),θi(t)]⊤.o^(in)_i(t)= [e_i(t), _i(t),f_i(t),\|q_i(t)\|_2,v_i (t), _i(t), _i(t) ] . (18) Thus, the complete observation vector provided to the controller is i(t)=[i(ex)(t),i(in)(t)]⊤.o_i(t)= [o^(ex)_i(t),o^(in)_i(t) ] . (19) 2.6 Controller model The agents control their motion using velocity controllers, which are modelled as continuous-time recurrent neural networks (CTRNNs). Specifically, we use a modern interpretation of the Wilson-Cowan model [20]. For an agent i, the hidden state dynamics are given by ˙i(t)=(tanh(Ji(t)+Ei(t)+)−i(t))⊙(cτσ()), z_i(t)= ( \! (Jz_i(t)+Eo_i(t)+b )-z_i(t) ) (c_τσ( τ) ), (20) [ai(1)(t),ai(2)(t)]⊤=Di(t). [a_i^(1)(t),\,a_i^(2)(t) ] =Dz_i(t). (21) Here i(t)∈ℝhz_i(t) ^h is the controller hidden state, interpreted as neuronal population activity [20]. The recurrent matrix is J∈ℝh×hJ ^h× h, the observation matrix is E∈ℝh×(2r+8)E ^h×(2r+8), ∈ℝhb ^h is a bias, ∈ℝh τ ^h contains learnable time-scale parameters, and D∈ℝ2×hD ^2× h maps hidden states to readouts. The controller learnable parameter set is ϕ()=J,E,,,Dφ^(A)=\J,E,b, τ,D\. 2.7 Agent controller optimisation Before ABM-level tuning, we optimise separate sheep and wolf controller parameters, ϕ()∗φ^(S)* and ϕ()∗φ^(W)*, using the same-rollout evaluation trick used in [3]. First, to train the wolves, we sample different controller parameter vectors ϕi()φ^(W)_i, for all i∈i , from a CMA-ES distribution. We then evaluate all wolves against sheep whose controller parameters are randomly initialised as ϕj()∼U(−1,1)φ^(S)_j U(-1,1), for all j∈j . The fitness for each wolf at the end of a rollout is given as fi=ei(T)−ei(0)f_i=e_i(T)-e_i(0). The fitness assigned to each CMA-ES sample is averaged across H scenarios having different initial conditions and ABM parameters ϕ(ABM)φ^(ABM). Finally, CMA-ES optimises the wolf controller parameters over G()G^(W) generations. During this optimisation phase, the number of agents is held constant, such that N()(t)=Nmax()N^(W)(t)=N^(W)_ and N()(t)=Nmax()N^(S)(t)=N^(S)_ . The sheep controller is optimised in the same way by training sheep against randomly initialised wolf controllers. We use random opponents to reduce overfitting to a single adversarial strategy. After optimisation, the mean controller parameters of the respective CMA-ES distributions are used as the fixed sheep and wolf controllers in the Lotka-Volterra tuning phase. 2.8 Simulation details The entire model was implemented using the JAX [1] based ABMax framework [2], which enables batched vectorisation across multiple instances of the model. Specifically, the rank-match update algorithm from the framework allowed different model instances to maintain different numbers of active agents while being evaluated in parallel. We adopted the CMA-ES algorithm implemented in the Evosax library [13]. The values of all parameters and initial states are summarised in Table 1. Unless otherwise stated, these values are used for all simulations. All simulations were run on an NVIDIA A100 GPU. The code for the model is available online.111https://github.com/corimandl/LV_ABM.git Table 1: Model parameters and initial conditions Symbol Value Symbol Value Symbol Value ABM and initial conditions Nmax()N^(S)_ 400400 Nmax()N^(W)_ 400400 N()(0)N^(S)(0) 150150 N()(0)N^(W)(0) 8080 Δt t 0.10.1 (qxmax,qymax)(q_x ,q_y ) (500,500)(500,500) R()R^(S) 55 R()R^(W) 55 dGd_G 200200 dsd_s 3030 emaxe_ 100100 i(0)q_i(0) ([−250,250]2)U([-250,250]^2) i(0)v_i(0) 0 θi(0) _i(0) [−π,π)U[-π,π) ωi(0) _i(0) 0 ϵε 0.050.05 k(s)k^(s) 2R/Δt2R/ t k(u)k^(u) 2.02.0 Lotka–Volterra objective m 5151 ℓ 600600 ϵ(C)ε^(C) 10−810^-8 (ρlow,ρup,ρceil)( _low, _up, _ceil) (0.05,0.88,0.97)(0.05,0.88,0.97) AminA_ 0.200.20 Δmin _ 0.120.12 n(cross)n^(cross) 33 δdrift _drift 0.120.12 λ(ϕ)λ^(φ) 3535 λ(e)λ^(e) 12001200 λ(c)λ^(c) 10001000 λ(a)λ^(a) 500500 λ(t)λ^(t) 800800 λ(cross)λ^(cross) 5050 λ(d)λ^(d) 500500 T(ABM)T^(ABM) 60006000 Optimised ABM parameter ranges k(G)k^(G) [0.02,0.50][0.02,0.50] k(P)k^(P) [0.45,0.90][0.45,0.90] μ()μ^(S) [0.006,0.090][0.006,0.090] μ()μ^(W) [0.008,0.090][0.008,0.090] eb()e_b^(S) [0.30,0.60]emax[0.30,0.60]e_ ed()e_d^(S) [0.005,0.10]emax[0.005,0.10]e_ tb()t_b^(S) [0.003,0.020]T[0.003,0.020]T td()t_d^(S) [0.001,0.040]T[0.001,0.040]T eb()e_b^(W) [0.30,0.60]emax[0.30,0.60]e_ ed()e_d^(W) [0.005,0.10]emax[0.005,0.10]e_ tb()t_b^(W) [0.003,0.040]T[0.003,0.040]T td()t_d^(W) [0.001,0.020]T[0.001,0.020]T Sensing and controller r 1313 dℛ()d_R^(S) 300300 dℛ()d_R^(W) 300300 τ −1,0,1\-1,0,1\ i(t)o_i(t) ℝ2r+8R^2r+8 h 6060 cτc_τ 1010 ϕ()φ^(A) J,E,,,D\J,E,b, τ,D\ i(0),i(0)z_i(0),a_i(0) ,0,0 CMA-ES settings M 1616 G(ABM)G^(ABM) 150150 M()M^(S) 5050 M()M^(W) 5050 G()G^(S) 40004000 G()G^(W) 40004000 H 1212 T(),T()T^(S),T^(W) 1000,10001000,1000 ϕrand()φ^(S)_rand (−2,2)U(-2,2) ϕrand()φ^(W)_rand (−1,1)U(-1,1) fif_i ei(T)−ei(0)e_i(T)-e_i(0) αelite,σ0 _elite, _0 default 3 Results The results of the ABM-level parameter estimation can be seen in Fig. 2a. During this stage, agent controllers are fixed and only the ABM parameters ϕ(ABM)φ^(ABM) are optimised. We optimise for two fixed-controller settings: random controllers ϕ()∼U(−1,1)φ^(A) U(-1,1), and evolved controllers ϕ()∗,ϕ()∗φ^(S)*,φ^(W)* obtained from the controller-optimisation stage. The loss decreases in both settings, indicating that the proposed objective can guide ABM parameters toward lower-loss regimes. However, the loss decreases faster and converges to a lower value for evolved controllers. To make optimisation possible for random controllers, the grass patch spans the entire arena, so that randomly moving sheep have a constant energy supply. Figure 2: Lotka-Volterra tuning in the predator-prey ABM. (a) ABM-level loss during optimisation of ϕ(ABM)φ^(ABM) with fixed random or evolved controllers. Lines show means across four seeds and bands show one standard deviation. (b) Population trajectories under random ABM parameters. (c) Population trajectories after ABM-parameter optimisation. Fig. 2b shows the population trajectories over a rollout of T=15000T=15000 steps using randomly sampled ABM parameters. In both random and evolved controller settings, the resulting trajectories miss key characteristics of Lotka-Volterra dynamics, in particular sustained coupled oscillations between sheep and wolves. After optimisation using the procedure in Section 2.4, the population trajectories acquire oscillatory predator-prey structure and remain active for the duration of the rollout (Fig. 2c). Interestingly, the optimised random-agent setting produces smoother and more regular oscillations in the upper panel of Fig. 2c, but these oscillations also approach the artificial population ceiling. In contrast, the evolved-agent setting produces noisier and less idealised oscillations, but the trajectories remain more bounded and avoid the most extreme saturation events. This suggests that the loss can induce LV-like macroscopic structure in both settings, while learned controllers produce a more behaviourally grounded and less idealised population dynamics. This distinction is consistent with the fact that the classical Lotka-Volterra model is highly idealised and omits several ecological features present in spatial, individual-based predator-prey systems, such as finite resources, spatial structure, behavioural heterogeneity, and individual-level stochasticity [22, 19, 16]. Next, we rendered the optimised ABM using the evolved controllers from the separate controller-optimisation stage (Fig. 3a). The snapshots in Fig. 3b and c show that learned controller dynamics induce clear spatial organisation. Wolves aggregate in and around the grass patch, where sheep density is higher, while sheep remain more broadly dispersed. This behaviour was already present before ABM-parameter optimisation, suggesting that ABM tuning mainly makes this behavioural ecology viable over longer rollouts by balancing energy intake, predation, reproduction, and death. A video is available online.222https://youtu.be/vBS8EfF8GZ0 Figure 3: Agent controller optimisation and spatial organisation. (a) Sheep and wolf fitness curves, shown as mean across three seeds with one-standard-deviation bands. (b) Optimised ABM rollout with evolved agents at t=0t=0. (c) Same rollout at t=1000t=1000, showing wolf aggregation near the grass patch and broader sheep dispersion. 4 Discussion In this work, we showed that a continuous agent-based predator-prey model can be tuned such that its population trajectories show signatures of Lotka-Volterra-like dynamics. The implication is not that the ABM discovers Lotka-Volterra dynamics automatically, nor that a unique microscopic model can be inferred from a macroscopic trajectory. Rather, the result shows that a desired macroscopic ecological regime can be used as a top-down constraint for searching in the ABM parameter space. In this sense, Lotka-Volterra dynamics are used not as a replacement for the ABM, but as a target that identifies parameter regimes in which the ABM remains dynamically well behaved. This also demonstrates a use case for ABMax, where many ABM instances with different numbers of active agents can be evaluated in parallel on hardware accelerators. The results further suggest a distinction between behavioural structure and population-level stabilisation. The evolved controllers produced spatial predator-prey behaviour, with wolves aggregating near the grass patch and sheep spreading more broadly. This behaviour was not created by the Lotka-Volterra loss itself. Rather, the ABM-level tuning made the consequences of such behaviour viable over longer horizons by balancing energy intake, predation, reproduction, and death. Thus, our goal is different from showing that predator and prey agents optimised only for local behaviour necessarily produce Lotka-Volterra dynamics. Instead, we ask whether, once local behaviour exists, ecological parameters can be tuned so that the resulting population dynamics satisfy a desired macroscopic constraint. A limitation is that the loss function is indirect and hand designed. However, a direct mean-squared error to a single Lotka-Volterra trajectory would also be difficult to justify because exact trajectory matching is unrealistic. Future work can make the model more ecological and adaptive. For instance, internal energy could affect mass, speed, turning radius, or metabolic cost, and the grass patch could become a depleted and regenerating resource rather than a fixed energy source. Another direction is to allow predator and prey controllers to adapt online, creating a moving target for both populations. More broadly, this work points toward a way of studying open-ended multi-agent systems through ecological constraints. Instead of allowing rich agent interactions to drift without structure, macroscopic ecological targets may help keep the system within viable bounds while still preserving non-trivial dynamics. credits 4.0.1 Acknowledgements This publication is part of the project Dutch Brain Interface Initiative (DBI2) with project number 024.005.022 of the research programme Gravitation which is (partly) financed by the Dutch Research Council (NWO). 4.0.2 The authors have no competing interests to declare that are relevant to the content of this article. References [1] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang (2018) JAX: composable transformations of python+numpy programs. Note: http://github.com/jax-ml/jax Cited by: §2.8. [2] S. Chaturvedi, A. El-Gazzar, and M. van Gerven (2025) Abmax: a jax-based agent-based modeling framework. External Links: 2508.16508 Cited by: §2.8. [3] S. Chaturvedi, A. El-Gazzar, and M. van Gerven (2025) Emergence of internal state-modulated swarming in multi-agent patch foraging system. ArXiv preprint arXiv:2510.18886. Cited by: §2.7, §2. [4] S. Chaturvedi, A. El-Gazzar, and M. van Gerven (2026) Role differentiation in a coupled resource ecology under multi-level selection. ArXiv preprint arXiv:2604.00810. Cited by: §2. [5] A. Clauset, C. R. Shalizi, and M. E. J. Newman (2009) Power-law distributions in empirical data. SIAM Review 51 (4), p. 661–703. External Links: Document Cited by: §1. [6] D. L. DeAngelis and W. M. Mooij (2005) Individual-based modeling of ecological and evolutionary processes. Annual Review Of Ecology, Evolution, And Systematics 36, p. 147–168. External Links: Document Cited by: §1. [7] C. Ericson (2004) Real-time collision detection. Crc Press. Cited by: §2.5. [8] N. Hansen (2016) The cma evolution strategy: a tutorial. ArXiv preprint arXiv:1604.00772. Cited by: §2.4. [9] D. Helbing (2012) Agent-based modeling. In Social self-organization: Agent-based simulations and experiments to study emergent social behavior, p. 25–70. Cited by: §1. [10] M. Hodzic, S. Selman, and M. Hadzikadic (2016) Complex ecological system modeling. Periodicals Of Engineering And Natural Sciences 4 (1), p. 44–50. External Links: Document Cited by: §1. [11] J. H. Holland (1992) Complex adaptive systems. Daedalus 121 (1), p. 17–30. Cited by: §1. [12] T. Idema (2005) The behaviour and attractiveness of the Lotka-Volterra equations. Ph.D. Thesis, Universiteit Leiden. Cited by: §1. [13] R. T. Lange Evosax: jax-based evolution strategies. In Proceedings of the Companion Conference On Genetic And Evolutionary Computation, p. 659–662. Cited by: §2.8. [14] J. Lee Rodgers and W. A. Nicewander (1988) Thirteen ways to look at the correlation coefficient. The American Statistician 42 (1), p. 59–66. Cited by: §2.4. [15] C. Monti, M. Pangallo, G. De Francisci Morales, and F. Bonchi (2023) On learning agent-based models from data. Scientific Reports 13, p. 9268. External Links: Document Cited by: §1. [16] L. K. Mühlbauer, A. D. Schulze, W. S. Harpole, and A. T. Clark (2020) GauseR: simple methods for fitting Lotka-Volterra models describing Gause’s “struggle for existence”. Ecology And Evolution 10 (24), p. 13275–13283. External Links: Document Cited by: §3. [17] J. Niemann, S. Klus, and C. Schütte (2021) Data-driven model reduction of agent-based systems using the koopman generator. PLOS ONE 16 (5), p. e0250970. External Links: Document Cited by: §1. [18] Q. Quan, X. Yu, Y. Li, and G. Qi (2026) Macroscopic modelling and analysis based on microscopic models for swarm systems. Scientific Reports 16, p. 10342. External Links: Document Cited by: §1. [19] T. Räz (2017) The volterra principle generalized. Philosophy Of Science 84 (2), p. 209–229. External Links: Document Cited by: §3. [20] D. Sussillo (2014) Neural circuits as computational dynamical systems. Current Opinion In Neurobiology 25, p. 156–163. Cited by: §2.6, §2.6. [21] V. Volterra (1926) Fluctuations in the abundance of a species considered mathematically. Nature 118 (2972), p. 558–560. External Links: Document Cited by: §1. [22] P. J. Wangersky (1978) Lotka-Volterra population models. Annual Review Of Ecology And Systematics 9, p. 189–218. External Links: Document Cited by: §3. [23] U. Wilensky (1997) NetLogo wolf sheep predation model. Center for Connected Learning and Computer‐Based Modeling, Northwestern University, Evanston, IL. External Links: Link Cited by: §1, §2.1. [24] Y. Y. Yaya (2023) A predator–prey model from a collective dynamics and self-propelled particles approach. Computer Sciences & Mathematics Forum 7 (1), p. 50. External Links: Document Cited by: §1.