Learning Individual Behavior in Agent-Based Models with Graph Diffusion Networks
Authors: Francesco Cozzi
Marco Pangallo
Alan Perotti
André Panisson
Corrado Monti
Venue: NeurIPS 2025 poster ()
Keywords: Agent-Based Models
Surrogate models
Diffusion Models
Graph Neural Networks
Probabilistic Machine Learning
Primary area: machine_learning_for_sciences
TL;DR: Developed a model to reproduce stochastic individual agent dynamics in agent-based models

Abstract: Agent-Based Models (ABMs) are powerful tools for studying emergent properties in complex systems. In ABMs, agent behaviors are governed by local interactions and stochastic rules. However, these rules are ad hoc and, in general, non-differentiable, limiting the use of gradient-based methods for optimization, and thus integration with real-world data. We propose a novel framework to learn a differentiable surrogate of any ABM by observing its generated data. Our method combines diffusion models to capture behavioral stochasticity and graph neural networks to model agent interactions. Distinct from prior surrogate approaches, our method introduces a fundamental shift: rather than approximating system-level outputs, it models individual agent behavior directly, preserving the decentralized, bottom-up dynamics that define ABMs. We validate our approach on two ABMs (Schelling's segregation model and a Predator-Prey ecosystem) showing that it replicates individual-level patterns and accurately forecasts emergent dynamics beyond training. Our results demonstrate the potential of combining diffusion models and graph learning for data-driven ABM simulation.

Learning Individual Behavior in Agent-Based Models with Graph Diffusion Networks
Francesco Cozzi
Marco Pangallo CENTAL, Turin, Italy marco marco. warrooppngalll@centtt.cco
Sapienza University, Rome, Italy CENTAL, Turin, Italy f. cozzi Suniromal....... it
Alan Perotti C CENTAI, Turin, Italy
André Panisson C CENTAI, Turin, Italy
Corrado Monti @ CENTAI, Turin, Italy
alan. wwt.perrooticentai.....
panisson@centai.cm
me@corradomonti. com
Abstract
Agent-Based Based Models (ABMs) are powerful tools for studying emergent properties in complex systems. In ABMs, agent behaviors are governed by local interactions and stochastic rules. However, these rules are, in general, non-diiferentiable, which limits the use of gradient-based methods for optimization, and thus integration with real-world data. We propose a novel framework to learn a differentiable surrogate of any ABM by observing its generated data. Our method combines diffusion models to capture behavioral stochasticity and graph neural networks to model agent interactions. Distinct from prior surrogate approaches, our method introduces a fundamental shift: rather than approximating system-level outputs, models individual agent behavior directly, preserving the decentralized, bottom-up dynamics that define ABMs. We validate approach on two ABMs (Schelling's segregation model and a Predator- Prey ecosystem) showing that it replicates individual-level patterns and accurately forecasts emergent dynamics beyond training. Our results demonstrate the potential of combining diffusion models and graph learning for data-driven ABM simulation.
Introduction
Agent-Based Models (ABMs) are computational frameworks in which autonomous "agents" interact with each other and their environment, leading to emergent collective behavior [43]. ABMs are typically characterized by: (i) well-defined network of interactions, where the state each agent influenced by the states of specific set of other agents, usually from the previous time step; (ii) stochasticity, meaning that agents' decisions incorporate a degree of randomness, producing probability distributions over multiple runs that capture real-world uncertainty and variation. ABMs have proven to be powerful tool for developing and refining theoretical understanding, particularly in identifying minimal of micro-level rules that generate realistic macro-level outcomes [34]. In this sense, they have been instrumental in modeling a diverse range of phenomena [6]. including structure formation biological systems, pedestrian traffic, urban aggregation, and opinion dynamics. More recently, ABMs have demonstrated their value as forecasting tools [32], such as in predicting the economic impacts of the COVID-19 pandemie [30].
However, this progress is occurring despite the absence of principled methods to systematically align ABMs with real-world data. While various approaches have been proposed for calibrating mecro-level parameters ABMs [31], there are still no established methods for tuning the micro-level behaviors
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Training
ABM
1+1
Simulate
ABM-Generated Simulation
Train
1-1
Differentiable Surrogate
Generation
state
Sample future states
z
7
Differentiable Surrogate
GDN 4.00
7,00
(Z ee
Figure I: Overview of the training and generation pipeline for differentiable surrogates of Agent- Based Models. The top-left panel illustrates the training process: we run simulations using the original ABM, and use the resulting trajectories to train the differentiable surrogate, The top-right panel shows the structure of the ABM-generatee data using the Predator-Prey model as an example: a state at time step 1-1 gives rise to multiple possible states at time I, one which chosen to generate further possible states at 1++1. Colored cells highlight the behavior of specific "prey" agent green "move," for "die," and pink for "reproduce." The bottom panel shows the generation phase: given new observed state, the trained surrogate simulates plausible future states.
of individual agents to match observed data. One potential approach is to manually construct probabilistic model that replicates the ABM and then use its likelihood function estimate individual state variables [25]. However, this method requires the manual development of an ad ad-hoc probabilistic framework that reproduces the original ABM. Thus, what is currently missing is fully automated method for deriving a learnable, differentiable model directly from an ABM.
In this work, we propose novel approach to address this challenge: combining a graph neural network with a diffusion model learn differentiable surrogate of ABM, from ABM ABM generated data. We refer to this method as Graph Diffusion Network (GDN). In this framework, graph neural network captures the interactions that govern the evolution of each agent's state response to other agents, while the diffusion model learns the distribution of possible state transitions, conditioned on these interactions. central aspect of our approach is its explicit modeling of individual agent behavior. Rather than treating the system independent entity, while also incorporating the influence of other agents on its decisions. This approach ensures that the emergent dynamics remain faithful to the decentralized nature of ABMs. By constraining the surrogate to model only micro-level rules, it cannot rely on shortcuts that predict only macro-level outcomes, preserving the distributed logic of the original ABM.
whole, we focus
how each agent acts
Our approach draws inspiration from previous work on using neural network models to emulate deterministic cellular automata [13]. However, we extend this idea to the broader domain of ABMs by introducing a crucial component: stochasticity. By incorporating stochasticity, our architecture can learn directly from ABM-ggnnrated data traces, making it adaptable to wide variety of agent-based

---

models across diverse real-world applications. Furthermore, since our method is trained on data traces, can seamlessly integrate empirical observations alongside simulated data, thus being potentially applicable to real-world scenarios. In this sense, our work represents first step toward developing comprehensive methodology for creating easy-to-use, learnable ABMs.
2
Background
From general perspective, an ABM can be represented as stochastic process Z. ~ Pe Zr<c), where Z, denotes the state variables at time I, θ is set of parameters, and is a probability distri- bution implicitly defined by the model structure and parameters. The index l represents discrete time. Typically, @ consists of a small number of parameters, remains fixed in dimension. is interpretable by domain experts, and serves the model's primary control mechanism. Conversely, each element in Z, captures an agent's 'ss state, leading to a nigh-dimensional state space.
To illustrate this structure, we consider two ABMs used throughout the paper. The first is the well- known model by Schelling [36], where Zμ captures agents' positions and colors, and indicates preference for same-color neighbors. Even with some tolerance for neighbors of different colors, agents often form segregated clusters [41]. This clear mismatch between individual preferences and aggregate outcomes is a classic example of emergence. The second model is a predator-prey model [39, 43], describing the ecological dynamics between two interacting species, with one acting as predator and the other as prey, similarly to the Lotka-Volterra equations. In this ABM, Z, includes agent position and type (prey-predde prrey-predator, while θ governs the probability to move, reproduce, or die. This model replicates the cyclical predator-pree -rrey population dynamics, typical of Notka-Volteres systems, while also capturing complex spatial patterns reminiscent of spatial evolutionary games [27]. Both the Schelling and predator-prey models are widely recognized as canonical ABMs and are frequently used as testbeds for the development of novel calibration and surrogate techniques [37, 16, 23, 39, 40, 26).
ABMs have traditionally been powerful for theory generation, but recent years, they have become increasingly data-driven [29]. To align ABM output with empirical data, most efforts focus on calibrating parameters θ so that model-generated summary statistics match observed ones [31, 33]. Less attention, however, has been paid to estimating agent states Z, which is key for matching time series further to summary statistics. Some researchers use data assimilation methods like particle filters [22] ensemble Kalman filters [28] to infer Z, A more principled alternative is to make ABMs differentiable, enabling the maximization of a likelihood function via gradient descent and automatic differentiation [24, 25]. While differentiability is straightforward for simple stochastic behaviors, such as those governed by Bernoulli trials [3], it becomes far more challenging for complex behaviors like those observed in the Schelling and predator- prey models.
To address this and other challenges in ABMs, researchers have increasingly turned to more tractable surrogates, also known meta- models emulators [17, 9, 11]. Surrogate models typically learn directly the mapping from parameters θ to static summary statistics, disregarding individual behavior and model dynamics. For instance, surrogate in Lamperti et al. [20] maps θ to the mean growth rate of the economy. More recent research has also explored the emulation of model dynamics. Grattarola et al. [13] use Graph Neural Networks to approximate cellular automata, which can be seen as a special case of ABMs with deterministic interaction rules. Dyer et al. [8] propose Ordinary Differential Equation emulators to construct interventionally consistent surrogates, ensuring that micro-state interventions produce results aligned with macro-state interventions. Casert et al. [5] employ Transformers to model the transition of physical systems from one configuration another, in terms of their transition rates rather than reproducing individual agent behavior. Their method is tailored physical systems in continuous time, making inapplicable for interacting agents in the general case, since it requires explicitly enumerating allowed transitions between configurations.
In contrast to these approaches, our work is the first to jointly emulate individual and stochastic interacting agents. This is particularly important, since ABMs are inherently stochastic and rely on individual-level interactions to produce emergent aggregate outcomes. Moreover, since our surrogate differentiable by design, it paves the way for methods that estimate both individual-level parameters and state variables.
To achieve this goal, our framework relies on a novel combination of graph neural networks and diffusion models. Diffusion models [14] were first introduced in the context of image generation,
3

---

where they demonstrated impressive generation capabilities [7], and were then applied to other domains [18]. A number of works addressed graph data [21], for example in molecule modeling [15] and protein structure generation [2]. However, these works focus on the generation of graphs, while our architecture learns to generate random samples that are conditioned on information found on a graph. To the best of our knowledge, our work is the first application of this generative framework to individual behavior modeling in simulation systems, such as ABMs.
3 Methods
Denoting the set of agents by A, let each agent E A at discrete time be described by a state vector T which may include both continuous and categorical features. Given the ABM parameters θ the update rule of Z follows stochastic transition process Pe given by
z 4+1 Pa (Z Z (Z SEN(1))
(1)
where N:(() is the set of agents interacting with agent at time f. inducing a (time-varying) interaction graph G2 (A, EE) that assume to be known. This formulation focuses on individual agents, capturing not the dynamics of the entire system, but the evolution of each agent over time. In this way, it makes the two core ingredients of ABMs explicit: (i) relational structure via local neighhorhoods
N (ii) stochasticity in the choice of next states.
To effectively model these components in the same individual-level view, we leverage respectively (i) message-passing GNNs [10], which model the relationship between the evolution of an agent's state and the state of its neighbors on the graph; (ii) conditional diffusion models [44], generative architectures well-suited learning complex, multimodal distributions, allowing us to capture the intrinsic stochasticity of agent behavior.
Our proposed method, dubbed Graph Diffusion Network (GDN), combines these two components into a single architecture. Together, these components let us learn both how any agent state is affected by its neighbors on the graph, and the inherent randomness driving agent dynamics, yielding surrogate that can both emulate the original ABM and be differentiated.
Overview. To learn the distribution Pe, the training phase requires observations of different outcomes given the same starting conditions. To do so, in our framework, we use the original ABM to generate
data set a ramification of possible states, namely (z,"), [z!"") }JEN())
(i) (see Figure 1).
Our Graph Diffusion Network then approximates the stochastic kernel Pe by integrating a Message- Passing GNN with a Conditional Diffusion Model, of learnable parameters w and á respectively.
The GNN aggregates each agent's state z,(0) and its neighbors' states (Z(") JE via permutation-
invariant message and readout functions to produce interaction embedding (1) This embedding acts as a compact representation of the information coming from i's neighbors at t, affecting the
distribution of possible states of agent i at time + 1. As such, it is passed to the diffusion model: conditioned on Z, and (8) the diffusion model learns to transform a sample of Gaussian noise into a possible instance of the next state 1 By minimizing the standard denoising loss over all observed transitions, this hybrid architecture captures both the graph-structured interactions and the inherent stochasticity of agent dynamics. The trained model GDN is therefore able to generate,
given& a graph G2 of interacting agents and the state of each one Z a sequence of possible next states z,0) The consecutive application of GDN allows for reproducing the behavior of the original model. We now describe in detail each of these components.
dlesagg-paassig GNN. The GNN operates the provided interaction graph G. (A, E.), that we assume to be known or to be computable from Z. (e.g., in the Schelling model, the position of the agents determines who interacts with whom). For each agent i, the GNN aggregates its state z together with each neighbor's 's state Z(5) via a concatenated result through an MLP fa [10], that is gi (i) fa Z JEN (z() ) The
innaaiaann operator D and then feeds the (i)
(i) resulting vector gr captures how i's local neighborhood influences its next-state distribution. In

---

practice, for the operator b we used sum aggregation for Predator-Prey and mean aggregation for Schelling, since the latter's dynamics depend on the degree. Minimal experimental evaluation can guide practitioners toward the most suitable operator.
Conditional diffusion model. The diffusion model then learns the distribution over future states given this output from the graph and the state of given agent. Diffusion models do so by reversing a fixed Gaussian noising process [14]. The obtained denoising process, indexed by E (Tmax...... 0}, starts from a sample of Gaussian noise Tuu ~ (O, 1) and in a sequence of denoising diffusion steps transforms it into possible next state X RJ Z t+1 In this setting, we denote the general latent X, as Z (T). Each step of this process receives as input (i) the agent's is current state Z (ii) its
interaction embedding gi (6) and (iii) a sinusoidal positional embedding of r. These inputs are first transformed by MLPs to form the condition vector c!") Then, feed-forward network á is trained to
predict the noise residual cØ. i.e., the change to apply to the input 7 (r) to continue the denoising process.
Ramification data set. Given these two components, our framework uses the original ABM to produce a ramifications data set (see Figure 1). Such data set follows one main branch that specifies the evolution of the ABM, and multiple alternative stochastic evolutions of each time step from time t to time t + 1. This method makes it possible to expose the model to multiple stochastic successors from identical conditioning context, while avoiding exponential growth in the number of histories. Starting from an initial configuration Zu = (Z)) ..... z)") we recursively simulate R+++ 1 child configurations at each time step I, yielding (Z++1 1 or....o.....R R2. We designate the branch r = as the main branch [0 from which we extract the conditioning tuples (z(", {z(") )
for all agents i. The remaining R sibling branches at each t supply the target next states 1+11 ensuring that each context yields multiple outcomes.
Learning procedure. Our framework uses these data sets to train the Graph Diffusion Network. It minimizes the expected de- noising loss over the outcomes observed in the ramification data (see Algorithm 1). At each training iteration, it uniformly sam- ples time index and extracts the condi-
Algorithm 1: Training Procedure
1:
repeat
2:
Uniform O....,,--.,
3:
z
(z(")) JEN( [z
4:
( (Z)
+
JEN (z 2(0)))
tioning pair (2))) {z(²))
}) from the main
5:
Uniform(1 (1, .... Tmax
branch Z, 0). We compute the interaction
6:
Temb SinusoidaPosittooEmbbeeeegggg
7:
c (8) MLP(Z, MLP(g!") + MLP(Temb)
embedding gi via Equation (3), then draw
Uniform (1,
diffusion step T to form the condition vec-
tor (6) and uniformly select one of the R next-state realizations to obtain the target 11:
9:
Z 2 + Zr.i(r]
10:
~ 1)
Optimizers step over all i E A
z Finally, we minimize the denoising
12:
Vowlle - ea (2 (T),
( ct )|²
loss in Equation (5) by backpropagating through both the diffuusion model and the
13:
until convergence
GNN. More details about the architecture, and a discussion on computational costs and scalability can be found in Supplementary Section A.
4 Experiments
In this section, we present different experiments to assess and demonstrate our framework's 's ability to learn micro-level agent behaviors and faithfully reproduce emergent system-level dynamics. We evaluate our Graph Diffusion Network on two canonical agent-based models: the Schelling's 's segrega- tion model and the Predator-Prey ecosystem, presented in Section and detailed in Supplementary Section B. We test both its micro-level and macro-level fidelity. At the micro level, we measure how well the surrogate reproduces the conditional next-state distribution of each agent under identical con- text on out-of out-of-training ramification data set. At the macro level, we assess whether the surrogate, once trained on the first T Ttrain 10 timesteps, can accurately reproduce the subsequent Tust = 25 timesteps of aggregate summary statistics.
5

---

Because no existing method directly accepts graph-strictured agent states and outputs per-agent state distributions, there are no directly applicable baselines for our approach. Existing surrogates typically operate at the macro level, learning mappings from parameters to aggregate outcomes rather than reproducing full system dynamics. However, such approaches, including standard time- series models such as AR(1), fail to capture non-monotonic or cyclic patterns (e.g., predator-prey oscillations) and do not generalize beyond in-sample dynamics, as we show in Supplementary Section C.S. Therefore, we evaluate against two ablated variants. The first replaces the GNN embedding with flat concatenation of all agent states, removing relational structure. The second keeps the GNN but removes the diffusion component, predicting deterministic next states instead. Both ablations trained on the same ramified datasets and under identical protocols.
In the remainder of this section, we first describe the experimental setup, including dataset generation, model variants, and evaluation metrics. We then present a qualitative analysis of emergent patterns, followed by a comprehensive quantitative comparison. All implementation and reproducibility details are provided in the Supplementary Materials. Full code to reproduce our experiments is available at http://github. http://github..cm//faaacca httprrrroaeiiiaaaaaaaaa--iffftia
4.1 Experimental design
Ablation. Our core hypothesis is that both relational structure and stochastic modeling are crucial for accurate ABM surrogates. We consider therefore two possible ablations. In the first, we remove the message-passing GNN and replace the interaction graph with flat concatenation of all agents' state vectors this isolates the impact of neglecting agent interactions. The second drops the diffusion component entirely. yielding a purely deterministic model similar in spirit to prior GNN-based approaches for deterministic automata [13]. In this ablated version, we predict the next agent state by minimizing the mean squared error (MSE) with respect to the true next agent state. Details about the architecture can be found in Supplementary Section A. These ablations allow us to measure the improvement achieved by combining relational and stochastic modeling.
Agent-based models. We evaluate our approach on the two ABMs described in Section 2 as case studies. The first is the Schelling segregation model, in which n agents occupy cells on a two- dimensional grid. Each agent has fixed binary "color" and position the grid. At each timestep, an agent is considered happy if the proportion of its (up to eight) immediate neighbors sharing its color exceeds tolerance threshold ş; otherwise, it is unhappy and relocates empty cell; thus, the interaction graph G, links each agent to its neighbors at time t. We adopt the standard NetLogo implementation of this model [42]. The second is predator-prey ecosystem model, where agents belong to one of two species (predator or prey), inhabit grid cells, and cycle through life phases-Unborn, Alive, Pregnant, and Dead. At each timestep. an Alive agent may move to neighboring cell, reproduce (becoming Pregnant), or die, with probabilities specified by a parameter matrix V and conditioned on the local presence of predators or prey [39, 43). Pregnant agents revert to Alive after giving birth; Unborn agents become Alive if their parent is Pregnant; and Dead agents remain inactive. Here, G. links Alive neighboring agents, with specific rules for Pregnant and Unborn agents. See Suppiementary Section B for more details. In both ABMs, each agent's ' full state at time & comprises its position, type (color or species), and, for the predator-pre prey ABM, its life phase. Together, these two models span both simple relocation dynamics and richer birth-death interactions, providing diverse testbeds for our surrogate.
randomly selected
Micro evaluation metries. To quantify how faithfully our surrogate captures individual agent behavior, we compare its predicted conditional next-state distributions against the ABM's true stochastic transitions using the Earth Mover's Distance (EMD) [35]. We extend the ramification dataset beyond the training horizon and generate corresponding datasets for both the surrogate and the ablation models. The EMD is then computed as the mean value across timesteps and individual agents. In the Schelling ABM, the EMD compares the distribution of agent positions. This directly measures the surrogate's 's ability to relocate unhappy agents correctly, while keeping happy agents stationary. In the predator-prey model, we treat the agent's categorical life phase as the random variable and compute the EMD over its four-state distribution. This metric captures both deterministic transitions (e.g., Unborn -> Alive, Dead + Dead) and stochastic, interaction-drive transitions (e.g.. Alive + Dead, Alive + Pregnant).

---

E1
E3
15
30
30
15
30
Ablation Surrogate Original
Figure 2: Evolution of the position of black and red agents in the Schelling model, for three simulation runs, one for each of the considered tolerance thresholds & = 0.625. ε = 0.75, E3 = 0.875 (left, center, and right panel). We compare the ground truth (top row) with our surrogate (middle row) and with the best-performing ablation (according to sMAPE, bottom row). For each panel and model, we show three time steps: t (initial conditions, same for each column but kept for elarity), t 15, and = 30.
Macro evaluation metrics. Next, test whether agent-level predictions translate into faithful reproduction of emergent, system-level behavior. For each model, we track a summary statistic time: the number of happy agents in Schelling, and the number of active (i.e. Alive and Pregnant) agents in the predator-prey ecosystem. Reusing the same ramification branches as in state
training would offer little new information, since different stochastic branches from the tend to produce very similar macroscopic trajectories. Instead, we generate a fresh ensemble of main-branch simulations (100 independent runs) beyond the training horizon. We then compute the symmetric mean absolute percentage error (sMAPE) between the mean ground-truth trajectory and the mean currogate-predicte trajectory across this ensemble, providing quantitative measure the surrogate's 's ability to capture oscillations and steady-state behavior truly out-of-sampll
Experimental set-up. We consider three parameter combinations for the Schelling ABM, each producing distinct segregation outcomes, and four V combinations for the predator-pres ABM, reflecting different oscillatory patterns in the population dynamics. For each ABM and parameter setting, we simulate Tirain yielding the training ramification in Figure 1. For macro-evaluatiin main-branch simulations to calculate sMAPE. For microoevaluation, we generate an out-of-sample ramification dataset of T = 25 timesteps. We train both surrogate and ablations for 100 epochs using Adam with learning rate 10-5 for the diffusion model and Adam with learning rate 10-5 for the GNN, batch size equal to number of agents, and diffusion aaper-parameters Tmax 100 (more information Supplementary Section A).
main-branch steps with
500 stochastic branches per step, we run 100 independent
4.2 Results
To build intuition, we first qualitatively compare the surrogate and its ablated variant on their ability reproduce key emergent patterns of agent-based dynamics. We then consolidate these insights with a comprehensive quantitative evaluation using the macro- and micro-level metrics introduced in the previous section. We report selection of results in this Section; more Supplementary Section C.
Reproducing emergent segregation. Let first consider the Schelling ABM, under the config- urations & = 0.625. & 0.75, E3 0.875. Figure illustrates how the ground-truth ABM (top row) progresses from randomized initialization structured, segregated communities for the first two configurations {1, &z, while it remains unsegregated for Es. At the first two tolerance levels. in fact, the agents gradually self-organize into distinct clusters, with segregated communities clearly emerging around / = 20 (see Supplementary Section C). The middle row represents the evolution of the system according to our surrogate model: we initialize the system with the same starting condition, and then we iterate giving the current state Z, to model, and using one sample of

---

1 Oscillations predators and preys
Number of ageents
Ground
Sumrogate model (Ours)
Diffusion-only
GNN-only
200
Preys Predators
10
20
30
10
20
30
10
20
30
50
Time
Time
Time
Time
V4: Oscillations predators, monotone preys
400
Ground truth
Surregate model (Durs)
Diffusion-oni
GNN only
Prays
Predators
Number
200
200
Time 30
Time
Time
Figure 3: Forecasting macro-level summary statisties (here, the number of alive preys and predators over time), starting from the last condition seen in training, for 100 independent simulation runs, under configuration v (oscillations for both predators and preys, top) and V (oscillations only for predators, bottom). From left to right: original ABM simulations. surrogate, diffusion-only ablation, GNN-only ablation. The dashed vertical line indicates the end of the training phase for the surrogate and ablation models.
the generated output as the next state Ze+1. We observe that our surrogate exhibits a qualitatively similar pattern of cluster formation over time, distinct for each of the three configurations, Instead, the ablation models fail to meaningfully relocate agents. The best- rest-performing one, according to sMAPE, the diffusion-only model, largely maintains random configuration.
Reproducing emergent oscillations in predator-prey dynamics. Next, we consider the Predator- Prey ecological model. Figure 3 overlays 100 trajectories of prey and predator populations starting from the same state at the end of training, comparing the stochastic trajectories from the ground- truth model with those obtained by surrogate and by the ablation. For both configurations, the surrogate and the ablation models are trained only with the initial time steps (up to the dashed line in the plots). Under the parameter set 1, the ground-truth ABM (top-left plot) exhibits classical Lotka-Volterra oscillations: rise in prey growth drives delayed increase in predators, which then triggers prey decline and a subsequent predator decline. Under 44, instead, only predators show a rise and decay, while preys only decline (bottom-left plot). We observe that the surrogate (second column) accurately captures both the phase lag and amplitude of these oscillations, while the diffusion-only ablation (third column) collapses to near monotonic trends. The GNN-only ablation (fourth column), besides following completely deterministic dynamic, completely diverges from the ground truth. The number of preys under Vi, for instance, quickly reaches a plateau at value of 750 (not shown), almost twice the real one. We perform the same analysis for alternative parameterizations 2/ V/ (included in Supplementary Section C) that show different types of dynamics, as the population of predators and/or preys may exhibit monotonic extinctions. In all cases, the surrogate faithfully reproduces monotonic declines single-peak dynamics and both ablations fail. We also observe (figures in Supplementary Section C) that the surrogate recreates the rich spatial patterns of predator-prey clusters, also seen in similar settings in evolutionary game theory [27].
Quantitative results. Now we present the results of quantitative analysis, systematizing the previous comparisons. Here, each comparison with the ground truth is quantified using of the metrics presented in the previous subsection, i.c. Earth Mover's Distance (EMD) for the micro- level comparisons, and sMAPE for the macro-level ones. Figure summarizes the results of our experiments: the left panel shows the microscopic evaluation of both our surrogate model and the ablated variant, while the right panel presents the macroscopic evaluation results.
For the Schelling model, we observe that, on the micro level, the surrogate's mean EMD is lower than the diffusion-only ablation's mean EMD in all cases. The differences between the surrogate and the diffusion-only ablation are less pronounced the thresholds & and E3. At & (few unhappy agents) behavior is almost entirely deterministic and agents rarcly move, while at & (almost all
8

---

Surrogate (ours)
Diffusion-one only
GNN only
Schelling
c
Low segregation
High segregation
segregation
10
0.0
0.5
1.0
1.5
2.0
2.5
Predator-
Osciliatoons Predators
41
Osciliations Preys
Predators Monolone Preys
Monotone Predators
Osciliations Preys
Oscillations Predators
Monotome Prays
0.02
0.04
0.00
0.25
0.75
1.00
2
1.50
EMD (Micro)
sMAPE (Macro)
Figure 4: Errors obtained by the proposed approach (Surrogate) and by the naive baselines (Diffussio- only and GNN-only ablation models) in four different tasks. In the first column. error is measured as the EMD between the true and predicted distribution of individual (micro-level) behavior, i.c. predicting the next state of each agent from the previous one. In the second column, error measured as the difference (sMAPE) in system-level quantities, i.e. comparing the true values of the number of agents with given state with the one predicted by our model when trained on a fraction of the initial time steps (as in Figure 3). In the first row, we test three configurations of the Schelling model; in the second row, we compare four configurations of the Predator- Predator-Prey model.
agents unhappy) behavior is uniformly random, even flat, "always-move" "never-move" rule yields near-optimal predictions these two cases. In contrast, at the intermediate threshold &z, where roughly half the agents unhappy, the difference between the surrogate's 's and the diffusion-only ablation's EMD is more pronounced. A similar pattern is observed in the macroscopic evaluation. The surrogate's 's sMAPE remains below 0.2, whereas the diffusion-only ablation fails to distinguish happy from unhappy cases, resulting in large macro-level errors. In the GNN-only ablated model, the micro-level EMD increases proportionally with the level of stochasticity introduced by the parameters, with EMD values similar to the surrogate for & (the most deterministic) and much higher & (the most stochastic). At the macro-level. we observe large errors, with sMAPE always higher than the surrogate. These gaps confirm that only the full model, with explicit graph-based interaction modeling and stochasticity, can learn the conditional relocation rule critical balanced regimes.
For the Predator Prey model, regarding micro-level behavior, our surrogate achieves a low EMD from the ground truth on average, and it consistently outperforms both ablation models. These results confirm that our model is able to faithfully reproduce the complex dynamics of this ABM even at the individual agent-level (thus explaining Figure 3). The most successful case is V2 where surrogate exhibits a near-zero difference from the ground truth. In fact, in this configuration, most agents follow deterministic update rules (c.g., dead - dead), which are perfectly recovered by our model, but not by the diffusion-only ablation - which also obtains worse results on stochastic rules, as shown in Supplementary Section C.4. Instead, the GNN-only ablation performs well only in those deterministic cases, but fails in all the others. At the macro level as well, the surrogate consistently outperforms the ablation, generally achieving low error. The best result is obtained with V1 the most complex dynamics, where the surrogate achieves an average sMAPE of approximately 0.08. This configuration produces two distinet population peaks, and the surrogate faithfully reproduces both their timing and amplitude (Figure 3). The worst result is obtained with V2 as this configuration is almost monotonic and dominated by long, near-zero tails that are noisy very small scales, making them difficult for any model to reproduce.

---

Discussion
We introduced Graph Diffusion Networks, differentiable surrogate for agent-based models that combines graph neural networks to model agent-level interactions with diffusion models to capture stochasticity. Our experiments the Schelling segregation model and Predator- ecosystem show that this approach not only accurately reproduces individual-level transition distributions, but also faithfully captures emergent, system-level dynamies beyond the training horizon.
Limitations. Our approach is limited by our assumptions about the characteristics of the ABM to emulate. First, the interaction graph is assumed to be fully known. Future work might remove this limitation by estimating such graph directly from available data. However, the estimation of a latent interaction graph is a follow-up challenge, for which our GNN-based approach represents a necessary first step. Second, highly sophisticated ABMs may include features not addressed in our framework such as all-to-all interactions, multiple rounds of decision decision-making, or sequential stochastic events within single time step. Capturing these dynamics may require extending our architecture to incorporate sequential or hierarchical components. While our method may not yet fully generalize to such settings, our findings demonstrate that building surrogates capable of replicating -edividual-level behavior is both feasible and effective, laying the groundwork for broader applications.
Future work. Building on this foundation, the differentiability of our surrogate opens up range of powerful applications. It enables the use of gradient-basee methods for any optimization task, such policy optimization [1]. It allows for efficient calibration of macro parameters by treating key parameters as additional inputs to the neural network. Most importantly, our approach naturally allows for the estimation of micro (i.e., agent) level variables a challenge for ABMs, that often requires the ad hoc development of handcrafted probabilistic models [24, 25]. In fact, our model
already contains such parameters expressed as agents' individual states (Z("), ), something typically not available in ABM surrogates [11]. Morcover, our method can in principle be applicd directly to real- world datasets whenever sufficient observations of comparable agent-contex. pairs and transitions available. In doing so, our framework helps make ABMs more data-driven and empirically grounded, with promising applications in several scientific domains, such economics [29], epidemiology [12), sustainability [19], urban science [4], and ecology [38].
Acknowledgments
The authors wish to thank Daniele Grattarola and Federico Cinus for insightful early discussions that supported the initial development of this work. We also thank Alberto Novati for his contribution to the early draft of the code for the original ABM of the predator-prey system.
References
[1] Akash Agrawal, Joel Dyer, Aldo Glielmo, and Michael Wooldridge. Robust policy design in agent-based simulators using adversarial reinforcement learning. In The First MARW: Multi- Agent Al in the Real World Workshop at AAAI 2025, 2025.
[2] Namrata Anand and Tudor Achim. Protein structure and sequence generation with equivariant denoising diffusion probabilistic models. NeurIPS 2022 Workshop on Machine Learning in Structural Biology, December 2022.
[3] Gaurav Arya, Moritz Schauer, Frank Schäfer, and Christopher Rackauckas. Automatic differ- entiation of programs with discrete randomness. Advances in Neural Information Processing Systems, 35:104355110444 2022.
[4] Mark Birkin, Patrick Ballantyne, Seth Bullock, Alison Heppenstall, Heeseo Kwon, Nick Malleson, Jing Yao, and Anna Zanchetta. Digital twins and ai for healthy and sustainable cities. Computers, Environment and Urban Systems, 120: 102305, 2025.
[5] Corneel Casert, Isaac Tamblyn, and Stephen Whitelam. Learning stochastic dynamics and predicting emergent behavior using transformers. Nature Communcationns, 15(1):1875, 2024.
[6] Claudio Castellano, Santo Fortunato, and Vittorio Loreto. Statistical physics of social dynamics. Reviews of modern physics, 81(2):591-646, 2009.
10

---

[7] Florinel-Alin Croitoru, Vlad Hondru, Radu Tudor Ionescu, and Mubarak Shah. Diffusion Models in Vision: A Survey IEEE Transactions on Pattern Analysis & Machine Intelligence. 45(00))108500-0865 September 2023.
[8] Joel Dyer, Nicholas Bishop, Yorgos Felekis, Fabio Massimo Zemnaro, Anisoara Calinescu, Theodoros Damoulas, and Michael Wooldridge. Interventionally consistent surrogates for complex simulation models. Advances in Neural Information Processing Systems, 37, 2024.
[9] Marian Farah, Paul Birrell, Stefano Conti, and Daniela De Angelis. Bayesian emulation and calibration of dynamic epidemic model for a/hInl influenza. Journal of the American Statistical Association, (19550))1398 - 2014.
[10] Justin Gilmer, Samuel S. Schoenholz, Patrick Riley, Oriol Vinyals, and George E. Dahl. Neural message passing for quantum chemistry. In Proceedings of the 34th International Conference on Machine Learning, volume 70, pages 1263-1272. PMLR, 2017.
[11] Aldo Glielmo, Marco Favorito, Debmallya Chanda, and Domenico Delli Gatti. Reinforcement learning for combining search methods in the calibration of economic abms. In Proceedings of the Fourth ACM International Conference on Al in Finance, pages 305-313, 2023.
[12] Nicolò Gozzi, Matteo Chinazzi, Jessica Davis, Corrado Gioannini, Luca Rossi, Marco Ajelli, Nicola Perra, and Alessandro Vespignani. Epydemix: An open-source python package for epidemic modeling with integrated approximate bayesian calibration. medRxiv, 2025.
[13] Daniele Grattarola, Lorenzo Livi, and Cesare Alippi. Learning graph cellular automata. Ad- vances in Neural Information Processing Systems, 34:20983-20998 2021.
[14] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. In Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS '20, Red Hook, NY, USA, 2020. Curran Associates Inc.
[15] Emiel Hoogeboom, Victor Garcia Satorras, Clément Vignac, and Max Welling. Equivariant diffusion for molecule generation in 3D. In Proceedings of the 39th International Conference on Machine Learning, volume 162, pages 8867-8887. PMLR, 2022.
[16] Ruhollah Jamali, Wannes Vermeiren, and Sanja Lazarova-Molnar Data-driven agent- based modeling: Experimenting with the schelling's 's model. In Procedia Computer Science, volume 238, pages 298-305, 2024.
[17] Mare C Kennedy and Anthony O' D'Haga.. Bayesian calibration of computer models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(3):425 464, 2001.
[18] Akim Kotelnikov, Dmitry Baranchuk, Ivan Rubachev, and Artem Babenko. Tabddpm: modelling tabular data with diffusion models. In Proceedings of the 40th International Conference on Machine Learning, ICML'23. JMLR.org. 2023.
[19] Francesco Lamperti, Giovanni Dosi, and Andrea Roventini. A complex system perspective on the economics of climate change, boundless risk, and rapid decarbonization. Technical report, LEM Working Paper Series, 2025.
[20] Francesco Lamperti, Andrea Roventini, and Amir Sani. Agent-based model calibration using machine learning surrogates, Journal of Economic Dynamics and Control, 90:366-389, 2018.
[21] Chengyi Liu, Wenqi Fan, Yunqing Liu, Jiatong Li, Hang Li, Hui Liu, Jiliang Tang, and Qing Li. Generative diffusion models on graphs: Methods and applications. In Proceedings of the T'hirtySSccond International Joint Conference on Artificial Intelligence, page 6702, 2023.
[22] Thomas Lux. Estimation of based models using sequential monte carlo methods. Journal of Economic Dynamics and Control, 91:391-400 2018
[23] A. J. McLane, C. Semeniuk, G. J. McDermid, and D. J. Marceau. The role of agent-based models in wildlife ecology and management. Ecological Modelling, 222(8):5544, 201
[24] Corrado Monti, Gianmarco De Francisei Morales, and Francesco Bonchi. Learning Opinion Dynamics from Social Traces. In ACM, KDD, pages 764-773, 2020.
11

---

[25] Corrado Monti, Marco Pangallo, Gianmarco De Francisci Morales, and Francesco Bonchi. On learning agent-based models from data. Scientific Reports, 13(1):9268, 2023.
[26] Kilian J. Murphy, Simone Ciuti, and Adam Kane. An introduction to agent-based models as an accessible surrogate to field-based research and teaching. Ecology and Evolution, (0(22):12482 12498, 2020.
[27] Martin A Nowak and Robert M May. Evolutionary games and spatial chaos. Nature, (5966398):826-822 1992.
[28] Yannick Oswald, Keiran Suchak, and Nick Malleson. Agent- models of the united states wealth distribution with ensemble kalman filter. Journal of Economic Behavior & Organization, 229; 106820, 2025.
[29] Marco Pangallo and R Maria del Cho-Chhnona. Data-driven economic agent-based models. In The economy as an evolving complex system IV. SFI Press, Santa Fe, N.M., 2025.
[30] Anton Pichler, Marco Pangallo, R Maria del Rio-Chanona, François Lafond, and J Doyne Farmer. Forecasting the propagation of pandemic shocks with dynamic input-output model. Journal of Economic Dynamics and Control, 144:104527, 2022.
[31] Donovan Platt. A comparison of economic agent-based model calibration methods. Journal of Economic Dynamics and Control, 11:100385,, 2020.
[32] Sebastian Poledna, Michael Gregor Miess, Cars Hommes, and Katrin Rabitsch. Economic forecasting with an agent- based model. European Economic Review, 151:104306, 2023.
[33] A. Quera-Bofarull G. Dyer, A. Calinescu, J. D. Farmer, and M. Wooldridge, Blackbirds: Black-box inference for differentiable simulators. Journal of Open Source Sofrware, 8(89), 2023.
[34] Steven F Railsback and Volker Grimm. Agent-based and sibividuul--basee modeling: a practical introduction. Princeton University Press, 2019.
[35] Yossi Rubner, Carlo Tomasi, and Leonidas Guibas. The earth mover's distance as a metric for image retrieval. International Journal of Computer Vision, 40(2):99-121, 2000.
[36] Thomas C Schelling. Dynamic models of segregation. Journal of mathematical sociology, 1(2):143-18, 1971.
[37] Squazzoni. The impact of agent-based models in the social sciences after 15 years of incursion. History of Economic Ideas, (V11012):100011377 2010.
[38] Michiel Stock, Olivier Pieters, Tom De Swaef, and Francis Wyffels. Plant science in the age of simulation intelligence. Frontiers in Plant Science, 14:1299208 2024.
[39] Daniel Tang and Nick Malleson. Data assimilation with agent-based models using markov chain sampling. arXiv preprint arXiv:2205.016/66 2022.
[40] Guus Ten Broeke, George A.K. Van Voorn, and Arie Ligtenberg. The use of surrogate models to analyse agent-based models. Journal of Artificial Societies and Social Simulation, 24(2), 2021.
[41] Ruta Ubarevičiené, Maarten van Ham, and Tiit Tammaru. Fifty years after the schelling's 's models of segregation: Bibliometric analysis of the legacy of schelling and the future directions of segregation research. Cities, 147:104838, 2024.
[42] Uri Wilensky. Netlogo, 1999. Center for Connected Learning and Computer- Based Modeling, Northwestem University, Evanston, IL.
[43] Uri Wilensky and William Rand. An introduction to agent-bused modeling: modeling natural, social, and engineered complex systems with NetLogo. Mit Press, 2015.
[44] Zheyuan Zhan, Defang Chen, Jian-Ping Mei, Zhenghe Zhao, Jiawei Chen, Chun Chen, Siwei Lyu, and Can Wang. Conditional image synthesis with diffusion models: survey. Transactions on Machine Learning Research, 2025.
12

---

Learning Individual Behavior in Agent-Based Models with Graph Diffusion Networks Supplementary Material
A Neural models and training details
In this section, we provide a detailed overview of the core components of the Graph Diffusion Network (GDN) and its methodology. We begin by introducing the diffusion process (A.1), which defines the diffusion process to be reversed to generate future agent states starting from a sample of Gaussian noise. Next, we detail the Graph Diffusion Network architecture (A.2) and its components. We then discuss the loss and optimization strategy (A.3), covering the training objectives and gradient flow between the diffusion model and graph components. Following this, we outline the generation algorithm (A.4), where the iterative denoising process generates future agent states. Finally, we provide a discussion on computational costs and scalability our methodology (A.5).
A.1 Diffusion process
Our diffusion model is designed to generate future agent states 2 by reversing a known Gaus- sian noising process (i.e. the forward process) through of latents 2 2(41(T) indexed by {Tmaxs ..., 0}. The forward process is fixed Markov chain that gradually adds Gaussian noise to the input z(() according to previously defined variance schedule B,. Each latent given by:
(1) (41 (T) Va (8) (+1 -
e~,
where 0r = 3r and a, = II Og. For our model, we selected a cosine variance schedule:
8, Sstart (Bend Batars) )1 cos(
x))
Tmaz
with Betart = 10-4 and Bend = 0.02. This choice ensures that 3, increases more gradually at the beginning of the forward process, retaining more of the original input information, and at the end of the forward process. We note that, in preliminary trials, it showed to be more stable in our scope with small input dimensions compared to the cosine variance schedule proposed by [48] in the context of image generation.
A.2
Graph Diffusion Network architecture
The primary input of the conditional diffusion model inside the Graph Diffusson Network is the latent 7 (i)
2441 (T). noised version of Z given by equation (2). In general, not all variables contained in
z 1+11 are time-dependent, and some remain stationary through time (e.g. color in Schelling and kind in Predator-Prey, see Supplementary subsections B. B.2). We only include the time-dependent features (or dynamical features) in the input of the diffusion model, they are the ones that evolve over time
and need to be predicted. The output the diffusion model is the denoising step e ( of introduced in Section 3, which has the same size as the input. Thus, the neural network follows symmetrical structure with hidden layers of increasing width the first half, and decreasing width
in the second. The condition vector c is applied to the hidden layers of the neural network by applying an activation function, performing linear operation to match the width of layer and summing element-wise with the hidden layer. To increase stability, hidden layers are first normalized and there is a residual connection after conditioning has been applied. All details the architecture of the conditional diffusion model are reported in Table 1.
13

---

The Message Passing GNN takes in input the entire agent state Z,") node features. The messages correspond to the node features and are aggregated by an aggregation function such as sum or mean value. The choice of the aggregation function depends the ABM to be reproduced. In general, sum a suitable choice, as the MLP fw will capture the behavior rules of the agents. However, for ABMs where the behavior of agents is influenced by the node degree, in the of Schelling, mean is more appropriate choice. All details of the architecture of the Message Passing GNN are reported in Table 1.
To make the network more stable, all features are scaled. particular, agent states Z can contain both numerical and categorical features. Numerical features are scaled to the interval [-1, 1]. In our experiments, we scaled numerical features with standard scaler. After generation, they are scaled back to their original domain and, for integer numerical features, binning function is applied afterwards. Categorical features are one-hot encoded.
Table I: Architecture and training details of the Graph Diffusion Network Details
Component
Conditional diffusion model
Input dimension
aaannicalfffaturee __dim
Hidden layers
[128, 256, 1024, 1024, 256, 128] dynamical_features.
Output dimension
Activation function
LeakyReLUU (slope 0.1)
MLP time embedding
Linear(256) > Act + Linear(256)
MLP current state
Linear(256) - Act Linear + Act Linear(250 (256)
MLP graph embedding Condition block in hidden layers
Linear(256) - Act Linear(256) + Act + Linear(256)
LayerNorm
Linear(dim _ottiimm_mtt____oo
(um((Lin(Ac(cconditioo))) + Linear(dim_out) - Residual
connection
Weights initialization
Xavier uniform
Optimizer
Adam 10-5
Learning rate
Message Passing GNN
Input dimension
X agent_state_da
Hidden layers
[32, 64, 128]
Output dimension
256
Aggregation function
sum mean
Activation function
LeakyRel.U (slope 0.1)
Message passing
Message [z)) Ij
Weights initialization
Kaiming uniform
Optimizer
Adam
Learning rate
X hateniiin_rtt__diiiiiiiee
Other details
Twan Batch size Number of epochs
100
Number of agents the system
100
GNN-only ablated model The only architectural change of the GNN-only ablated model is in the MLP following the aggregation pass of GNN: output dimension set to the agent-state dimension rather than the embedding dimension. While the original GNN uses hidden layers of size [32, 64, 128], the ablated model adopts symmetrical structure: [32, 64, 128, 128, 64, 32]. We trained for 100 epochs with the Adam optimizer (learning rate 2 10 10-5) and mini-batch size of 16.
14

---

A.3 Loss and optimization
The learning objectives of the Graph Diffusion Network are the noise residuals 6 (T), (t of
(0)
the denoising diffusion process in Tmoz: 0 used to generate ( 1++1 given (Z,") (2
(1)) (EN(())
The generative process is conditioned by the condition vector of which is learned by the network and is given by:
ML.((Z))) (T. MM.((Temb)
(4)
where Temb is the sinusoidal positional embedding of and gi the embedding produced by the Message Passing GNN of parameters w. More details on the three MLPs that form c, are given
in Table 1. The loss function for the denoising diffusion steps ca c2 (T). cf is calculated as the expected value over agents in the system E A, all ABM timesteps I, all and ~ N/ I):
((,w) E,,,, [|| | e ( (Z of
(5)
The parameters to optimize are the parameters of the diffusion model & and the parameters of the GNN w. At each training step, the loss in equation (5) is calculated over the batch made of all agents i E A and gradients are backward propagated. First, the optimizer for Φ is applied, and then the optimizer for w. Thus, the GNN is trained by inheriting the gradients from the loss of the conditional diffusion model, through the learned condition representation c,
A.4 Generation
The generation of 2 starts from the last latent of the denoising diffusion process, which is a sample of Gaussian noise ( (Tmax) ~ (N, I). The Message Passing GNN takes in input the current state JJEN
z (i) and the states its neighbors (Z(3)) { and forms embedding g:. Then, iteratively over Tmaz, ..... 1, the conditional diffusion model takes in input the sinusoidal positional embedding Temb: the current agent state z and the embedding vector gf and forms the condition vector c,
Lastly, the previous latent 2+11 (T 1) is calculated given 2 and the output of the Graph Diffusion Network c (2 ( (T), of (0) as in lines 7-8 Algorithm 1.
(
Algorithm 1 Generation
Z (z
(jEN(0
Data
2: gt (Z,
JEN(O (z(") z(())
3: [ (Tmaw) ~ (0.1)
4:
for Tians .... 1 do
5:
Temb = Sinussidapssttoonbbeddnggg
(1) MLP(Z + (LP(g.") + MLP (Temb))
6:
c
7:
2 C~, if else z =
8:
(T- - (Z,2, (2,
1-ar co (Z (Z1 o,z
VI-ar
9:
end for
10: return
2 (0) (0) 1411
We set I- 1-ū, Br. This choice of a, is optimal for deterministically set points [14], which the case for some update rules in ABMs (e.g. happy agents in Schelling and dead agents in (redator-Prey), Alternatively, one can also choose a, v3,. which is more optimal for normally distributed points.
A.5 Computational costs and scalability.
In the following subsection, we provide additional details and discussion on the computational costs and scalability of our proposed methodology. In particular, details the hardware employed to run
15

---

our experiments, performance details with respect to the number of ramifications provided during training, as well as discussion over the scalability of Graph Diffusion Networks to larger systems of agents, and on the interplay between the number of agents and the number of ramifications.
Runtime characteristics. All experiments were run in cloud-based server with 15 vCores, 180 GB of RAM, and an NVIDIA A100 80GB PCIe GPU. Execution times depend on the size of the datasets and GPU occupancy by other processes. In our experiments, where the training datasets were made of R 500 ramifications over T = 10 timesteps with 2048 agents in the system for Predator-Prey and 1950 agents for Schelling, the training time over 100 epochs typically lasted around I hour (around 37 seconds per epoch). Generating a simulation of 25 timesteps for the entire system of agents takes roughly 7.5 seconds for both Predator- Prey and Schelling, around 0.3 seconds per timestep.
Scalability with respect to number of ramifications. During training, each training step processes the entire system of agents for one timestep and one ramification. Therefore, the runtime grows linearly with respect to the number of ramifications or timesteps. We present results from 10 experiment runs trained on one of the Predator-Prey Prey datasets of parameter 1/1 to compare training time, micro-level metric (EMD), and macro-level metric (sMAPE) as the number of ramifications provided during training increases in Figure 5 and Table 2. As expected, training time increases linearly with the number of ramifications included in the dataset used to train the model. Furthermore, micro-level metrics decrease the number of ramifications increases, showcasing the benefit of providing a higher number of stochastic outcomes to the model during training to better fit the stochastic rules that govern the original ABM. Macro-level metrics are subject higher variance and do not show a trend as clear as with the micro-level metrics, still the best results are reached with the highest number of ramifications.
Scalability with respect to number of agents. The experiments covered in Section 4.2 are per- formed on mid-sized ABMs, with the total agent count ranging in the thousands (for Schelling there are 1950 agents in total and 2048 for Fredator-Prey) Increasing the number of agents in the system (for example by increasing the grid size and keeping the density of agents constant) increases the training time sub-linees libearinrrl. In fact, with mid-sized systems, we can train the network with batches that correspond to the whole set of agents. Thus, the number of training iterations per epoch does not change, but rather the size of the batch in input to the network. It should be noted that increasing the number of agents also increases the inference time (time required by the trained surrogate to produce a simulation for the whole system). We present results from 5 experiment runs trained on one of the Predator-Pei Prey datasets of parameter v. 1. to compare training time, inference time, micro-level metric (EMD) and macro-level metric (sMAPE) as the total number of agents in the system increases by increasing grid size and keeping agent density constant in Figure and Table 3. From our results, it can be noted that micro-level metrics slightly decrease as agent count increases, as well as macro-level metrics reach lower values in some of the experiments with higher agent count. We argue that increasing the number of agents in the system naturally increases the number of agent transitions available in the dataset, expanding the statistical coverage provided to the model to learn, and yielding slightly better results.
Scalability interplay between number of agents and ramifications. Generating hundreds of futures per state can be costly for many-agent long-horizon systems, but larger systems naturally provide more independent samples of similar conditions; thus, fewer ramifications can suffice for comparable statistical coverage. To evaluate the trade-off between data generation cost and predictive quality, we performed 5 experiment runs and trained our model on one Predator-Prey dataset (parameter set V )) with increasing agent counts and decreasing ramifications, keeping the training time per epoch approximately constant, and the total number of agent transitions roughly similar. All other hyperparameters match the main experiments. We observe from Table 4 that the errors (both sMAPE and EMD) remain close to the original results even when the number of agents grows substantially, keeping training time low thanks a lower number of ramifications needed.
16

---

Table Effect of ramifications R on training time and performance.
Ramifications
Time/epoch
sMAPE
EMD
50
3.71 + 0.02
0.076 = 0.019
0.0161 + 0.0003
100
7.44 + 0.02
0.071 = 0.015
0.0156 + 0.0004
250
18.61 + 0.06
0.083 = 0.039
0.0115 0.011
500
37.17 + 0.12s
0.062 = 0.030
0.0088 + 0.0007
Table Scaling with grid size, number agents with fixed number of ramifications (R 500).
*: Mean time produce simulation of the whole system for t = 25 timesteps.
Grid size
Agents
Time/epoch
Inference time*
sMAPE
EMD
64
8192
108.58 +: 0.30
35.32 + 2.03
0.051 + 0.026
0.0082 + 0.0003
48
4608
63.78 + 0.32 s
16.45 0.72 s
0.053 0.026
0.0083 0.0005
32
2048
37.20 + 0.12 s
7.44 + 0.01
0.065 + 0.033
0.0092 + 0.0007
Table 4: Scaling with grid size, number of agents, and ramifications while keeping training time and transitions roughly constant.
Grid size
Agents
Ramifications
Time/epoch
s.APE
EMD
8192
180
38.88
0.18
0.066 0.027
0.0095
0.0005
48
4608
300
38.41 + 0.25
0.073
0.024
0.0092 + 0.0005
2048
500
37.20
0.12
0.065 = 0.033
0.0092 + 0.0007
17

---

Training
Micro
40
0.0175
0.125
stepoch
30
0.0150
SMAPE
0.100
20
0.0125
0.075
10
0.0100
0.050
0.0075
0.025
100
200
300
400
500
100
200
200
400
500
100
200
300
400
500
ramifications
ramifications
ramifications
Figure 5: Training time, micro and macro metrics with respect to the number of ramifications provided during training for experiments on one dataset of Predator- Prey with parameter set 1 Points indicate the mean value and error bars standard deviation over the 10 experiments.
Training time
Micro
Macr
0.010
s/epoch
90
0.009
SMAPE
60
0.008
30
0.007
4000
6000
8000
2000
N. agentss
N. agents
N. agants
Figure 6: Training time, micro and macro metrics with respect to the total number of agents in the system for experiments on one dataset of Predator- with parameter set V Points indicate the mean value and error bars standard deviation over the 5 experiments.
Training
Micro
Macro
40
0.10
slepoch
39
38
EMD
0.009
SMAPE
0.06
0.08
37
0.04
36
0.008
0.02
2000
4000
E000
8000
2000
4000
6000
8000
2000
4000
6000
8000
N. agents
N. agents
N. agants
Figure 7: Training time, micro and macro metrics with increasing the total number of agents in the system and decreasing the number of ramifications to yield similar training times. Results are shown for 5 experiments on one dataset of Predator-Pre Prey with parameter set Vr Points indicate the mean value and error bars standard deviation over the 5 experiments.
18

---

B
ABM Case Studies
B.1
Schelling model
The Schelling model of segregation is a classic example to showcase the emergence property of ABMs. Agents i A placed in 2-dimensional grid I. x L. Their state is given by their color (a binary variable such as black and white) and their position on the grid.
z = (C(+)
c E (C1, C2).
E 0,
Each agent has a fixed color CO), which remains constant over time, while their position on the grid may change. The set of agents that interact with agent i, denoted E Wi),, includes those in the eight adjacent cells (Moore neighborhood) (x) 3/
The ABM mostly depends on a parameter & E [0, 1J. representing the intolerance of the agents. If the fraction of neighbors j (i)) that share the same color as agent is greater than or equal to &. agent i is considered happy and remains in its current position:
(4 Bi41) = (z, ))
Conversely, if the fraction is less than agent is considered unhappy and moves to randomly chosen empty cell on the grid. Thus, the update rule is deterministic when agents are happy, and stochastic when they are unhappy.
Algorithm 2 presents the pseudo-code of the ABM. Of particular interest are lines 15-25, which describe how agents relocate by searching for an empty cell on the grid. It is clear that this search process is not simple draw from a probability distribution, as in the framework by Arya et [3]. but a much more complex trial and error process. The pseudocode makes it clear that agents are more likely to relocate to nearby positions rather than distant ones, due to the way direction and distance are sampled. This spatial bias will be quantitatively confirmed in Figure 16.
B.2 Predator-Prey model
The Predator-Prey ABM is simulation model that captures the dynamics of interacting populations over time. We use a slightly adapted version of the model introduced in Ref. [39] (Algorithm presents the pseudo-code of our ABM). Agents E A occupy two-dimensional grid of size I. x T.. Each agent's 's state at time t is given by its kind (either Prey or Predator), its life phase (Unborn, Alive, Pregnant, or Dead), and its position on the grid:
Z (K(())
(3) (1) Jt 24
k( (Prey, Predator).
(1) 1 (Unborn, Alive, Pregnant, Dead), I ,r, E L
The agent's kind (G) is fixed over time, while the phase and position evolve. The set of interacting agents j E A/()) consists of those located the four cardinally adjacent cells to (xx(0) yE (()
(Von Neumann neighborhood).
In addition, each Unborn agent is assigned parent agent of the same kind, provided is either Alive or Pregnant. This parent-child relationship governs the birth mechanism.
The dynamics more complex than in Schelling's 's model of segregation. An Alive agent can transition to one three states: remain Alive, become Pregnant, or become Dead. These transitions stochastically determined. If the agent remains Alive, it moves random to cardinally adjacent cell. If it becomes Pregnant, it remains in place. it becomes Dead, it loses its position on the grid. Dead agent remains dead and off-grid. A Pregnant agent returns being Alive in the same cell. An Unborn agent becomes Alive only if its assigned parent currently Pregnant; otherwise, remains Unborn.
This gives rise to a set of deterministic update rules:
Alive Unborn
if parent is Pregnant otherwise
Dead Dead, Pregnant Alive, Unbor
19

---

Random direction, distance up dmax Algorithm 2 Schelling Model of Segregation
Require: Agent set A, grid size L, tolerance ç, max steps T, max distance dmax, max trials
1: Initialize: For cach E A, sample ( Uni,orm(Ci,Ca,, (x 36 " ~ UniformGrid(L)) and set - = (C(()) (6) 3o ) 2: for t 0,
do
3: unhappy +
Reset unhappy list
4:
for all A do
5:
N1 < {j A (xt) ) E Moore(a (0) .2i y
Moore neighborhood
6:
(jENs (b) H. 1
Same-color neighbors, all neighbors
s/n n>0
7:
P Similarity ratio
n=0
8:
If < then
Agent is unhappy
9:
unhappy + unhappy (tt)
10:
end if
11:
end for
12:
if unhappy then break
Convergence
14:
end
15:
for all i unhappy do
16:
(x M + (a (1) ( for k ..... K do
(8)
17:
18:
Uniform(0, 2x), d~ Uniform(0, dmax)
Ax [dcos 0J. Ay <- [dsin
Convert to grid movement
x= +- (xt + Ar) mod (L), g² (v, + Ay) mod (L)
Wrap around border
if 3j + i: (xx .14 "e","") then
Check if cell empty
(a Me41) (x", ,8"
Move to new location
23:
break
D Stop searching once valid position is found
24:
end if end for end for
25:
26:
for all
i E A unhappy do
z + 7,5
Happy agents stay put
29:
end for
30:
end for
And set of stochastic update rules:
Alive
(move)
Alive
Pregnant
(reproduce)
Dead
(die)
All transitions are governed by matrix V, which specifies deterministic rules as probabilities equal and defines stochastic transitions through probabilities that depend spatial inter-species interactions. We report the values of the matrix Vi that define our four experimental settings in Tables 5-8. For instance, Prey that interacts with Predator more likely die compared to Prey that does not interact with Predator. Similarly, a Predator that does not interact with a Prey is more likely to die compared a Predator that does interact with Prey, is more likely starve. These spatial interactions are defined by proximity: an agent interacts with others located in its Von Neumann neighborhood (i.e., the neihhooooooo..
20

---

Table 5: Transition matrix V1
Die
Move
Turn pregnant
Turn alive
Stay dead
Stay unborn
Alive Pred Prey
0.15
0.45
0.40
0.00
0.00
0.00
Alive Pred + No prey
0.25
0.55
0.20
0.00
0.00
0.00
Alive Prey Pred
0.30
0.45
0.25
0.00
0,00
0,00
Alive Prey + No pred
0.15
0.40
0,45
0.00
0.00
0.00
Pregnant Unborn child
0.00
0.00
0,00
1.00
0.00
0.00
Dead
0.00
0.00
0.00
0.00
1.00
0.00
Unborn + Npp"
0.00
0.00
0.00
0.00
0.00
1.00
Table 6: Transition matrix W2
Die
Move
Turn pregnant
Turn alive
Stay dead
Stay unborn
Alive Pred + Prey
0.35
0.45
0.20
0.00
0,00
0,00
Alive Pred + No prey
0,25
0.60
0.15
0.00
0,00
0.00
Alive Prey Pred
0.45
0.50
0,05
0.00
0,00
0,00
Alive Prey + No pred
0.35
0.35
0.30
0.00
0.00
0.00
Pregnant Unborn child
0.00
0.00
0.00
1.00
0.00
0,00
Dead
0.00
0.00
0.00
0.00
1.00
0.00
Unborn Npp"
0.00
0.00
0.00
0.00
0.00
1.00
Table 7: Transition matrix 43
Die
Move
Turn pregnant
Turn alive
Stay dead
Stay unborn
Alive Pred + Prey
0.15
0.30
0.55
0,00
0.00
0.00
Alive Pred No prey
0.30
0.55
0.15
0.00
0,00
0,00
Alive Prey + Pred
0.70
0.20
0.10
0.00
0.00
0.00
Alive Prey No pred
0.10
0.40
0.50
0.00
0.00
0.00
Pregnant Unborn child
0.00
0.00
0.00
1.00
0.00
0.00
Dead
0.00
0.00
0.00
0.00
1.00
0.00
Unborn + Npp"
0.00
0.00
0.00
0.00
0.00
1.00
Table 8: Transition matrix V
Die
Move
Turn pregnant
Turn alive
Stay dead
Stay unborn
Alive Pred + Prey
0.15
0.35
0.50
0.00
0.00
0.00
Alive Pred + No prey 0.25
0.45
0.30
0.00
0,00
0,00
Alive Prey + Pred
0.45
0.40
0.15
0.00
0.00
0.00
Alive Prey + No pred
0.30
0.40
0.30
0.00
0,00
0,00
Pregnant Unborn child
0.00
0.00
0.00
1.00
0.00
0,00
Dead
0.00
0.00
0.00
0.00
1.00
0.00
Unborn + Npp"
0.00
0.00
0.00
0.00
0.00
1.00
Not pregnant parent.
21

---

end for Algorithm Predator- Prey model
Require: agent A (1.,......) n). grid size L, transition matrix w. max steps
......,,,,,,,...................
1: for all E A do
2: K()) ~ Uniform({ Prey, Predator}) })
b Assign kind Initial phase
3:
(()) ~ Uniform( (Alive, Unborn})
4:
if <()) = Unborn then
5:
parent(i) ~ Uniform (1j E A A: A:j = k(D)))
Assign parent
6: end if
7:
Zo
(K))) (())
50 )
Agent state
8:
end for
9:
fort=
do for all i E A do
T ......
10:
11:
switch f( case Alive
Alive agent dynamics
13:
N19) Von Noumann(() (0) 2)
Von Neumann neighborhood
14:
neighbor 3j E A: (xl) () E ENO a KG) K(0 Find opposite-kind neighhors
Jet1 Categorical (K((, f neighbor))
Random life phase update
f((] = Alive
17:
(z (8)
Uniform(A7.)
Move randomly
18:
else if f( Dead then
19:
(x Mi+1 + 0)
Agent dies
20:
else f1-1 Pregmant then (1 De+1) - (a/ v()
21:
(6)
b Stay in place
22:
end if
23:
Pregnant
Birth
24:
p()) Alive
1et1::Be41 (x(()) w/.
Stay place
26:
case Dead
27:
s++1 Dead
p(i)
Remain dead
28:
(0 Mi41) (8,0)
(2) (8)
29:
Unborn
30:
j + (1)
Get parent
(()) Pregnant then
((t) Alive
b Born if parent pregnant
33:
N1 < VonNeumann (1) (2
34:
(51+1- ~ Uniform((V) (NV (1)
(4) (6)
D Place near parent
35:
else
36:
8441 Undorn
(1)
b Remain unborn
37:
(x Bet1) + (0,0)
38:
end
39:
end for
40:
if Va. A, f( Dead, Unborn) then
41:
break
Simulation ends
end
43:
22

---

C Further experimental details and results
In this section, we provide additional details on the experiments and present further results. In the subsection on experimental design (C.1), we describe the micro and macro metrics used, including the number of points over which these metries are computed. We then present additional qualitative results on reproducing emergent segregation in the Schelling model (C.2) and emergent oscillations predator-pre dynamies (C.3). Then, we provide further quantitative results for both models (C.4). Finally, we provide discussion on macro-level baseline (C.5), specifically an autoregressive model of order one, or AR(1), a standard time-series model.
C.1 Experimental design
Experiment details In our experiments, we considered three combinations of the parameter & for the Schelling ABM and four combinations of the matrix V for the Predator-Prey ABM. For each parameter setting, we generated 8 training datasets obtained with different initial seeds and trained surrogate model and two ablated models for each. In total, we trained 168 models, 96 (8 x 4 x 3) for Predator-Prey, and 72 (8 x 3 x 3) for Schelling. All our evaluations done across these models per parameter configuration.
We fixed some of the ABM parameters across experiments, which are reported in Table 9. For Predator-Prey, density refers to the density agents that are initialized as Alive, whereas the number of agents refers to the total number of agents in the system (Alive, Dead, Pregnant, and Unborn agents). Color Ratio indicates the ratio between the number of black and white agents; Kind Ratio indicates the ratio between the number of predators and the number of preys.
Table 9: ABM parameters in our experiments
Component
Schelling
Predator-Prey
Grid size
51 51
32 32
Density
0.75
0.3
Agents
1950
2048
Color/Kind Ratio
1:1
1:1
Micro evaluation. We evaluate how well the surrogate captures individual behavior of agents on a future ramification dataset of T = 25 timesteps. We generate this initial condition the last system configuration Zy-1/r i = 0] from the training ramification dataset.
dataset by giving as
Thus, for each agent E A have 24 initial conditions (z(²) (z(")
JEN(() and 500 outcomes L+1 to build 24 ground truth probability distributions, such as equation (1). Then, we use our surrogate model to generate 500 outcomes 2 given as condition (z(") { Jev())), producing 24 probability distributions such as (1) for each agent i A, which we compare to the probability distributions from the ground truth. For Schelling, we measure the EMD on the marginals of the
coordinates x and y. For Predator-Prey we measure the EMD on the distributions of the phases f(() and fix the distance between the different phases to 1, since they represent a categorical variable. For each of the 8 experiments, we calculate the mean EMD over all agents and all timesteps. For Schelling, we evaluate the EMD on 1950 distributions (one per agent) for the x coordinate and 1950 distributions for the y coordinate, over 24 timesteps, yielding 93600 EMD entries. For Fredator-Prey we evaluate the EMD on 2048 distributions (one per agent) over 24 timesteps, yielding 49152 EMD entries. The box plots in Figure 4 have entries the 8 mean EMD values for the surrogate, 8 mean EMD values for the diffusion-only ablated model, and 8 mean EMD values for the GNN-only ablated model.
Macro evaluation. We evaluate how well the surrogate model reproduces system- level behavior by tracking summary statistics over time. We generate 100 independent simulations of 25 timesteps beyond the training horizon for each model, and compute the symmetric mean absolute error (sMAPE) between the mean ground-truth trajectory and the mean surrogate- -eeeertteeeee trajectory across these 100 independent simulations. In the case of Schelling, where we only track the number of happy agents, the sMAPE definition is straightforward. Let A, and F; be the mean number of happy agents
23

---

across the 100 independent simulations from the ground-truth and the surrogate model respectively. Then. the sMAPE computed as:
T A Fd A
sM APE
(6)
T r=1
For Predator-Prey, we track the number of predators and preys on the grid, which follow two distinct trajectories. Thus, for each kind we apply Formula 6 separately and get sM Epreys and sM AP preedtorrs. Then, work with a single value, we calculate the mean value:
sMAPE prodprcy
ppedppeyy
(sMAPE preys sMMPPE PPPPPPeerrrssss predators)
For each the 8 experiments, we calculate the sMAPE over 25 timesteps and then compute its mean. The box plots in Figure 4 have entries the mean sMAPE values for the surrogate, 8 mean sMAPE values for the diffusion-only ablated model, and 8 mean sMAPE values for the GNN-only ablated model.
C.2 Reproducing emergent segregation
Figure in the main text illustrated how, for three selected time steps (t = 0, = 15, and = 30) and three different values of the intolerance threshold Ç. our surrogate model successfully reproduced the qualitative dynamics of the original ABM. contrast, the ablated models failed to capture these dynamics.
Figures 8, 9, and 10 provide more detailed view of the system's evolution for the values &1, {2. and Sa. respectively. In addition to the previously shown snapshots, these supplementary figures include intermediate time steps (t 5, 10, 20, and = 25), offering a more complete picture dynamics.
Original A81
Surrogate
Diffusion only
GNN-only
Figure 8: Evolution of the position of black and red agents in the Schelling model, for tolerance thresholds { = 0.625. We compare the ground truth (top row) with our surrogate (second row) and the ablations (third and fourth row).
Figure illustrates that the original ABM rapidly evolves toward a segregated configuration, which becomes visible as early as ( = 5. However, the resulting clusters remain relatively small, indicating that low intolerance threshold allows for a moderate level of social mixing. This motivates the label "Low Segregation" for the & parameter setting in Figure 4. The surrogate model accurately reproduces
24

---

these qualitative patterns, albeit with slightly larger clusters. In contrast, the diffusion-only ablation fails to capture the emergent spatial structure, remaining close the initial random configuration across time steps. Instead, the GNN-only ablation ends up overlapping most agents on the same coordinates.
1=0
Original
Surrogate
Diffusion only
GNN--nnly
Figure 9: Evolution of the position of black and red agents in the Schelling model, for tolerance thresholds & 0.75. We compare the ground truth (top row) with our surrogate (second row) and the ablations (third and fourth row).
Figure 9 shows that, for E E2, the original ABM also converges toward segregated state, but the convergence occurs more slowly than in the & case. The resulting clusters are larger and more distinct, reflecting a "High Segregation" scenario. The surrogate model again successfully replicates these dynamics, while the ablation models continue to exhibit no structured behavior.
30
Original
Surrogate
Diffusion only
GNN-only
Figure 10: Evolution of the position of black and red agents in the Schelling model, for tolerance thresholds { = & 0.875. We compare the ground truth (top row) with our surrogate (second row) and the ablations (third and fourth row).
25

---

Finally, Figure 10 demonstrates that, for = E3, the original ABM does not converge to a stable configuration. Instead, the high intolerance threshold causes agents continuously relocate, pre- venting the emergence of segregated clusters. In this degenerate case, both the surrogate and the diffusion-only ablation model correctly reproduce the persistent disorder of the original dynamies, while the GNN-only ablation still fails.
C.3 Reproducing emergent oscillations in predator-prey dynamics
Figure 3 in the main text compares population trajectories from the ground-truth Predator- Prey ABM, our surrogate model, and the ablated model. For both parameter sets shown (L1 and 1 the surrogate accurately reproduces the stochastic dynamics beyond the training window, while the ablation fails to capture the key oscillation patterns.
Figure 11 extends the previous analysis to parameter sets 2/2 and V. which produce, respectively, monotonic decline in prey and predator populations, and oscillatory dynamics for preys but not for predators. The case of V particularly illustrative: despite its apparent simplicity, the ablated models fail to reproduce the monotonic trends, showing an unrealistic spike in both populations immediately after the training window in the diffusion-only ablated model, and a roughly constant number of preys and predators in the GNN-only ablated model. In contrast, the surrogate model accurately captures the expected decay. Under v3, the surrogate successfully replicates the oscillations in the prey population and the stable predator trend, while the diffusion-only ablation once again outputs generic dynamics, largely insensitive to the underlying regime, and the GNN-only generates an oscillating dynamic for both. This highlights the ablation models' inability to distinguish between qualitatively different behaviors.
Beyond aggregate population counts, it is instructive to examine the spatial distribution of predators and preys over time. Prior studies [45, 46] have shown that similar predator-prey models give rise rich spatial dynamics, characterized by the spontaneous emergence of structured patterns (see, e.g., Figure 1.3 in [47]). Starting from initially random configurations, the interactions between agents give rise to both short- and long-range range spatial correlations, with predators and preys organizing into dynamic clusters and propagating waves. These patterns are reminiscent of those observed in spatially extended faction-ddifusion systems and bear a strong resemblance to spatial chaos phenomena in evolutionary game theory [27].
We observe similar spatial patterns in our predator-pre preyy ABM. Figures 12, 13, 14, and 15 show the evolution of the positions of predators and preys on the grid.
Figure 12 illustrates the spatial dynamics under parameter set = V,, which is characterized by oscillations in predator and prey populations. The ground-truth model displays rise in population densities mid-simulation, followed by a near-extinction phase toward the end, in line with the temporal trends shown in Figure 3. Notably, the ground-truth also exhibits distinct spatial patterns, with predators and preys forming dynamic clusters. The surrogate model successfully reproduces both the population dynamics and the emergent spatial structures, whereas the ablated models fail to capture any meaningful spatial organization.
Similar results are observed in Figures 13, 14, and 15, which capture different dynamic regimes for predators and preys. In all cases, the surrogate model accurately reproduces both spatial and temporal patterns, while the ablated models consistently fail to capture the underlying dynamics spatial structure.
26

---

22 Monotone predators and preys
Number aggnnts
200
Ground truth
Sumrggate model (Ours)
Diffusio-only
GNN--nnl
Preys
150
100
50
3: Oscillations preys, monotone predators
Ground
Surrogate
Diffusion-only
GNN-only
800
Predators
Number of
200
20
30
10
20
30
10
20
30
30
Time
Tirne
Time
Time
Figure 11: Forecasting macro-level summary statistics (here, the number of alive preys and predators time), starting from the last condition in training, for 100 independent simulation runs, under configuration 2/2 (monotonic dynamics for both predators and preys, top) and 3/3 (oscillations only for predators, bottom). First colummn: original ABM simulations. Second column: surrogate. Third and fourth column: ablations. The dashed vertical line indicates the end of the training phase for the surrogate and ablation models.
GNN-only Diffusion-only Surrogate Original ABM
30
Figure 12: Evolution of the position of preys (black) and predators (red) the predator-prey ABM. for parameters V v We compare the ground truth (top row) with our surrogate (second row) and ablations (third and fourth row).
27

---

Original ABM
t=0
t=5
t= 10
t= 15
t= 20
25
30
GNN-only Diffusion- only Surrogate
Figure 13: Evolution of the position of preys (black) and predators (red) in the predator-prey ABM, for parameters T = 1/2. We compare the ground truth (top row) with our surrogate (second row) and the ablations (third and fourth row).
Original ABM
t=0
5
t= 10
t= 15
t= 20
== 25
30
GNN-only Diffusion- only Surrogate
Figure 14: Evolution of the position of preys (black) and predators (red) in the predator-prey ABM, for parameters W V We compare the ground truth (top row) with our surrogate (second row) and the ablations (third and fourth row).
28

---

GNN-only Diffusion- only Surrogate Original ABM
t= 5
t= 10
t= 15
t 20
= 25
30
Figure 15: Evolution of the position of preys (black) and predators (red) in the predator-prey ABM, for parameters V = v We compare the ground truth (top row) with our surrogate (second row) and the ablations (third and fourth row).
29

---

C.4 Quantitative results
In this section, we present additional details regarding the quantitative evaluation discussed in Section 4.2 (page 8). In particular here we detail further the micro-level evaluation; that is, how much our surrogate model can reproduce the distribution of possible states of an agent at (that is, Z, given its past state Z To better understand how this comparison works, Figure 16 (left side) shows the distribution for single agent of the position z; one of the components of the state Z for each of the two possible past conditions in the Schelling model, happy and unhappy. We see that when the agent is happy, its position remains fixed, while in the case the agent is unhappy, moves randomly with certain distribution peaked around the starting point (implicitly defined by the ABM, see Algorithm 2). Our surrogate model aims at reproducing this distribution, without knowing the original ABM and without access to the latent variable happylunhappy, but only observing the sequence of states and the graph of interactions. The distributions obtained by our surrogate in the same starting conditions reconstruction, we quantify it as the Earth Mover's Distance between the distributions, for each agent in each timestep, and then aggregate these measures.
shown the right side of the figure. To evaluate the quality of this
The result of this comparison shown in Figure of the main text. Here, we present additional experiments showing how these results change across different conditions.
Schelling model. Figure 17 shows the distribution of the EMD scores obtained by our surrogate model and by the ablation model only for the stochastic rules of the Schelling model. In this ABM, the stochasticity lies in the random movement of the unhappy agents (Algorithm 2, lines 16-25). We see that our surrogate can capture these distributions well (EMD averages below 5, considering that the scale given by the length of the grid, L 50) and better than the surrogate model.
Predator-pres prey model. Figure 18 instead shows the distribution of the EMD scores only for the stochastic rules of the Predator- Prey model. Here, the stochasticity lies in the Alive phase, where the agent might transition to another phase depending on its interaction with its neighbors. We see that also in this model our surrogate able to capture the stochasticity well, with EMD averaging below 0.1 for each of the four configurations 1 By contrast, the diffusion-only ablation model's EMD for these same stochastic rules averages above 0.12, up to values around 0.2. Here, the scale is in (0, 1] since this eтror is measured on the life phase binary vector. We further investigate these results in Figure 19 by dividing them by agent state each configuration. Here, the stochastic transitions happen only in the first column (Alive). This figure confirms that the EMD values are quite low, and lower than those obtained by the ablation model. Moreover, here we observe that also the deterministic transitions are recovered with precision (error is always below 10 10-3)) by our surrogate model.
30

---

Ground truth
Surrogate model
1.00
1.00
Probability
0.75
0.75
0.50
0,50
0.25
0.25
0.00
0.00
-20
10
10
(a) Happy agent
Ground truth
Surrogate model
0.06
0.06
Probability
0.04
0.04
0.02
0.02
0.00
0.00
-20
-10
-20
-10
20
(b) Unhappy agent
Figure 16: Example distributions of the position of a single agent the Schelling model, in two different conditions: (a) happy agent; (b) unhappy agent. Histograms the left represent the original ground-truth ABM, those on the right the ones obtained by our surrogate model. The dashed red vertical line indicates the initial position of the agent. Note that the coordinates have been rescaled to [-25, 25), compared to the [0, 50] interval used for ease of exposition in Algorithm 2.
Surrogate
Surrogate &
Surrogate Ft
0.0
0.4
Density
Density
0.3
Density
0.2
0.4
0.0
0.0
0.0
10
15
20
10
15
20
10
15
20
&a
&s
1.00
Dersity
0.15
Dersity
0.75
0.10
0.50
0.05
0.25
0.0
0.00
EMD
EMD
EMD
Figure 17: Distribution of errors (EMD) obtained by our surrogate model (top row) and by the diffusion-only ablation model (bottom row) only for the stochastic rules the Schelling ABM for each considered configuration of the parameter & of the ABM. Different colors indicate independent experiments that only differ by the random seed used to generate the ground truth data.
31

---

Surrogate
Sumogaate
Surrogace
Surnogate
2.7
0.3
2.4
0.6
6.8
100
0.20
02
0.4
9.6
0.8
1.0
0.2
0.4
2.5
0.8
1.0
0.2
0.6
0.6
9.8
1.0
Ablacion
Abbidon W
Abladion %
Ablatton
10
Dansty
Dendly
Dendyy
0.0
0.6
0.2
1.0
Figure 18: Distribution of errors (EMD) obtained by our surrogate model (top row) and by the diffusion-only ablation model (bottom row) only for the stochastic rules of the Predator-Prey ABM for each considered configuration of the parameter matrix V of the ABM. Different shades indicate independent experiments that only differ by the random seed used to generate the ground truth data.
Surrogate
Ablation
0.35
0.064
8.7e-06
0.00013
0.11
1.2e-07
1.3e-05
0.34
0.30
0.25
0.061
4.3e-05
7e-07
0.12
1.3e-07
5.4e-05
0.095
0.20
0.15
0.084
1.4e-05
6.4e-05
0.2
2.6e-08
5.7e-06
0.23
0.10
0.058
2.4e-06
5e-05
0.074
8.4e-08
5.1e-06
0.05
0.00
Alive
Dead
Pregnant
Unborn
Alive
Dead
Pregnant
Unborn
Figure 19: Mean error (EMD) obtained by surrogate (on the left) and by the diffusion-only ablation model (on the right) for each configuration of the Predator-Prey Preutor-Prry ABM (on the rows, 11-1) and for each initial state of an agent (on the columns).
C.5 Macro baseline (AR1)
In this subsection, show that simple yet commonly used time-series baseline the autoregressive model of order one (AR(1)) - completely fails to reproduce the macro-level dynamics of our system. We focus on the predator-prey model, although qualitatively similar results hold for the Schelling model and for higher-order time-series models (e.g., AR(2)).
The AR(1) model assumes that the value of a time series step depends linearly its previous value:
d0r Ct.
where & is the autoregressive coefficient and _N00, o2) is Gaussian noise. We fit the model separately to each macro variable using least squares estimation on the training portion of the simulated data. The fitted coefficients and Ô are then used to generate out-of-sample forecasts, producing trajectories that can be directly compared with those generated by our surrogate model.
As shown in Figure 20, the AR(1) baseline simply extrapolates the last observed trend in the training data, completely failing to reproduce the non-monote non-monotonic, oscillatory patterns characteristic of the underlying dynamics - and even performing poorly in capturing monotonic trends. This limitation arises because the system's 's behavior is defined at the micro level, while macro-level variables are
32

---

only aggregate summaries of those micro interactions. Consequently, surrogate model that learns from micro-level states, ours does, is inherently better positioned to capture both the micro and the emergent macro dynamics.
Ground truth
Surrogate model (Ours)
AR(1) model
Number of agents
500
Preys Predators
400
300
200
100
0
20
30
10
20
30
10
20
30
Time
Time
Time
(a) 1
Ground truth
Surrogate model (Ours)
AR(1) model
Number of agents
200
Preys
150
Predators
100
50
0
Time
Time
Time
(b) V/2
Ground truth
Surrogate model (Ours)
AR(1) model
Number of agents
400
Preys
300
Predators
200
100
0
20
30
10
30
10
20
Time
Time
Time
(c) V
Ground truth
Surrogate model (Ours)
model
Number of agents
400
Preys
300
Predators
200
100
%
10
20
30
10
20
30
10
20
30
Time
Time
Time
(d) V
Figure 20: Comparison of - under the surrogate GDN model and an AR(I) baseline.
SM References
[45] André De Roos, Edward McCauley, and William Wilson. Mobility versus density-linme -imited predator-prey prry dynamics different spatial scales. Proceedings of the Royal Society of London. Series Biological Sciences, 246(1316):117-12) 1991.
[46] Douglas Donalson and Roger M Nisbet. Population dynamics and spatial scale: effects of system size on population persistence. Ecology, (M8)):442-25077 1999.
[47] Volker Grimm and Steven Railsback. Individual-basee modeling and ecology. Individual- based modeling and ecology. Princeton university press, 2013.
[48] Alexander Quinn Nichol and Prafulla Dhariwal. Improved denoising diffusion probabilistic models. In Proceedings of the 38th International Conference on Machine Learning, volume 139, pages 8162-8171. PMLR, 2021.
33

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Yes, the claims made in the abstract and introduction accurately reflect the paper's contributions and scope. We provide the methodology and framework in $3 and empirical results in $4.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Limitations are discussed in the Discussion $5.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, c.g., if the approach was only tested few datasets with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-tex system might not be used reliably provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
34

---

Answer:
[NA]
Justification:
The paper does not include theoretical results.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: The paper describes the architecture clearly and fully in $3, and the steps to reproduce the results are described in the Supplementary Material.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, descríbing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is new model (e-g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
5. Open access to data and code
35

---

Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Code and data is available in the Supplementary Material, and instructions reproduce the results are given.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.ce
httpesuamisssiooo11iooooo Polliyyy) for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: / /ttpps.ccpppliiiiim wwtppccppccc///ccccccccccccccclt htpollc/ggiieee http//uuddeee/eee
for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The experimental settings are described in a higher level of detail in $4.1 and in full detail in the Supplementary Material.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: The box plots of Figure 4 include the statistics of all quantitative results presented the paper.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
36

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
OK to report 1-sigma error bars, but one should state The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: $4.1 provides the computational resources needed reproduce the experi- ments. Additional details about the time needed to run experiments are included in the Supplementary Material.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https: / /ttp//nuuupuppeee cttpsscccppbblie hetpubli//tthhssuudeeeeee
Answer:
[Yes]
Justification:
We followed the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g.. if there special consid- eration due to laws or regulations in their jurisdiction..
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: There no societal impact the work performed. As discussed in $2, ABMs have traditionally been powerful for theory generation, and recently, they have become increasingly data-driven, but our work focuses solely on the methodological advance of learning differentiable ABM surrogates.
37

---

Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied partieular applications, let alone deployments. However, if there is direct path any negative applications, the authors should point it out, For example, it is legitimate to point that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI...
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [NA]
Justification:
The paper does not use existing assets.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
38

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, paperswi thcode. com comddataees has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mmodel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
39

---

Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if' applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: LLMs were used only for editing/formatiing.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy ptttps:/neeurpp. cc///oneeeeecccs///2//////////// for what should or should not be described.
40


Reviews:
[{"review_id": "JBxjviPiqk", "reviewer": "Reviewer_BiVF", "summary": "This paper addresses an important yet challenge problem making Agent-Based Models (ABMs) differentiable, particular by making the behavior of each agent differentiable. The authors propose Graph Diffusion Networks (GDNs), a novel combination of GNN and diffusion models to learn differentiable ABM surrogates. The GNN capture the interactions of agents, and the diffusion model for stochasticity behavior. The method is evaluated on the Schelling segregation model and a Predator-Prey ecosystem, showing strong fidelity in reproducing both micro-level (individual) and macro-level (aggregate) behaviors.", "questions": "See weaknesses point 2 for why not compare using other existing methods.", "limitations": "The author discussed several key limitations. See weakness point 1 for limitation of scalability", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "The paper is very well written and presents a clear and comprehensive review of related work. It is well structured, with a strong motivation and background. The proposed method is novel and demonstrates strong experimental performance compared to the baseline. The limitations of the current approach are also thoughtfully discussed.\n\n### weaknesses:\n\n(1) The paper does not clearly state how large the ABMs used in the experiments are — there's no mention of the number of agents, grid size, or how dense the interactions are. This information is missing from the main text and, if only available in the supplement, should really be included in the main paper. From the figures and setup, it seems these are small-scale ABMs with less than 500 agents. \nThat’s important, because the proposed method involves several computationally intensive components: diffusion models with many sequential steps, a GNN applied per agent, and a ramification process that generates hundreds of alternate futures per timestep. These choices likely make scaling to larger systems quite costly.\nSince the primary motivation for making agent behavior differentiable is to enable learning and data fitting in complex, realistic ABMs, a method that scales only to small systems may fall short of this goal.  A discussion of computational efficiency is needed to understand where and how this method could realistically be used.\n\n(2) While the authors focus on individual-level modeling, both ABMs used (Schelling and Predator–Prey) have well-understood emergent dynamics that could also be modeled by existing surrogate approaches. The decision not to compare against other surrogate models -- as discussed in the background section -- is not fully justified. Even if these baselines do not model per-agent transitions, they are appropriate for comparison on the macro-level forecasting tasks Including such baselines would strengthen the empirical claims and clarify the added value of modeling individual behavior explicitly.", "was_revised": true, "final_justification": "I appreciated the feedback from the author, I will keep my original rating", "rebuttal": "We thank the reviewer for the thoughtful and constructive feedback. We reply by numbering the weaknesses by W1, W2, … and the questions by Q1, Q2, …, grouping weaknesses and questions where relevant. For each reply, we start with a one-sentence TL;DR and then develop our argument.\n\n---\n# W1: Scaling\n**TL;DR**: Our surrogate model has been tested on mid-sized ABMs, the training time scales linearly with number of ramifications and timesteps, and sublinearly with number of agents.\n\nWe agree that the main manuscript should state the concrete dimensions of the testbeds and clarify how the method scales.  We will move the numbers that are currently in the supplement to the main text (§4) and add the following discussion.\n\n| Component           | Schelling | Predator‑Prey |\n| ------------------- | --------- | ------------- |\n| Grid size           | 51 x 51   | 32 x 32       |\n| Initial density     | 0.75      | 0.30          |\n| Number of agents    | 1,950     | 2,048         |\n| Colour/Kind ratio | 1:1     | 1:1         |\n\n\nThese settings are in the range of mid-sized ABMs: each training epoch processes about $\\sim 9\\times10^{6}$ agent transitions (9 time steps x 2000 agents x 500 ramified futures).  On a single RTX-A6000, one epoch takes ≈ 37 seconds.\n\nThe reviewer is correct that an ABM containing millions of agents can stress any learning algorithm. As we have mentioned in the discussion with Reviewer McsJ, the training time of the model grows linearly with the number of ramifications provided to the network during training, as well as the number of timesteps. On the other hand, increasing the number of agents in the system (for example by increasing the grid size and keeping the density of agents constant) increases the training time sublinearly. In fact, with mid-sized systems, we can train the network with batches that correspond to the whole set of agents. Thus, the number of training iterations per epoch does not change, but rather the size of the batch in input to the network. Considering the Predator-Prey model, if we increase the grid size from 32 to 64, going from 2048 agents to 8192 (4 times as many agents), we have an increase from 37.09 s/epoch to 108.37 s/epoch, which is 2.92 times larger. Similarly, by increasing from grid size 32 to 48, we increase by a factor of 2.25 the number of agents, but the training time increases by a factor of 1.72. Regarding inference, the computational cost to produce the next agent state is always the same. However, it will take longer for larger systems to simulate the evolution of the entire system. Below we show a table with the training time and inference time for the 3 settings of grid sizes we have discussed, where inference time refers to the time necessary to generate 1 timestep of a simulation of N. agents with the trained GDN. All other hyper-parameters are the same as in the paper experiments. We include the evaluation metrics at the macro and micro level used in the paper, to assess that even with a larger system the dynamics are learned by our methodology.\n\n| Grid size |  N. agents | Time/epoch |Inference time |  sMAPE   |    EMD    |\n|-------------|---------------|---------|--------------|-----------|-------------|\n|      64     |      8192     | 108 s |  1.6 s/it. |0.041  |  0.0080 |\n|      48     |      4608     | 64 s   |  0.7 s/it. |0.044  | 0.0084  |\n|      32     |      2048     |  37 s  |  0.3 s/it.|0.079*|  0.0091*|\n\n*: Mean values from the paper micro and macro evaluations. \n\nTraining the model with a ramified dataset of R ramifications is a choice made to allow the model to be exposed to many agent transitions ($\\sim 9\\times10^{6}$ in our experiments). Increasing the number of agents in the system increases the number of agent transitions available in the dataset, expanding the statistical coverage provided to the model to learn. In our previous experiment, the models with higher numbers of agents performed better both at the macro and micro level than in the original experiments. The number of ramifications R could be decreased to potentially get the same result quality of the paper experiments. We ran an experiment by adjusting the number of ramifications R according to the number of agents so that we have a similar number of agent transitions in the training set as in the original experiments. Given the time constraints of the rebuttal period, we were able to train one model for each grid size configuration on one dataset of Predator-Prey with parameter $\\Psi_1$ and evaluate the model as done in the paper. The results are shown in the table below. We find that the sMAPE is around the value found in the original experiments (actually lower for grid size 48), and the mean EMD is slightly higher, but close to the mean value from the original experiments. We will add the clarifying discussion above and additional scaling experiments with more hyperparameter values to the revised manuscript and supplement material.\n\n\n| Grid size |  N. agents |  N. ramifications  |  Time/epoch  |    sMAPE   |    EMD    |\n|-----|----|-----|--------|----|----|\n|      64     |      8192     |    125     |   27 s |     0.076     |  0.0121   |\n|      48     |      4608     |     250      |   32 s |     0.033     |  0.0104   |\n|      32     |      2048     |     500      |   37 s |     0.079*    |  0.0091* |  \n\n*: Mean values from the paper micro and macro evaluations. \n\nInference time refers to the time needed to produce 1 timestep of a simulation of N agents with the trained GDN. \n\nRegarding scaling to even larger systems (e.g. in the range of $10^4$ to $10^6$), it would be more appropriate to train the model on batches made of subgraphs of the system configuration, as it would become troublesome to optimize the model on batches of tens of thousands or millions of samples. However, apart from batch optimization, the methodology would remain unchanged. \n\n---\n# W2, Q1: Macro-level baselines\n**TL;DR**: While surrogate models for aggregate dynamics do exist, they struggle to reproduce trends that are absent from the training data. In contrast, our approach (based on learning individual behavior) can generalize to these unseen dynamics.\nMost existing surrogate approaches do not actually learn system dynamics. Instead, quoting from our manuscript, surrogate models typically learn directly the mapping from parameters $\\Theta$ to static summary statistics, disregarding individual behavior and model dynamics.  \n\nThere are a few exceptions. For example, Fadikar et al. [1] use functional Gaussian Processes to capture the shape of epidemic curves, while Kieu et al. [2] employ LSTMs to forecast time series of pedestrian counts. However, these approaches share two key limitations: (i) they model aggregate variables rather than micro-level behavior, and (ii) they rely on statistical models fitted to in-sample data.  \n\nAs a result, these models perform well when out-of-sample dynamics resemble those seen during training, but their performance deteriorates when the dynamics change. For instance, Fadikar et al. focus on cumulative case counts in epidemics, which follow a monotonic trajectory. In such cases, using early time steps to predict future ones is a reasonable strategy. But this assumption breaks down in more complex settings.  \n\nConsider the predator-prey model shown in Figure 3 (top panels), where both predator and prey populations grow and then decline. These dynamics are far from monotonic. A surrogate trained on aggregate counts during the early phase (when populations are rising) is unlikely to capture the subsequent decline.\n\nTo illustrate this point, we implemented a simple AR(1) process as a representative of aggregate time-series-based models. We trained it on just 10 time steps (matching the surrogate’s training window). Given this limited horizon, using more complex models like LSTMs would be inappropriate.\n\nThe AR(1) surrogate substantially underperforms our individual-level surrogate. For example, in scenario $\\Psi_1$—characterized by oscillatory predator-prey dynamics—the AR(1) either predicts indefinite growth or convergence to a steady state. It fails to capture the true cyclical behavior, resulting in a significantly higher sMAPE. Even in the simpler scenario $\\Psi_2$, where dynamics are monotonic, the AR(1) model incorrectly predicts convergence to a steady state (e.g., a predator and prey population of 60), again leading to large errors.\n\nThe table below reports the sMAPE (mean and standard deviation) across four scenarios. As shown, the AR(1) surrogate performs poorly for all scenarios, with errors far exceeding those of our micro-founded surrogate reported in Figure 4.\n\n| Scenario | sMAPE (mean) | sMAPE (std) |\n|---|---|----|\n| $\\Psi_1$ | 0.82 | 0.08 |\n| $\\Psi_2$ | 0.82 | 0.06 |\n| $\\Psi_3$ | 0.44 | 0.06 |\n| $\\Psi_4$ | 0.96  | 0.12|\n\nIn summary, we deliberately chose not to include comparisons with macro-level dynamic surrogates in the main paper because their limitations are structural: they rely entirely on aggregate time series observed during training. As such, they are fundamentally incapable of extrapolating to qualitatively different out-of-sample dynamics—precisely the type of generalization our setting requires. The AR(1) results shown here support this point, illustrating that even the simplest non-monotonic patterns are missed. For this reason, we maintain that comparing our surrogate to these macro-level approaches would not provide a meaningful benchmark. \nWe will include this discussion and the results on the AR(1) in the revised paper.\n\n[1] Fadikar, A., Higdon, D., Chen, J., Lewis, B., Venkatramanan, S., & Marathe, M. (2018). Calibrating a stochastic, agent-based model using quantile-based emulation. SIAM/ASA Journal on Uncertainty Quantification, 6(4), 1685-1706.\n\n[2] Kieu, M., Nguyen, H., Ward, J. A., & Malleson, N. (2024). Towards real-time predictions using emulators of agent-based models. Journal of Simulation, 18(1), 29-46."}, {"review_id": "y31hQiHQG8", "reviewer": "Reviewer_McsJ", "summary": "The authors propose a novel framework, the Graph Diffusion Network (GDN), that  combines graph neural networks and diffusion models to learn a fully differentiable surrogate model directly from ABM-generated simulation data. The core innovation is a shift in focus from emulating macro-level system outputs to directly modeling the behavior of individual agents. The framework is validated on two canonical ABMs: Schelling's segregation model and a Predator-Prey ecosystem. The experiments demonstrate that the learned surrogate not only replicates the conditional next-state distributions at the individual agent (micro) level but also accurately reproduces the emergent, system-level (macro) dynamics, even when forecasting beyond the training period. The results show a significant improvement over an ablation variant that does not model the interaction graph, highlighting the importance of the proposed architecture.", "questions": "- The paper should include a discussion on the computational costs and scalability of the proposed framework. How the method scales with increasing agent count, time horizon, or stochastic branches per step. A more detailed discussion would be valuable.\n- To assess the importance of modeling stochasticity, the authors should include the diffusion-removed version of the model in their experiments. Furthermore, comparisons against other differentiable simulators—such as Transformer-based architectures or temporal GNNs—would provide stronger evidence for the effectiveness of the proposed GDN framework.\n- The paper does not clearly state the type of GNN employed (e.g., GraphSAGE, GCN, GAT). Since the choice of message-passing scheme can significantly impact expressiveness and inductive bias, the authors should clarify the GNN variant used and its implications for learning agent interactions.", "limitations": "The paper does not present any apparent risks of negative societal impact.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "## Strengths\n- The proposed GDN architecture is sound, well-motivated, and elegantly combines the strengths of GNNs and diffusion models to tackle the specific challenges of ABMs.\n- The paper is clear and well written\n- The ablation study that removes the GNN component clearly and effectively demonstrates the value of explicitly modeling the graph-structured interactions.\n- The authors provide detailed architectural choices and experiment details, supporting reproducibility. \n\n\n## Weaknesses\n- The proposed framework requires full access to an existing ABM simulator to generate a ramified dataset. This restricts its use  and limits its immediate applicability to real-world applications where the simulator may be unavailable or partially known.\n- The empirical comparison is restricted to a single ablation (which removes the GNN component). The paper does not compare against alternative differentiable surrogates such as Transformer-based simulators [Casert et al. 2024] or neural emulators for cellular automata [Grattarola et al.], making it harder to contextualize the contribution within the growing literature on learnable simulators.\n- The diffusion-removed variant of the model is discussed but not included in the experimental evaluation. Including this variant would help isolate the contribution of the diffusion component. Moreover, in this ablation, the authors should consider comparing against temporal GNN baselines [1], which are a natural alternative for modeling sequential agent behavior. Such baselines would offer a stronger baseline for GDN.\n- The framework requires generating a ramified simulation dataset, which can be expensive. This could limit applicability in large-scale ABMs with many time steps, high agent counts, or other real-time constraints.\n\n[1] Gravina et al. Deep learning for dynamic graphs: models and benchmarks. In IEEE TNNLS 2024.", "was_revised": true, "final_justification": "Given the clarifications on all points and the additional experiments, I am now convinced the proposed method does have specific technical advantages.", "rebuttal": "We appreciate the reviewer’s detailed and constructive comments. We reply by numbering the weaknesses by W1, W2, … and the questions by Q1, Q2, …, grouping weaknesses and questions where relevant. For each reply, we start with a one-sentence TL;DR and then develop our argument.\n\n# W1: Need access to the ABM \n**TL;DR:** ABMs cover lots of different systems, but in principle, the method can be used on data with no existing ABM too, granted that sufficient agent/context pairs and transitions are observed in the data.\n\nAgent-based models (ABMs) are widely used to represent individual behavior and interactions across disciplines, and after more than three decades of development, it is likely that for many real-world datasets, an ABM exists that closely approximates the underlying data-generating process. The central challenge lies in aligning such ABMs with observational data.\n\nThat said, our method is not limited to synthetic data. In principle, it could be applied directly to real-world data, provided that the data includes sufficient repeated instances of similar individuals facing similar interaction structures. Taking the Schelling model as an example: in our synthetic setup, we generated 500 alternative futures for an agent moving to a new location given the characteristics of its neighbors. If we had real-world data from, say, 500 households with matching characteristics and neighborhood compositions, their observed behaviors could serve as a real-world ramifications dataset.\n\nWhile this paper focuses on synthetic data as a necessary first step, we are actively exploring applications to real-world datasets as part of ongoing and future research.\n\n# W2, W3, Q2: Baselines\n**TL;DR:** GNN models completely fail to learn stochastic behavior, and thus lead to much lower micro- and macro performance, while alternative baselines are not general enough for ABMs.\n\nThe reviewer proposes three further possible choices for ablations/baselines:\n\n1. Perform an ablation that removes the diffusion component while keeping the GNN.  \n2. Consider alternative GNN models (e.g., Grattarola et al., temporal GNNs).  \n3. Compare to Casert et al. (2024).\n\nWe tackle each of these one by one.\n\n**1. Ablation without diffusion**\n\nWe conducted additional experiments using the ablated model suggested by the reviewer, which excludes the diffusion component while retaining the GNN. This modification is non-trivial, as the original surrogate relies on gradients from the diffusion model to train the GNN. In the ablated version, we instead minimize the MSE between the true and predicted next agent states. The only architectural change is in the MLP following the aggregation pass: its output dimension is set to the agent state dimension rather than the embedding dimension. While the original GDN uses hidden layers of size [32, 64, 128], the ablated model adopts a symmetrical structure: [32, 64, 128, 128, 64, 32]. We trained for 100 epochs with the Adam optimizer (learning rate 2×10⁻⁵) and mini-batch size of 16.\n\nOur results clearly show that the GNN successfully learns to predict agents’ future states when the underlying dynamics are deterministic. However, it performs poorly in the presence of stochastic rules, resulting in significantly lower overall performance compared to the GDN.\n\nThe results for the predator-prey model are summarized in the tables below (mean±standard deviation over 8 runs). For each parameter, the second and third columns reproduce the results from Figure 4 (GDN and diffusion-only ablation), while the last column reports results from the GNN-only ablation as suggested by the reviewer. We observe that both macro-level and micro-level errors are substantially higher in the GNN-only ablation, and often worse than the diffusion-only model. This is expected, as stochastic rules are dominant in most parameterizations, making the GNN's deterministic forecasting approach inadequate.\n\nPredator-Prey macro-level:\n\n|Parameter |GDN |Diffusion-only |GNN-only|\n|-|-|-|-|\n|\\$\\Psi\\_1\\$ |0.079±0.030 |0.266±0.046 |0.940±0.082 |\n|\\$\\Psi\\_2\\$ |0.323±0.149 |1.204±0.135 |0.917±0.198 |\n|\\$\\Psi\\_3\\$ |0.206±0.073 |0.508±0.084 |0.744±0.235 |\n|\\$\\Psi\\_4\\$ |0.324±0.115 |0.768±0.074 |0.999±0.254 |\n\nPredator-Prey micro-level:\n\n|Parameter |GDN |Diffusion-only |GNN-only |\n|-|-|-|-|\n| \\$\\Psi\\_1\\$ |0.0091±0.0005 |0.0521±0.0024 |0.0750±0.0014|\n| \\$\\Psi\\_2\\$ |0.0013±0.0002 |0.0654±0.0102 |0.0122±0.0022|\n| \\$\\Psi\\_3\\$ |0.0096±0.0008 |0.0812±0.0028 |0.0540±0.0034|\n| \\$\\Psi\\_4\\$ |0.0041±0.0003 |0.0784±0.0030 |0.0386±0.0019|\n\nA slightly different picture emerges in the Schelling model for parameterization $\\xi_1$. Here, nearly all agents remain stationary during the test period, leading to effectively deterministic dynamics. Consequently, the GNN-only surrogate achieves a micro-level EMD of 0.73±0.07, which is comparable to the GDN’s performance (0.71±0.06). However, from preliminary results (mean and std over 6 experiments, instead of 8 due to time) on parameterization $\\xi_3$ , the most stochastic configuration, the GNN-only surrogate performs poorly with EMD (9.73±0.06) against GDN’s performance (2.81±0.11). \n\nIn the final version of the paper we will include the full results from this additional ablation study.\n\n**2. Alternative GNN models**\n\nWe argue that all GNN-based models share the same structural limitations observed in the GNN-only ablation examined above, and therefore exhibit similar performance. This class includes the approach by Grattarola et al., which in fact was applied to deterministic cellular automata, as well as the temporal GNNs discussed in the reference suggested by the reviewer  (Gravina et al. 2024)\n\n**3. Casert et al.**\n\nThe reviewer suggested using the Transformer-based model by Casert et al. as a baseline. While it is the only stochastic emulator we are aware of, it is not applicable to the ABMs we consider. Their method is tailored to physical systems in continuous time and relies on explicitly enumerating allowed transitions between configurations. This is feasible when transitions involve small local changes, as in particle systems, but not in ABMs where agents can make large, unconstrained moves (e.g., in Schelling) or where the number of agents varies over time (e.g., predator-prey). In contrast, our ramifications approach samples from the space of possible configurations without requiring predefined transitions, making it suitable for ABMs but incompatible with the Casert et al. methodology.\n\n# W4, Q1: Scaling\n**TL;DR:** Our surrogate model training cost scales linearly w.r.t. time-horizon and number of ramifications and sublinearly w.r.t. number of agents. Increasing time-horizon and number of agents increases stochastic coverage, so the number of ramifications can be reduced.  \n\nWe acknowledge the need to discuss computational cost and scalability. Training our GDN on 9 timesteps, 500 ramifications, and ~2 000 agents (1 950 Schelling, 2 048 Predator–Prey) takes ≈ 37 s per epoch. Because each step processes one timestep and one ramification, runtime grows linearly with the number of ramifications or timesteps and only sub-linearly with the agent count (see reply to Reviewer BiVF). We use multiple stochastic branches during training to increase the statistical coverage of the ABM micro-dynamics. A new experiment on the predator–prey model (parameter set $\\Psi_1$​) confirms that increasing R improves EMD and sMAPE until they match the main-paper results, while training time scales linearly.\n\n|Ramifications |EMD |sMAPE |Time/epoch |\n|-|-|-|-|\n|50 |0.0159 |0.179 |3.7 s |\n|100 |0.0138 |0.108 |7.4 s |\n|250 |0.0096 |0.088 |18.5 s |\n|500 |0.0091\\* |0.079\\* |37.1 s |\n\n*: Mean values from the paper micro and macro evaluations. \n\nWe also agree that generating hundreds of futures for each state may be prohibitive in systems with many agents or longer time horizons. However, larger systems naturally provide more independent samples of similar conditions. In such cases, fewer ramifications may suffice to achieve similar statistical coverage. To explore the trade-off between data generation cost and predictive quality, we trained GDN on three Predator-Prey datasets (parameter $\\Psi_1$) with increasing agent counts and decreasing ramifications, keeping the number of agent transitions roughly constant. All other hyperparameters remain as in the main experiments. Results are shown below.\n\n|Grid size |Agents |Ramifications |Time/epoch |sMAPE |EMD |\n|-|-|-|-|-|-|\n|64 |8192 |125 |26.94 s |0.076 |0.0121 |\n|48 |4608 |250 |32.32 s |0.033 |0.0104 |\n|32 |2048 |500 |37.09 s |0.079\\* |0.0091\\* |\n\n*: Mean values from the paper micro and macro evaluations. \n\nWe find that sMAPE remains close to the original results (and is notably lower for the 48-grid case), while EMD increases slightly but stays near the reported average. We will include these scaling results and the discussion above in the revised manuscript and supplement to better communicate the computational footprint and scalability of our framework.\n\n# Q3: GNN details\n\nWe use a simple MessagePassing GNN, where each agent’s state is passed as the message, aggregated, and fed into an MLP. As the reviewer has correctly pointed out, the choice of the MessagePassing aggregation scheme can impact expressiveness and inductive bias. The most commonly used aggregations for this setting are sum, mean and max. In our experiments we used “sum” for Predator-Prey and “mean” for Schelling since there is dependency on the degree in its dynamics. Minimal experimental evaluations allow to direct the practitioner towards the one which seems promising for the task at hand. \nIn the future we would like to explore ‘aggregation-agnostic’ approaches, in order to reduce inductive bias. Promising options are relying on Attention weighting or leveraging PowerMeanAggregation or MultiAggregation operators from Pytorch.\n\n---\nWe welcome any follow-up questions and would be happy to provide further details or clarifications."}, {"review_id": "Tnt2FKuaMu", "reviewer": "Reviewer_yYgi", "summary": "This paper proposes a method to predict the behavior of Agent-Based Models (ABMs). Specifically, the authors trained a differentiable surrogate of an ABM, which is a GNN, using ABM-generated data.", "questions": "1. Please take a look at the Weaknesses point 1. One unclear description example is: \n(L30) However, this progress is occurring despite the absence of principled methods to systematically align ABMs with real-world data. While various approaches have been proposed for calibrating macro-level parameters of ABMs [29], there are still no established methods for tuning the micro-level behaviors of individual agents to data.\nQ: Why does \"align ABMs with real-world data\" relate to \"tuning the micro-level behaviors of individual agents to data\"? Doesn't \"calibrating macro-level parameters of ABMs\" belong to \"systematically align ABMs with real-world data\"? Also, this paper also doesn't use real-world data, so what is the advantage of the proposed method compared with previous ones?\n\n2. Could the authors provide the related works and also add some baselines for experimental comparison?", "limitations": "Yes. There are limitation discussions in the Section 5.", "rating": 4, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths\nThis paper explained the experimental details very clearly. \n\nWeaknesses\n\n1. The paper is clear in how the experiments are conducted, however, it is unclear in the motivation, related works and baseline selection. \nMy suggestion is that, the authors can pay more attention to explain 1/ why do we need to model individual agent behaviors, what the benefit is. 2/ explanation of terminologies used to make the paper self-contained, e.g., macro-level and micro-level in the context of Agent-Based Models. 3/ the evidence that the selection of Schelling and predator-prey models are representative enough. 4/ systematic comparison with previous works. \n\n2. In this paper, all the baselines are variants of the proposed models. Also, the authors do not mention the well-recognized baselines in this area, and how previous works are evaluated.", "was_revised": true, "final_justification": "The authors directly and clearly answer my questions. I'm very unfamiliar with the area of ABM, so I can only review from a general machine learning and GNN perspective.", "rebuttal": "We appreciate the reviewer’s detailed and constructive comments. We reply by numbering the weaknesses by W1, W2, … and the questions by Q1, Q2, …, grouping weaknesses and questions where relevant. For each reply, we start with a one-sentence TL;DR and then develop our argument.\n\n---\n# W1, Q1: Clarity\n**TL;DR**: We clarify our terminology and explain why the predator-prey and Schelling models are representative of ABMs.\n\nWe realize that the benefits of **modeling individual behaviors** were not adequately explained.\n\nABMs are based on defining rules at the level of individual agents – a perspective known as the micro-level. These rules deal only with the characteristics and behaviors of single agents, such as their type (e.g., predator or prey) or positions. Micro-level dynamics are governed by these rules, that act only on these individual features, such as updating an agent’s position over time. The purpose of ABMs is often to observe how these simple micro-level rules give rise to complex, emergent properties at the system-wide level, referred to as the macro-level. Examples of macro-level outcomes include summary statistics, such as the total number of surviving agents or the total number of agents that are happy of their neighbors, associated to segregation patterns. \n\n\nWith this background, we can now explain the value of modeling individual behaviors directly. In fact, we constrain our surrogate to model only micro-level rules. This ensures that the surrogate remains structurally faithful to the original ABM and does not rely on shortcuts such as directly predicting macro-level outcomes. To illustrate this, consider the predator-prey model. A macro-level surrogate might emulate the system using equations like the Lotka–Volterra model, which predicts population counts over time without modeling individual agents. In contrast, our surrogate expresses rules in terms of agent-level attributes, such as the type and position of each agent. This approach leads to two major benefits. First, by modeling only micro-level rules, the surrogate stays aligned with the logic of the original ABM. Second, the surrogate is expressed directly in terms of micro-level variables, allowing for real-data alignment with respect to individual agent attributes – something standard surrogates cannot do.\n\nTraditional surrogates, in fact, can only perform macro-level calibration (e.g., tuning a global threshold like $\\epsilon$ in the Schelling model). Our model, however, is differentiable with respect to micro-level inputs, such as the type or position of each individual agent. Thanks to this property, our approach offers practical advantages in integration with real-world data. For instance, when micro-level data is available for a subset of agents (e.g., demographic or spatial attributes), our surrogate can use this information to simulate possible future states. \n\nTherefore, in response to Q1, to faithfully \"align ABMs with real-world data,\" we mean to align the micro-level properties of individual agents to match observed data. Calibrating macro-level parameters only allows to calibrate some key macro-level quantities of the ABM, but its agents remains totally disconnected from any entity observed in the real-world. Consider again the predator-prey model: if we have a real-world data source enumerating some real-world animal positions, aligning this data with the ABM at the micro-level would mean associating one of such real positions to a simulated agent, instead of simply tuning a global parameter for the whole system. In this sense, aligning micro-level states with real-world observations remains an open challenge in the ABM literature, yet it is essential for ensuring that ABMs can be reliably applied in real-world contexts. Our model is the first differentiable surrogate capable of expressing learned stochastic rules in terms of such micro-level attributes. As such, it represents a crucial step toward building tools that can calibrate the micro-level components of ABMs, an often overlooked but essential aspect of systematically aligning ABMs with real-world data.\n\nRegarding the **representativeness of our choice of ABMs**, both the Schelling model and predator-prey models are widely recognized as canonical examples of agent-based models. They are commonly used to introduce the field and are frequently employed in studies that explore calibration techniques and surrogate modeling.\n \nThe Schelling model is often cited as one of the earliest ABMs in the social sciences [1], and has served as the sole testbed for the development of data-driven model extraction methods [2]. Similarly, predator-prey models form the foundation for a broad range of more complex models in ecological research [3, 6], and they have played a key role as testbeds for the development of novel calibration and surrogate techniques. For example, in \"Data assimilation with agent-based models using Markov chain sampling\" [4], the authors test their method for an MCMC-based sampler of possible trajectories on a predator-prey ABM, very similar to ours. It is their only case study. Likewise, in \"The use of surrogate models to analyse agent-based models\" [5], one of the three case studies used to compare macro-level surrogates is a predator-prey model. It is worth noting that in all of these works, the predator-prey model is used as a representative of a broader class of ABMs, likely due to its simplicity, intuitive appeal, and ease of explanation to readers across different disciplines.\n\nWe will add these considerations and further clarifications to the final version of the manuscript.\n\n[1] Squazzoni, F. (2010). The impact of agent-based models in the social sciences after 15 years of incursion. History of economic ideas: XVIII, 2, 2010, 1000-1037.  \n\n[2] Jamali, Ruhollah, Wannes Vermeiren, and Sanja Lazarova-Molnar. \"Data-driven agent-based modeling: Experimenting with the Schelling’s model.\" Procedia Computer Science 238 (2024): 298-305.\n\n[3] McLane, Adam J., et al. \"The role of agent-based models in wildlife ecology and management.\" Ecological modelling 222.8 (2011): 1544-1556.\n\n[4] Tang, Daniel, and Nick Malleson. \"Data assimilation with agent-based models using Markov chain sampling.\" Open Research Europe 2 (2022): 70.\n\n[5] Ten Broeke, Guus, et al. \"The use of surrogate models to analyse agent-based models.\" Journal of Artificial Societies and Social Simulation 24.2 (2021).\n\n[6] Murphy, Kilian J., Simone Ciuti, and Adam Kane. \"An introduction to agent‐based models as an accessible surrogate to field‐based research and teaching.\" Ecology and evolution 10.22 (2020): 12482-12498.\n\n---\n# W2, Q2: Baselines\n**TL;DR**: Comparisons to baselines is difficult given the novel focus of our approach.\n\nAll reviewers raised concerns about the lack of baseline comparisons. We acknowledge this and have expanded our evaluation accordingly. However, our work introduces the first individual-level stochastic surrogate for ABMs—a direction that, to our knowledge, has not been explored before. As such, there are no directly comparable baselines. This is similar to the situation faced by Grattarola et al. (NeurIPS 2021), who introduced the first GNN-based surrogate and were also unable to compare against prior methods.\n\nBelow, we summarize our responses to the baseline-related concerns raised by the other reviewers.\n\nReviewer BiVF suggested we focus on macro-level surrogates. However, most of these surrogates do not model system dynamics directly, but rather map parameters to summary statistics. While a few exceptions exist (e.g., Fadikar et al., Kieu et al.), they still model aggregate variables and struggle to generalize beyond in-sample dynamics. For instance, an AR(1) time-series model fails to reproduce even simple non-monotonic patterns like predator-prey cycles, significantly underperforming compared to our individual-level approach.\n\nIn response to reviewer McsJ suggestions, we conducted an additional evaluation, training a GNN-only ablation. While this version performs adequately/decently in deterministic settings, it struggles with stochastic dynamics, showing much higher errors than the full GDN, particularly in the predator-prey model. In the Schelling model, where behavior is more deterministic under certain parameterizations, the GNN-only performance was closer to GDN. Reviewer McsJ also suggested that we considered alternative  GNNs (e.g., Grattarola et al., Gravina et al.), but these alternatives  suffer from the same structural limitation of our ablation study: they predict deterministic outcomes and cannot model stochastic transitions. This is supported by the underperformance of the GNN-only ablation. Another proposed baseline is the work from Casert et al., a promising Transformer-based stochastic emulator. However, their method is not applicable to our setting. It assumes fixed-size state spaces, continuous time, and explicit enumeration of all allowed transitions — constraints incompatible with ABMs like Schelling or Predator-Prey, where state spaces are large and dynamic.\n\nWe will include the new experiments, along with the AR(1) baseline and an extended discussion, in the revised paper and supplement.\n\n\nFinally, the reviewer raised a concern regarding the lack of related work, requesting that we “provide the related works.” While we will revise the manuscript to include the expanded discussion presented above, we believe that Section 2 already covers the key contributions of the most relevant prior work in this area. That said, we are fully open to including additional references and discussion if the reviewer can point to specific missing works or directions that merit inclusion.\n\n---\nWe welcome any follow-up questions and would be happy to provide further details or clarifications."}]

Author rebuttal:

Meta-review: This paper proposes a novel framework, which combines graph neural networks and diffusion models to learn a fully differentiable surrogate model directly from ABM-generated simulation data. The reviewers pointed out that the paper is clear and well-written and the proposed architecture is sound. The proposed method also shows strong performance compared with the baselines. 

The reviewers also pointed out some concerns, including missing details, scalability, and a lack of discussion on related works, which are well addressed by the authors during the rebuttal. I encourage authors to include all the rebuttal contents in their camera-ready version.
Decision: Accept (poster)