Paper deep dive
Collective cooperation without individual fidelity in LLM agents
Henrique Ferraz de Arruda, Carlos Gracia LĂĄzaro, Alberto Aleta, Yamir Moreno
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 96%
Last extracted: 7/5/2026, 3:59:43 AM
Summary
The paper investigates the fidelity of Large Language Model (LLM) agents as proxies for human decision-making in social simulations, specifically using a networked Prisoner's Dilemma benchmark. The researchers compared nine open-weight LLMs against human data from a large-scale experiment. The study identifies a 'macro-micro dissociation': while LLM populations can successfully reproduce macro-level features like the temporal evolution and stabilization of cooperation, they fail to capture micro-level nuances, such as individual-level heterogeneity and specific conditional cooperation decision rules. The findings suggest that aggregate agreement in social simulations is insufficient for validating LLMs as faithful human surrogates.
Entities (8)
Relation Signals (4)
LLM Agents â exhibits â Macro-Micro Dissociation
confidence 100% ¡ These findings reveal a macroâmicro dissociation in LLM-based social agents
LLM Agents â testedagainst â Prisoner's Dilemma
confidence 100% ¡ Here we test LLM agents against a direct empirical benchmark: a large-scale networked Prisoner's Dilemma experiment
Lattice Network â usedin â Prisoner's Dilemma
confidence 100% ¡ we fix the interaction structure to a regular lattice network
llama4:16x17b â reproducesmacrofeaturesof â Gracia-LĂĄzaro et al.
confidence 90% ¡ Among the models tested, llama4:16x17b yields the closest aggregate match to the empirical data.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Large language models (LLMs) are increasingly used as agents in simulations of social systems, yet it remains unclear when their behavior can be interpreted as a faithful proxy for human decision-making. Here we test LLM agents against a direct empirical benchmark: a large-scale networked Prisoner's Dilemma experiment with human participants. Using the same interaction protocol, payoff structure, and network topologies, we compare nine open-weight LLMs with the human data. The selected model reproduces several macro-level features of cooperation dynamics, including the early decline and later stabilization of cooperation. This aggregate agreement, however, does not extend uniformly to finer levels of behavior. LLM populations underestimate individual-level heterogeneity and generate conditional cooperation patterns that differ from those observed in humans. Adding a fraction of random agents improves some aspects of micro-level agreement, but does not remove the mismatch in decision rules. These findings reveal a macro--micro dissociation in LLM-based social agents: collective outcomes can appear human-like even when the underlying behavioral distributions and mechanisms are not. They suggest that validating LLM agents as human surrogates requires comparisons across aggregate dynamics, individual heterogeneity, and context-dependent decision rules, rather than outcome-level agreement alone.
Tags
Links
- Source: https://arxiv.org/abs/2606.30454v1
- Canonical: https://arxiv.org/abs/2606.30454v1
Trouble viewing inline? Open PDF directly â
Full Text
104,216 characters extracted from source content.
Expand or collapse full text
Collective cooperation without individual fidelity in LLM agents Henrique Ferraz de Arruda 1,2 , Carlos Gracia L Ěazaro 3 , Alberto Aleta 2,4 , Yamir Moreno 2,4 1 ARAID Foundation, Zaragoza, Spain. 2 Institute for Biocomputation and Physics of Complex Systems, University of Zaragoza, Zaragoza, Spain. 3 Universidad San Jorge (USJ), Zaragoza, Spain. 4 Department of Theoretical Physics, University of Zaragoza, Zaragoza, Spain. Abstract Large language models (LLMs) are increasingly used as agents in simulations of social systems, yet it remains unclear when their behavior can be interpreted as a faithful proxy for human decision-making. Here we test LLM agents against a direct empirical benchmark: a large-scale networked Prisonerâs Dilemma exper- iment with human participants. Using the same interaction protocol, payoff structure, and network topologies, we compare nine open-weight LLMs with the human data. The selected model reproduces several macro-level features of cooperation dynamics, including the early decline and later stabilization of coop- eration. This aggregate agreement, however, does not extend uniformly to finer levels of behavior. LLM populations underestimate individual-level heterogene- ity and generate conditional cooperation patterns that differ from those observed in humans. Adding a fraction of random agents improves some aspects of micro- level agreement, but does not remove the mismatch in decision rules. These findings reveal a macroâmicro dissociation in LLM-based social agents: collective outcomes can appear human-like even when the underlying behavioral distri- butions and mechanisms are not. They suggest that validating LLM agents as human surrogates requires comparisons across aggregate dynamics, individual heterogeneity, and context-dependent decision rules, rather than outcome-level agreement alone. Keywords: Large Language Models; Networked Social Dilemmas; Prisonerâs Dilemma; Conditional Cooperation; Behavioral Heterogeneity 1 arXiv:2606.30454v1 [physics.soc-ph] 29 Jun 2026 1 Introduction Large language models (LLMs) are increasingly used as proxies for human agents in simulations, evaluation settings, and collective decision-making environments [1â3]. This development is part of a broader shift in which language models are no longer treated only as systems that generate text, but also as tools to build agents whose behavior can be elicited, measured, and compared with that of humans. Recent work has extended this approach to experimental and agent-based settings, where LLM agents are used to simulate human-like behavior in social and strategic contexts [1, 4]. The appeal of this approach is clear. If LLM agents could reproduce human behavior at scale, they would provide a powerful tool for computational social sci- ence, behavioral experimentation, and the construction of social simulations. Yet this promise also raises a basic validation problem. Agreement with human outputs does not necessarily imply that a model reproduces the behavioral patterns or decision processes that generated those outputs. Several recent studies point in this direc- tion. LLMs may reconstruct users through compressed behavioral representations that introduce systematic biases in social simulations [5]; they may align with human eval- uative judgments while relying more heavily on lexical associations and statistical priors than on contextual reasoning [6]; and LLM-generated text can display statisti- cal signatures, such as higher structural regularity and compressibility, that differ from human language production [7]. These observations suggest that apparent behavioral realism may be shallow if it is assessed only at the level of visible outcomes. This issue is especially important in social simulations and game-theoretic environ- ments, where LLM agents are increasingly used as substitutes for human participants. LLMs can participate effectively in repeated strategic interactions, particularly in self-interested settings such as the iterated Prisonerâs Dilemma family, while show- ing weaker performance in coordination problems [8]. Other work has shown that LLM agents can develop collective conventions and shared behavioral patterns with- out explicit programming [9]. At the same time, these agents remain distinguishable from humans in interactive settings: in online debates, for example, LLM agents can remain on topic and blend into conversations, while still being perceived as less con- vincing and less confident than human interlocutors [4]. Related studies show that LLM outputs can approach human benchmarks in some evaluative tasks and can be systematically shaped through prompting and fine-tuning [10, 11]. A further complication is that current LLMs are not primarily optimized to repro- duce human behavior in experiments. They are generally trained and aligned to function as helpful assistants in human-oriented tasks [12]. This distinction matters because a model that is useful, compliant, or socially desirable need not be a faithful model of human decision-making. Alignment procedures may also compress heteroge- neous preferences into narrower sets of high-probability behaviors, as shown by Slocum et al. [13], with related empirical evidence of mode collapse in aligned language mod- els [14]. These concerns are consistent with recent commentaries emphasizing the lack of empirical validation against real human behavioral data in LLM-based agent mod- els [15], and with calls to calibrate scientific claims from LLM social simulations to the strength of robustness checks and empirical benchmarks [16]. 2 Here, we evaluate whether LLM agents can function as behavioral surrogates for humans in repeated networked social dilemmas. We use as a benchmark the large- scale Prisonerâs Dilemma experiment of Gracia-L Ěazaro et al. [17], in which human cooperation converged to similar levels across network topologies, challenging the the- oretical expectation that heterogeneous networks should promote cooperation. Using the same network structures and interaction protocol, we compare nine open-weight LLMs operating as autonomous agents in repeated Prisonerâs Dilemma games. Our aim is not to ask whether LLMs play the game better than humans, but whether they reproduce the empirical temporal dynamics, individual heterogeneity, and conditional behavioral patterns observed in the experiment. This design allows us to test behavioral similarity at several levels. At the macro level, we ask whether LLM agents reproduce the aggregate evolution of cooperation. At the micro level, we ask whether they reproduce the distribution of individual cooperation propensities. At the decision-rule level, we ask whether they reproduce conditional cooperation, namely, the way in which participants adjust their actions to the previous behavior of their neighbors. We find that the selected LLM reproduces several macro-level regularities, including the early decay and later stabilization of cooperation. However, this aggregate agreement does not fully extend to micro-level signatures: simulated populations underestimate behavioral heterogeneity and gener- ate conditional cooperation rules that differ quantitatively from the empirical data. These results reveal a macroâmicro dissociation in LLM-based social agents, with direct implications for the validation of machine behavior and for the use of LLM agents as human surrogates in experimental social science. 2 Results We present the results in three steps. We first describe the empirical benchmark and the simulation framework. We then assess whether LLM agents reproduce aggre- gate cooperation dynamics, before turning to individual heterogeneity and conditional cooperation. This ordering follows the central validation question of the study: whether agreement at the collective level extends to finer behavioral levels. 2.1 Empirical benchmark and simulation framework To evaluate whether LLM agents reproduce human cooperation dynamics, we simu- lated a repeated multi-agent Prisonerâs Dilemma on network structures derived from a human behavioral experiment. The simulation framework consists of an orchestra- tion layer that manages autonomous agents, each controlled by an LLM. This work builds on the experimental design of Gracia-L Ěazaro et al. [17], who examined human cooperation in networked Prisonerâs Dilemma games. The original study found that human cooperation converged to similar steady-state levels across network topolo- gies, challenging theoretical expectations that heterogeneous networks should promote cooperation. Using the same network structures and temporal dynamics as in the orig- inal study, we test whether LLM-based agents reproduce the empirical cooperation patterns or instead follow behavior closer to classical game-theoretic predictions. 3 The simulation environment is a repeated multi-agent game in which each net- work node is occupied by an LLM-driven agent. The simulation proceeds in discrete rounds that mirror the original experiment. In each round, after the introductory one, agents receive a prompt summarizing the previous round, including their neighborsâ choices and normalized payoffs. Agents then return their choice for the new one which effectively could be to cooperate or to defect. The underlying incentive structure is a Prisonerâs Dilemma where payoffs are computed pairwise. The total payoff for an agent is the sum of these pairwise interactions across all their network connections. The specific values used to drive the agentsâ goal of maximizing âECUsâ (Experimen- tal Currency Units) are detailed in Table 1 and follow Grujic et al (2010) [18] as in the original experiment. In this specific configuration, the rewards are designed such that cooperation is theoretically expected to reach a high level [17, 18]. Implementation details, model descriptions, and statistical procedures are provided in the Methods. Table 1: Pairwise payoff matrix (ECUs per neighbor). Your Choice / Neighbor ChoiceCooperate Defect Cooperate70 Defect100 2.2 Macro-level cooperation dynamics We first ask whether different LLMs produce comparable aggregate cooperation dynamics when placed in the same experimental environment. For this initial bench- mark, we fix the interaction structure to a regular lattice network (625 nodes, average degree 4). Figure 1(a) shows the resulting cooperation trajectories across all models, while Fig. 1(b) shows their average cooperation levels. The tested LLMs span several model families, parameter scales, and training paradigms. Despite identical incentives and network structure, they generate markedly different cooperation regimes. Notably, qwen3:32b produces the highest cooperation levels across the simulation horizon. Ablated variants derived from the same base checkpoints also display systematically altered trajectories relative to their originals, indicating that collective outcomes are sensitive to model choice and alignment. Among the models tested, llama4:16x17b yields the closest aggregate match to the empirical data. We therefore use this model as the primary agent in the remaining analyses, while treating the initial benchmark as evidence that LLMs should not be regarded as interchangeable behavioral surrogates. With this LLM fixed for the remaining analyses, we next evaluate the two network structures considered in [17]: a regular lattice and a heterogeneous network (604 nodes, degrees ranging from k = 2 to k = 16). Fig. 2(a) shows the temporal evolution of cooperation for one realization of the dynamics. As in the empirical experiment, both structures yield similar cooperation levels that overlap closely over time. Compared to the empirical data, the simulated trajectories align more closely in later iterations, though they remain noisier overall. 4 Empirical 1 2 3 6 5 8 4 7 9 a)b) Fig. 1: Comparison of cooperation outcomes across LLM models. The numbers shown on the plot correspond to the following models: (1) BlackHillsIn- foSec llama-3.1:8b-abliterated, (2) deepseek-llm:67b, (3) deepseek-llm:7b, (4) krith meta-llama-3.1:70b-instruct-abliterated IQ4 XS, (5) llama3.1:70b, (6) llama3.1:8b, (7) llama4:16x17b, (8) smolLM2.1:7b, and (9) qwen3:32b. Experiments performed for the Lattice network. Panel (a) shows the mean cooperation rate per round (empirical) or per iteration (simulation), and the error bars represent the standard error. Panel (b) shows the average values of the simulation curves, with error bars representing the standard deviation. The dashed horizontal line represents the empirical average, and the shaded area represents one standard deviation around it. To perform a thorough comparison, we also test the control cases, which con- sist of simulations on a time-varying network structure. Specifically, connections are randomly rewired at each iteration while preserving the original degree sequence. To match the original study, we do not generate a new rewiring sequence. Instead, we use the exact sequence of network connections reported in [18], which was originally gener- ated at random. Henceforth, we refer to the fixed-network simulations as the networked condition, and to the time-varying simulations as the control condition. Figure 2(f) shows the resulting cooperation trajectories over time. Notably, the heterogeneous control case more closely matches the empirical data, whereas in the lattice control case, agents tend to cooperate less. As in the networked simulations, the LLM-based trajectories remain noisier overall. Upon closer inspection of the cooperation histograms, we observe that the empirical experiments yield substantially broader distributions than those of the simulations under the networked condition. Specifically, the empirical lattice case (Fig. 2b) has mean cooperation Îź = 0.368 with standard deviation Ď = 0.170, and the empirical heterogeneous case (Fig. 2c) has Îź = 0.405 and Ď = 0.158. By contrast, the simulated lattice network (Fig. 2d) produces a more concentrated distribution (Îź = 0.291, Ď = 0.069), whereas the simulated heterogeneous network (Fig. 2e) yields Îź = 0.314 and Ď = 0.148. While the lattice simulation more closely reproduces the empirical mean, the heterogeneous simulation better captures the observed dispersion and therefore provides a closer overall match. 5 1020304050 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation a) Emp.: Lattice Emp.: Heterogeneous Sim.: Lattice Sim.: Heterogeneous 1020304050 Iteration 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation f) 0 50 b) Emp.: Lattice c) Emp.: Heterogeneous 0 100 d) Sim.: Lattice e) Sim.: Heterogeneous 0 50 g) Emp.: Lattice h) Emp.: Heterogeneous 0.00.51.0 0 200 i) Sim.: Lattice 0.00.51.0 j) Sim.: Heterogeneous Cooperative Rounds / Total Rounds Number of Players Fig. 2: Evolution and distribution of cooperation. Panels (a) and (f) show the temporal evolution of the mean cooperation rate per iteration for the networked and control experiments, respectively. Dashed lines represent empirical results, solid lines represent simulation results, and error bars indicate the standard error. Panels (b)-(e) and (g)-(j) display the histograms of individual cooperation propensities (cooperative rounds / total rounds). Specifically, Panels (a)-(e) show the results for the networked experiment, and Panels (f)-(j) show the results for the control case. To test the sensitivity of our results to agent variability, we repeated the lattice simulation with LLM temperatures sampled uniformly at random. The resulting dis- tributions were similar to those reported in the main text, showing that our findings are not affected by temperature variation (see Supplementary Material, Fig. S13). We also test three different âpersonasâ for the LLMs. However, since the results are con- sistent across personas, we omitted them from the paper (see Supplementary Material, Section S3.2.1, for more details). Furthermore, to evaluate whether the reduced disper- sion observed in the simulated cooperation distributions is specific to llama4 16x17b, we repeated the analysis using qwen3:32b, the second-best-performing model accord- ing to Fig. 1. As shown in Supplementary Fig. S14, qwen3:32b reproduces the same qualitative pattern, in which the simulated cooperation distributions remain narrower than the corresponding empirical distributions. These results suggest that the tendency of LLM-based agents to generate less variable outcomes than human participants is not restricted to a single model. A similar pattern is observed in the control condition. The empirical lattice control (Fig. 2g) exhibits Îź = 0.292 and Ď = 0.209, and the empirical heterogeneous control (Fig. 2h) shows Îź = 0.322 and Ď = 0.209, again indicating wide distributions. In contrast, the simulated lattice control (Fig. 2i) remains narrowly concentrated (Îź = 0.238, Ď = 0.047), while the simulated heterogeneous control (Fig. 2j) produces higher cooperation (Îź = 0.361) but still substantially reduced variability (Ď = 0.089). Thus, although the heterogeneous control simulation better matches the empirical mean 6 cooperation level, both simulated control cases underestimate the dispersion observed in human behavior. To evaluate the fidelity of the simulated temporal dynamics, we used metrics to compare simulated and empirical cooperation trajectories. In addition to measuring point-by-point error with Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE), we used an autoregressive model (AR(1)) to evaluate the persistence of cooperation. The difference in persistence coefficients (Diff Ď ) indicates whether sim- ulated and empirical trajectories exhibit similar temporal persistence (i.e., a form of aggregate âmemoryâ in the cooperation series). The difference in implied mean incre- ments (Diff Îź ) identifies systematic biases in the average rate of change. Finally, to measure the synchrony of cooperation trends, we calculated the Pearson correlation coefficient for the raw levels (r) and the first differences (r â ), ensuring that the simu- lated fluctuations were synchronized with the empirical observations. To ensure that these results were not due to random chance, we established statistical significance using a Monte Carlo null model that preserves marginal distributions while removing temporal structure. Table S1 shows that the simulations reproduce the overall shape of the empirical cooperation trajectories, but still deviate systematically in temporal structure. Point- wise errors (RMSE and MAE) are lowest in the heterogeneous control condition and highest in the heterogeneous networked condition. The AR(1) persistence compari- son shows that lattice simulations differ strongly from humans in decision âmemory,â with large and significant Diff Ď values (0.217 and 0.464), whereas the heterogeneous networked case shows the closest match in persistence (Diff Ď = 0.071). Differ- ences in implied long-run means remain small across all conditions (Diff Îź ⤠0.004), suggesting limited bias in average cooperation levels. Finally, correlations in raw cooperation levels are consistently moderate and significant, but correlations in first differences are weaker, indicating that the simulations capture broad trends better than round-to-round fluctuations. As an exploratory test of whether non-strategic behavior contributes to human- like variability, we introduced a fraction of random agents. This choice is motivated by qualitative evidence from the original experiment suggesting that some participants may have played without a stable strategy, or even at random [17]. To reproduce this situation, we assign a fraction Ď of nodes to act as random agents. These agents do not follow the LLM policy, but choose between cooperate and defect uniformly at random in every round. The remaining fraction (1â Ď) of nodes is controlled by the LLM- driven decision rule. Based on a systematic sweep over values of Ď (Supplementary Information, Table S3), we use Ď = 0.2 as an illustrative hybrid condition. This value balances improved temporal alignment with the empirical cooperation trajectories, particularly in terms of fluctuation synchrony (r â = 0.469), with preservation of the overall trajectory shape (r = 0.690), while avoiding the stronger distortions observed at higher random-agent fractions. Alternative trajectories for other values of Ď are shown in Supplementary Information, Fig. S15. Figure 3 shows that introducing a fraction of random agents (Ď = 0.2) preserves the overall qualitative shape of the empirical trajectories, with cooperation rapidly decay- ing in early rounds and subsequently stabilizing. However, the inclusion of random 7 ExperimentRMSEMAEDiff Ď Diff Îź r â lattice0.098**0.083**0.217***0.0020.612***0.276* lattice control0.089***0.074*0.464***0.002***0.583***0.214 heterogeneous0.107**0.0960.0710.0020.622***0.151 heterogeneous control0.079***0.0510.3590.004*0.440***0.362*** Table 2: Trajectory errors, autoregressive comparison, and correlation. RMSE and MAE measure point-by-point deviations between empirical and simulated cooperation rates. Diff Ď denotes the absolute difference between the AR(1) persistence coefficients, and Diff Îź denotes the absolute difference between the implied long-run means of the simulated and empirical series. r denotes the Pearson correlation coefficient between empirical and simulated trajectories. r â represents the Pearson correlation coefficient using the first differences. For RMSE, MAE, Diff Ď , and Diff Îź , statistical sig- nificance is assessed using the Monte Carlo null model described in the text. For r and r â , significance refers to the null hypothesis of zero correlation (H 0 : r = 0). â denotes p < 0.01, â denotes p < 0.05, â denotes p < 0.1, and no symbol indicates that the result is not statistically significant. 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation a)b) 1020304050 Iteration 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation c) 1020304050 Iteration d) Empirical Hybrid Simulation (80% LLM + 20% Rand.) Full LLM Simulation LLM Subset (within Hybrid) Fig. 3: Evolution of cooperation with a fraction of random agents (Ď = 0.2). Panels show cooperation over iterations for (a) lattice, (b) lattice control, (c) heterogeneous, and (d) heterogeneous control network structures. Solid lines denote the mean cooperation rate per round (empirical) or per iteration (simulation), and error bars represent the standard error. agents alters the quantitative agreement with human data in a condition-dependent manner. In the lattice network, adding random agents increases the correlation with 8 the empirical trajectory and substantially improves synchrony in short-term fluctua- tions (r â ), although pointwise errors increase slightly. A similar pattern holds for the heterogeneous network, where the hybrid case yields stronger correlation (r = 0.743 vs. r = 0.622) and higher agreement in first differences (r â = 0.343 vs. 0.151), while maintaining comparable error levels (RMSE = 0.109 vs. 0.107). In the control scenarios, the effect of random agents is more mixed. In the lat- tice control case, the hybrid simulation slightly reduces RMSE (0.087 vs. 0.089) and markedly lowers MAE (0.046 vs. 0.074), but decreases correlation with the empirical trajectory (r = 0.509 vs. r = 0.583). In the heterogeneous control case, the hybrid configuration yields the lowest RMSE across all conditions (0.070), but also produces a weaker correlation (r = 0.325), compared to r = 0.440 in the scenario without random agents. See Supplementary Material, Table S4, for more details on all the trajectory errors. Together, these results indicate that random agents can improve alignment with human temporal fluctuations in the networked condition, while in the control condition, they primarily reduce absolute error without consistently improving temporal synchrony. 2.3 Micro-level heterogeneity and decision rules We next ask whether aggregate agreement is accompanied by individual-level fidelity. To do so, we compare the distribution of cooperation propensities across agents, rather than only the population mean. This distinction is important because two pop- ulations can display similar average cooperation while differing substantially in the diversity of individual behaviors. We quantify the distance between empirical and sim- ulated propensity distributions using the first-order Wasserstein distance, W 1 (E,S). Lower values indicate closer agreement in behavioral heterogeneity, whereas larger values indicate that the simulated population differs from the empirical one in how cooperation is distributed across individuals. Table 3 summarizes the results across experimental scenarios. Table 3 shows that the simulations match empirical individual-level heterogeneity with varying accuracy across experimental scenarios. The smallest distance is obtained in the lattice network with mixed populations, where W 1 (E,S) = 0.081, indicating the closest match to the empirical distribution of cooperation propensities. In contrast, the lattice control condition with a fully LLM-driven population yields the largest discrep- ancy, W 1 (E,S) = 0.136. Introducing random agents generally reduces the distance relative to the corresponding fully LLM-driven case, with the strongest improvement observed in the lattice control condition (from 0.136 to 0.095). In the heterogeneous network, however, the mixed-population case slightly worsens the match under the networked condition (from 0.091 to 0.095). For each network condition (lattice, lattice control, heterogeneous, heterogeneous control), we quantified behavioral heterogeneity in the empirical data as the distribu- tion of agent-level cooperation propensities (fraction of cooperative actions per agent). We then generated a parametric null model in which each agent follows an independent Bernoulli process with success probability equal to the network-level empirical mean cooperation rate, while preserving each agentâs number of rounds. Across (B=10,000) 9 ExperimentW 1 (E, S)W 1 (E, Null) lattice (100% LLM + 0% Rand.)0.1010.085 [95% CI: 0.082, 0.089] lattice (80% LLM + 20% Rand.)0.0810.085 [95% CI: 0.082, 0.089] lattice control (100% LLM + 0% Rand.)0.1360.125 [95% CI: 0.123, 0.128] lattice control (80% LLM + 20% Rand.)0.0950.125 [95% CI: 0.123, 0.128] heterogeneous (100% LLM + 0% Rand.)0.0910.074 [95% CI: 0.071, 0.077] heterogeneous (80% LLM + 20% Rand.)0.0950.074 [95% CI: 0.071, 0.077] heterogeneous control (100% LLM + 0% Rand.)0.1090.126 [95% CI: 0.123, 0.129] heterogeneous control (80% LLM + 20% Rand.)0.0930.126 [95% CI: 0.123, 0.129] Table 3: Wasserstein distance between empirical and simulated dis- tributions of individual cooperation propensities across experimental scenarios. Lower values indicate better reproduction of behavioral heterogeneity. The null model assumes that each agent follows an independent, identically dis- tributed, Bernoulli process with the networkâs empirical mean cooperation rate, all while maintaining the same number of rounds per agent. null realizations, we computed W 1 (E, Null) between empirical and null propensity dis- tributions. We report the mean null distance and a 95% confidence interval of this Monte Carlo distribution, providing a baseline for whether observed heterogeneity departs from what is expected under homogeneous random behavior. To contextualize these values, we compare the results of Table 3 with the null-model baseline. In the lattice network, the null model yields W 1 (E, Null) = 0.085 (95% CI: [0.082, 0.089]), which is close to the best-performing simulation (W 1 (E,S) = 0.081), indicating that the empirical heterogeneity in this condition can largely be explained by variation consistent with independent behavior around the population mean. In the lattice control condition W 1 (E, Null) = 0.125 (95% CI: [0.123, 0.128]), which is not comparable to the full-LLM case (W 1 (E,S) = 0.136), while the mixed simula- tion (W 1 (E,S) = 0.095) falls below the null baseline. In the heterogeneous network, the null distance is smaller (W 1 (E, Null) = 0.074, 95% CI: [0.071, 0.077]) than both simulation cases (W 1 (E,S) = 0.091 and 0.095), suggesting that the simulated populations exhibit less empirical-like heterogeneity than expected under the null benchmark. Finally, in the heterogeneous control condition, the null distance is again high (W 1 (E, Null) = 0.126, 95% CI: [0.123, 0.129]), whereas both simulation scenarios yield lower distances (W 1 (E,S) = 0.109 and 0.093), indicating improved agreement with empirical heterogeneity compared to the null expectation. We then examine a stricter behavioral signature: conditional cooperation. This measures the probability that an agent cooperates as a function of the fraction of cooperating neighbors in the previous round, and according to whether the agent cooperated or defected in the previous round. Conditional cooperation is important because it probes the local decision rule, rather than only the resulting level of coop- eration. As in the empirical study, we first analyze the networked condition, where neighborhoods remain fixed. In the human experiment, this relationship is well approx- imated by a linear trend with distinct responses following cooperation and defection. 10 Across most simulated scenarios, we observe a qualitative dependence on neighbor- hood cooperation, but the responses are quantitatively different and often steeper than in the human data. The âafter defectionâ response deviates most strongly from the empirical pattern in the lattice network, exhibiting an inverted slope in both the full-LLM and hybrid simulations. Moreover, unlike the empirical data, the sim- ulated conditional-cooperation points do not consistently follow an approximately linear relationship, making the linear-fit comparison used in the original study less informative here. We therefore report the full curves in the Supplementary Material, Section S3.2.2, highlighting an additional discrepancy between LLM-based agents and empirical behavioral patterns. To extend this analysis, we examine conditional cooperation across network struc- tures and experimental conditions. Figure 4 shows conditional cooperation across experimental scenarios for a single execution of the dynamics. In each case, the curves represent the estimated conditional probabilities of cooperation as a function of neigh- borhood cooperation in the previous round. We repeated this conditional-cooperation analysis over 20 independent executions of the dynamics and obtained consistent curves with narrow confidence intervals, indicating that a single execution is suffi- cient to capture the qualitative conditional-cooperation trends of the simulations (for details see Supplementary Material, Figs. S22 and S23). The quantitative deviations reported in Table 4 are consistent with the qualita- tive patterns shown in Fig. 4. Across most scenarios, errors are larger in the z = 1 case (after cooperation) than in the z = 0 case (after defection), reflecting stronger discrepancies in how simulated agents condition their behavior on neighborhood coop- eration following a cooperative move. The largest deviations are observed in the lattice network under the networked condition, where the full-LLM population yields RMSE z=1 = 0.397 and MAE z=1 = 0.388, matching the visibly distorted and non- linear response in panel A. Introducing random agents reduces these discrepancies in the lattice networked case (RMSE z=1 decreases to 0.255), although a substantial mismatch remains. In contrast, the lattice control condition shows comparatively small error after defection (RMSE z=0 = 0.046 in the full-LLM case), consistent with the flatter con- ditional patterns visible under rewiring, while errors after cooperation remain large (RMSE z=1 = 0.341). For the heterogeneous network, the mixed population reduces deviations relative to the full-LLM case, particularly in the z = 1 regime (RMSE z=1 decreases from 0.355 to 0.232), which aligns with the closer agreement between the empirical and simulated curves observed in the corresponding panels. Overall, these results confirm that while simulations capture the qualitative dependence of cooperation on local neighborhood behavior, the inferred conditional rules remain quantitatively distinct from the empirical patterns, particularly following cooperative actions. 3 Discussion This study tested whether autonomous LLM agents can serve as behavioral surro- gates for humans in a repeated networked social dilemma. The answer is mixed in a 11 After CooperationAfter Defection a)b) c)d) e)f) g)h) Fig. 4: Conditional cooperation across network structures and experimen- tal conditions. Each panel shows the estimated probability of cooperation as a function of the fraction of cooperating neighbors in the previous round, computed from a single execution of the dynamics. Panels (a) and (b) correspond to the lattice network under the networked condition, (c) and (d) to the lattice network under the control condition, (e) and (f) to the heterogeneous network under the networked con- dition, and (g) and (h) to the heterogeneous network under the control condition. In each case, empirical conditional cooperation patterns are compared with simulations under different agent compositions, including fully LLM-driven and mixed populations with random agents. Error bars indicate 95% confidence intervals. revealing way. The selected model, llama4:16x17b, reproduces several macro-level fea- tures of the human experiment, including the early decline and later stabilization of cooperation observed by Gracia-L Ěazaro et al. [17], as well as the absence of network 12 ExperimentRMSE z=0 MAE z=0 RMSE z=1 MAE z=1 lattice (100% LLM + 0% Rand.)0.237**0.217***0.397***0.388* lattice (80% LLM + 20% Rand.)0.1800.165***0.255***0.246* lattice control (100% LLM + 0% Rand.)0.046***0.038***0.341***0.337*** lattice control (80% LLM + 20% Rand.)0.090*0.086*0.221***0.211*** heterogeneous (100% LLM + 0% Rand.)0.177***0.1540.355***0.320*** heterogeneous (80% LLM + 20% Rand.)0.161***0.1400.232**0.214 heterogeneous control (100% LLM + 0% Rand.)0.107**0.091***0.2170.196 heterogeneous control (80% LLM + 20% Rand.)0.1900.1240.1730.161** Table 4: Deviation between empirical and simulated conditional cooperation rules across experimental scenarios. Errors are reported separately for cases in which the previous action was defection (z = 0) or cooperation (z = 1). â denotes p < 0.01, â denotes p < 0.05, â denotes p < 0.1, and no symbol indicates that the result is not statistically significant. reciprocity reported in the original study. Yet this aggregate agreement does not imply individual-level fidelity. The simulated populations underestimate important aspects of human heterogeneity and generate conditional cooperation patterns that differ from the empirical decision rules. This pattern is the central result of the paper. LLM-based agents can reproduce collective cooperation dynamics without fully reproducing the distribution of indi- vidual behaviors or the context-dependent mechanisms through which humans adapt their decisions. In other words, the model captures part of the phenomenology of the experiment, but not all of its behavioral structure. This macroâmicro dissociation mat- ters because many proposed uses of LLM agents in social simulation rely implicitly on the assumption that matching aggregate outcomes is sufficient evidence of behavioral realism. Our model benchmark also shows that LLMs cannot be treated as interchangeable agents. Different model families produce markedly different cooperation regimes, and alignment-perturbed variants derived from related checkpoints can generate system- atically different collective dynamics. Model choice is therefore a substantive modeling assumption, not a technical detail. This sensitivity reinforces the need for empirical benchmarking before drawing conclusions from LLM-based social simulations [15]. The micro-level results are especially informative. Human participants display broad distributions of cooperation propensities, whereas the simulated populations often produce narrower or differently shaped distributions. Conditional cooperation provides an even stronger test. Human behavior in the original experiment is charac- terized by structured responses to the local social environment, with distinct patterns after cooperation and after defection. The LLM agents reproduce some qualitative dependence on neighborhood cooperation, but the inferred rules differ quantitatively, and in some cases qualitatively, from the empirical patterns. Thus, the model does not simply implement a noisy version of the human rule. It appears to reach partially similar aggregate outcomes through different local response functions. 13 The hybrid simulations with random agents help clarify this point. Introducing a fraction of random decision-makers improves agreement with the empirical data in some dimensions, particularly fluctuation synchrony in the networked condition and distributional agreement in selected cases. This should not be interpreted as a tuned correction that makes the model human-like. Rather, it suggests that stochasticity, inconsistent behavior, or non-strategic play may be necessary ingredients in models of human cooperation. Real participants are not homogeneous optimizers, and a faithful behavioral surrogate may need to represent this diversity explicitly. At the same time, the hybrid model does not remove the mismatch in conditional cooperation, indicating that noise alone is not sufficient to recover human decision mechanisms. These results have broader implications for the validation of machine behavior. A model may be useful for reproducing an aggregate trajectory while still being inadequate for causal explanation or intervention design. This distinction is particu- larly important for digital twins, policy simulations, and agent-based models that use LLMs as substitutes for human actors. If the aim is only to forecast a limited aggre- gate outcome under conditions close to the benchmark, macro-level agreement may be informative. If the aim is to understand how individuals respond to incentives, social context, or interventions, then micro-level heterogeneity and decision-level rules become essential. Several limitations should be kept in mind. The empirical benchmark comes from a specific experiment, participant pool, and institutional context, and the same compar- ison should be repeated across other social dilemmas and populations. The simulations also rely on a particular prompt structure and on open-weight models available at the time of the study. Although persona variants and repeated executions did not qualita- tively change the results, broader prompt and model robustness remain an important direction for future work. More broadly, our findings suggest that LLM simulations cannot currently replace human experimentation, which remains necessary to establish the empirical ground truth used for model validation. However, when benchmarked against human data, LLM agents may still be valuable for scaling simulations and exploring conditions that are difficult to study experimentally. Finally, current LLMs are typically optimized to follow human instructions and preferences, not to repro- duce the full diversity of human behavior in experiments [12]. This may contribute to the reduced behavioral diversity observed here, consistent with forms of model homogenization and mode collapse reported in recent work [13, 14]. More generally, our findings support the view that claims based on LLM social simulations should be calibrated to the strength of empirical validation and robustness checks [16]. Overall, llama4:16x17b provides a promising but incomplete surrogate for human behavior in networked social dilemmas. It reproduces robust macro-level cooperation dynamics, but not the full heterogeneity and conditional structure of human strategic decision-making. The main lesson is therefore not simply that LLM agents succeed or fail as human surrogates, but that surrogacy is level-dependent. Validating machine behavior requires asking not only whether artificial agents produce similar outcomes, but also whether they do so through similar behavioral distributions and decision mechanisms. 14 4 Methods 4.1 Simulation framework All LLM-based agents were implemented using the CrewAI orchestration framework 1 , which handles prompt construction, model querying, response parsing, and payoff updates at each iteration. Model inference was performed locally using Ollama 2 as the serving backend, enabling self-contained execution of all agents without external API dependencies. Each agent was mapped to a node in the network and interacted only with the information in its local prompt, which contained its neighborsâ previous-round actions and the corresponding normalized payoffs. Models were queried independently for each agent at each round, ensuring decentralized decision-making consistent with the experimental protocol. We simulate a repeated Prisonerâs Dilemma game on empirical social networks (using the same network topologies and temporal dynamics as in Gracia-L Ěazaro et al. [17]). Each node (agent) chooses between GREEN and BROWN in each round, corresponding to cooperation and defection, respectively. We use color labels rather than the canonical game actions in the agent-facing prompts to reduce the likelihood that the LLM recognizes the task from familiar Prisonerâs Dilemma terminology and responds on that basis rather than from the local interaction history and payoffs. The game proceeds in discrete rounds. In each round, each agent observes its neighborsâ previous choices and (degree-)normalized payoffs, then chooses an action for the cur- rent round. Payoffs were computed pairwise as follows: if both agents choose GREEN, each receives 7 ECUs; if one chooses BROWN and the neighbor chooses GREEN, the BROWN agent receives 10 ECUs, and the GREEN agent receives 0 ECUs; if both choose BROWN, both receive 0 ECUs. An agentâs total payoff in a round is the sum of these pairwise payoffs with all neighbors. In the last experiments, a fraction Ď of nodes were designated random (i.e., choosing actions uniformly), while all other nodes used the LLM-driven policy. We also tested conditions with optional backstory cues in the agent prompts (e.g., labeling an agent as a student or specifying a gender, or both); details of these vari- ants are omitted here. Full prompt templates, backstory variants that define the persona, and implementation details, including parsing routines, are provided in the Supplementary Information, Section S1. 4.2 LLM models and framework We evaluate nine open-weight large language models spanning multiple model fam- ilies, parameter scales, and training paradigms, including variants of LLaMA [19], DeepSeek [20], Qwen [21], and SmolLM [22]. All models were selected to represent a broad spectrum of open-source LLMs commonly used in agent-based simulations. The evaluated models belong to four open-weight families. LLaMA models are decoder-only transformer language models developed by Meta, available in multiple generations and sizes. In this study, we use both standard instruction-tuned variants (i.e., llama3.1:8b, 1 https://w.crewai.com/, accessed March 2026 2 https://ollama.com/, accessed March 2026 15 llama3.1:70b, and llama4:16x17b) and alignment-modified derivatives (i.e., Black- HillsInfoSec llama-3.1:8b-abliterated and meta-llama-3.1:70b-instruct-abliterated IQ4 XS), allowing us to assess whether differences in model scaling and instruction tuning affect emergent cooperative behavior in multi-agent simulations. In particular, we com- pare standard instruction-tuned LLaMA models with modified âabliteratedâ variants of the same base architecture, which serve as alignment-perturbed counterparts to eval- uate the sensitivity of the results to behavioral steering. SmolLM models are compact open-weight language models designed for efficient inference and deployment. Here we test smolLM2.1:7b. DeepSeek models are open-weight language models released in multiple parameter scales and variants, and we tested two different sizes: deepseek- llm:7b and deepseek-llm:67b. Qwen models are multilingual instruction-tuned large language models, with recent variants (e.g., Qwen3) showing strong performance on reasoning and general language understanding benchmarks, here we use qwen3:32b. 4.3 Quantitative analysis We evaluate whether the simulations reproduce the empirical experiments at two complementary levels: (i)macro dynamics (i.e., the time evolution of the fraction of cooperators in the network), and (i) micro behavior (i.e., heterogeneity across indi- viduals and context-dependent decision rules). Full definitions, robustness checks, and statistical testing procedures are reported in the Supplementary Information (Section S2). 4.3.1 Macro-level agreement For each condition and network, we summarize behavior at iteration t by the cooperation rate X D (t) = 1 N N X i=1 1 a D i (t) = C ,(1) where D â E,S denotes empirical (E) or simulated (S) data, X D (t) is the frac- tion of cooperative agents, and C means cooperation. We then compare empirical and simulated series using two complementary notions of similarity. First, we com- pute pointwise trajectory errors (RMSE and MAE). These measures quantify how close simulations are in level at each iteration, in which smaller values indicate closer agreement, with RMSE penalizing large deviations more. Second, we compute Pearson correlation between the two trajectories, both on the raw series (r) and on first differ- ences (r â ). Correlation captures whether simulations reproduce the temporal pattern, where co-movement of increases and decreases, even when small level shifts persist. 4.3.2 Macro-level temporal dependence Trajectory similarity in level or correlation does not guarantee that the simulated process evolves over time in the same way as the empirical one. To compare tem- poral dependence, we fit a first-order autoregressive model (AR(1)) [23] to the first-differenced series âX D (t) (stationarity check [24] and the differenced specifi- cation are detailed in the Supplementary Information, Section S2.2.4). For a series 16 âX D c,g,r (t), the AR(1) model is given by âX D c,g,r (t) = Îş + Ď âX D c,g,r (tâ 1) + Îľ(t),(2) where Îş is a constant, Ď captures persistence, and Îľ(t) is a white-noise error term. We estimate the model separately for the empirical series X E c,g,r (t) and the simulated series X S c,g,r (t). The coefficient Ď summarizes how strongly changes in cooperation propagate from one iteration to the next: values closer to zero indicate rapid adjustment, whereas larger magnitudes indicate more persistent (or anti-persistent) dynamics. We quantify dynamic similarity using the absolute difference|Ď S âĎ E |, and we additionally compare the implied mean increment of the differenced process, which captures systematic drift in cooperation over time. 4.3.3 Significance via null models Because the absolute magnitudes of these metrics depend on the scale and length of the series, we assess whether the observed agreement exceeds what would be expected by chance. We construct Monte Carlo null time series that preserve the empirical marginal distribution of cooperation rates while destroying temporal dependence, and we compute two-tailed p-values by comparing the observed metric to its null distribu- tion (see Supplementary Information, Section S2.2.5). For correlations, we test against the standard null hypothesis of zero correlation. 4.3.4 Micro-level mechanisms We quantify individual behavioral heterogeneity by computing each agentâs coop- eration propensity over the full time horizon T . For dataset D, this is defined as P D i = 1 T T X t=1 y D i (t),(3) where y D i (t) = 1 a D i (t) = C â 0, 1 indicates whether agent i cooperated at round t. To test whether simulations reproduce individual-level diversity, we compare the empirical and simulated distributions ofP i using the first-order Wasserstein distance (W 1 ). Lower values indicate closer agreement in behavioral heterogeneity, rather than merely similar mean cooperation levels (see Supplementary Information, Section S2.3). To probe context-dependent strategies, we analyze conditional cooperation as a func- tion of (i) an agentâs own previous action and (i) the fraction of cooperating neighbors in the previous round (for details, see Supplementary Information, Section S2.3.4). Acknowledgements. H.F.A. acknowledges ARAID for its financial support. A.A. acknowledges support from the grant RYC2021-033226-I funded by MICI- U/AEI/10.13039/501100011033 and the European Union NextGenerationEU/PRTR. A.A. and Y.M. were partially supported by the Government of Arag Ěon, Spain, and ERDF âA way of making Europeâ through grant E36-23R (FENOL), and by Grant No. PID2023-149409NB-I00 from Ministerio de Ciencia, Innovaci Ěon y Universidades, Agencia Espa Ěnola de Investigaci Ěon (MICIU/AEI/10.13039/501100011033) and ERDF âA way of making Europeâ. 17 Author contributions. H.F.A.: Conceptualization, Data curation, Software, For- mal analysis, Validation, Investigation, Visualization, Methodology, Writing-original draft, Writing-review and editing. C.G.L.: Conceptualization, Data curation, Inves- tigation, Methodology, Writing-review and editing. A.A.: Conceptualization, Investi- gation, Methodology, Writing-review and editing. Y.M.: Conceptualization, Formal analysis, Investigation, Methodology, Writing-original draft, Writing-review and edit- ing. Competing interests. The authors declare no competing interests. 18 Supplementary information S1 Technical details S1.1 Overview We simulate a repeated multi-agent Prisonerâs Dilemma on an input network. Each node is controlled either by an LLM-based agent or by a randomized policy. The CrewAI orchestration layer is used to orchestrate LLM agents, tasks, and responses, as well as an LLM served via an Ollama-compatible endpoint. In summary, in each round of the dynamics: 1. Each agent receives a brief description of the previous round, including neighborsâ choices and normalized payoffs. 2. Each agent replies with a JSON object containing a choice (âGREENâ or âBROWNâ) and a free-text reasoning. 3. The simulation computes payoffs for every node based on the payoff matrix, logs the answers, and proceeds to the next round. As for the graph input and experimental setup, we used the same graphs used in the original paper (i.e., [17]). This includes the network structures and the temporal changes of edges for the control cases, as well as numbers of iterations. S1.2 Agent creation and random agents For each node, we create either: ⢠an LLM-driven Agent (CrewAI), or ⢠a random agent (chooses uniformly between GREEN and BROWN). A fraction Ď of nodes are forced to be random. These nodes remain random throughout the dynamicsâ execution. S1.3 Payoff computation Payoffs are pairwise and defined by the following matrix (per neighbor): You / Neighbor GREEN BROWN GREEN70 BROWN100 The nodeâs payoff in a round equals the sum over its neighbors of the payoff according to each neighborâs choice. Implementation details: ⢠It iterates over neighbors of node i and sums the correct payoff depending on the nodeâs choice and the neighborâs choice. ⢠For comparison purposes, a normalized neighbor payoff is computed when exposing neighbor payoffs to agents. The normalization rescales a neighborâs payoff by the ratio of the focal nodeâs degree to the neighborâs degree. S1.4 LLM orchestration and prompts Important fields passed to each Agent instance include: ⢠role: short role description shown to the model. ⢠goal: short task: âmaximize your own ECUs.â ECU is a monetary unit. ⢠backstory: personalized per agent (degree and optional backstory fields). ⢠system template: a brief system message explaining the experiment. ⢠prompt template: in the implementation set to: Answer the task below in valid JSON format: input 19 The placeholder input is a dynamic variable that is automatically replaced at runtime with the task-specific description passed to the Task object. ⢠The CrewAI Task is created with expected output describing the required JSON format. The main LLM configuration values used by the script are: ⢠llm model = âollama/llama4:16x17bâ (default, but we also tested several LLMs). S1.4.1 Prompts and textual templates The textual templates used are reproduced below. These are sent to the agent together with the task description. Role and goal role = "You are participating in a social experiment as a volunteer. Your objective is to maximize your total ECUs, the monetary unit. An ECU is a monetary unit" goal = "Objective: **Your goal is to maximize your own ECUs (ECU is a monetary unit).** " System prompt system_template = "You are not expected to behave in any particular way. Whatever you do will determine the amount of money you can earn. You have a written version of this direction, which you can check at any stage of the experiment." Prompt template "Answer the task below in valid JSON format: input" Expected response format All agents must reply with a JSON object in the following format (this is enforced using the expected output string in the code): String in a Json format ONLY with choice (GREEN or BROWN) and reasoning (plain text): "choice": ..., "reasoning": "..." First-iteration description (displayed to the agent at round 0) A human-readable explanation of payoffs and examples is provided. Possible payoffs per neighbor. In the following table each row corresponds to the decision you can make and each column correspond to one of your neighborsâ decision. Your choice / Neighbor choice GREEN BROWN GREEN 7 0 BROWN 10 0 Consider that: you and each of your neighbors (which are humans) will globally earn more if you both choose GREEN (7 ECUs you/7 ECUs your neighbor); you will earn more if you choose BROWN and your neighbor chooses GREEN (10 ECUs you/0 ECUs your neighbor); but if both you and your neighbor choose BROWN you both will earn less (0 ECUs you/0 ECUs your neighbor) than if you both chose GREEN. This is the screen you will be seeing during the experiment (note that each participant actually sees the graph corresponding to his/her connectivity). Each round you must choose one of them clicking the corresponding button. 20 These are some examples of what you could earn in a round. Example 1: Imagine you choose GREEN, 3 of your neighbors choose GREEN and 1 chooses BROWN. In that round you will earn 3 x 7 + 1 x 0 = 21 ECUs. Example 2: In another round you choose BROWN, 2 of your neighbors choose GREEN and 2 choose BROWN. In that round you will earn 2 x 10 + 2 x 0 = 20 ECUs. Round iteration. Remember that each part will consist of an undetermined number of rounds. Each round you will have up to 20 s to choose a color. After these 20 s, if you didnât choose, the system will choose for you. Whatever happens it will not affect the behavior of the system in the next rounds: you will be able to make your subsequent choices normally. (Donât worry: 20 s are more than enough to make a choice). The round will not end until all participants have made their choice. At the end of each round you will see a screen like this one. Your choice (as given by the color) and your earning in this round. Also your [NUMBER_OF_NEIGHBORS] neighborsâ choices (represented by their colors) and their respective earnings in that round. Your neighborsâ earnings are given with respect to your number of neighbors. For example, you have 5 neighbors and one is Ferdinand (fictitious name). Ferdinand in turn has two neighbors: one is you and the other a stranger. If Ferdinand has won 10 ECUs in the last round, the gain of Ferdinand that you are shown is: (10 ECUs/2 neighbors of Ferdinand) x 5 neighbors of you = 25 ECUs. Note that what each of your neighbors has won depends on what you have chosen and also on what the neighbors of your neighbors have chosen. Immediately after finishing a round there will be a new one, and then another one, and so on until you see a screen warning you about the end of that part of the experiment. Note that the prompt placeholder [NUMBER OFNEIGHBORS] is replaced by the number of neighbors. In addition, the prompt explicitly stated that these neighbors corresponded to human participants (âwhich are humansâ), thereby indicating to the LLM that it was interacting with people in a real experiment. This framing substantially affected the resulting behavior. Backstory variants The code supports three backstory type values: ⢠STUDENT: adds âYou are a high school student.â ⢠GENDER: adds a âGender: man/woman.â line (based on vertex property is man from the data of [17]). ⢠STUDENTANDGENDER: combines the above. Round-to-round description For rounds r ⼠1, the agent receives a per-neighbor summary. The description contains: ⢠The round index. ⢠Each neighborâs previous choice and the neighborâs normalized payoff. ⢠The focal nodeâs previous choice and payoff. ⢠The instruction: âBased on this information, make your choice for this round.â Now, you have to answer again: This is round [ROUND_NUMBER]. In the previous round, your neighbors made the following choices: Neighbor choices in the previous round: Neighbor [NEIGHBOR_ID_1]: choice = [NEIGHBOR_CHOICE_1], payoff = [NORMALIZED_NEIGHBOR_PAYOFF_1]. Neighbor [NEIGHBOR_ID_2]: choice = [NEIGHBOR_CHOICE_2], payoff = [NORMALIZED_NEIGHBOR_PAYOFF_2]. Neighbor [NEIGHBOR_ID_3]: choice = [NEIGHBOR_CHOICE_3], payoff = [NORMALIZED_NEIGHBOR_PAYOFF_3]. Your choice and payoff in the previous round was: choice = [YOUR_PREVIOUS_CHOICE], payoff = [YOUR_PREVIOUS_PAYOFF]. Based on this information, make your choice for this round. Notice that in the simulations, all prompt placeholders were replaced by their corresponding values. The prompt shown here illustrates an agent interacting with three neighbors. 21 S1.5 Parsing and error handling ⢠The script expects the agent response to be valid JSON. Because LLMs occasionally output non-strict JSON, the code contains helpers to sanitize and robustly parse responses: â Normalize json string(s): trims text outside the outermost..., removes control characters, attempts to fix common issues (missing braces, extra whitespace between braces, etc.), and finally ensures the string looks like a JSON object. â Convert string to json: wraps json loads and, when necessary, raises a descriptive error on failure. â String to choice output (json string): uses the normalizer, parses the JSON, and returns a Python dict with keys choice and reasoning. Note that this is not the standard implementation. However, we implemented it this way to avoid minor issues common to less sophisticated LLMs. ⢠If parsing or the Crew kickoff call fails, the code retries up to max tries times (default maxtries = 10). On repeated failure, the implementation falls back to: â adefaultresponse "choice": "GREEN", "reasoning": "Error: ... defaulting to GREEN."; or This error message was implemented as a precaution to ensure robustness across different LLMs. In practice, however, such parsing failures were not observed for llama4 16x17b, the model used in the main experiments. S1.6 Simulation: main loop and execution logic The simulation proceeds for roundnum = 0, ..., maxrounds-1. Each round consists of the following steps: 1. Description preparation. ⢠For round 1, agents receive a full explanation of the payoff matrix, examples, and the experimental structure. ⢠For subsequent rounds (r > 1), agents receive: â each neighborâs previous choice, â each neighborâs normalized payoff, â their own previous choice and payoff. 2. Task construction. A Task object is created for each non-random agent. The task description is injected into the prompt template through the input placeholder, ensuring that the round-specific information is embedded into a fixed JSON-format instruction. 3. Decision generation. For each node i: ⢠If the node is random, a uniform random choice between GREEN and BROWN is generated. ⢠Otherwise, the CrewAI pipeline is executed: response = Crew(...).kickoff() parsed = string_to_choice_output(response.raw) 4. Retry mechanism and robustness. LLM calls may fail for two reasons: (a) Server or connection failures (e.g., temporary unavailability of the Ollama endpoint). (b) Formatting errors (the LLM does not produce valid JSON). To ensure robustness, each agent call is retried up to max tries = 10 times. This retry mechanism is important because: ⢠Server-side failures occasionally occur in long simulations. ⢠We empirically observed that when formatting errors occur, they are significantly more frequent in the first round. In later rounds, once the model has produced a correctly formatted JSON response, formatting mistakes become much less frequent. If all retries fail, the system assigns a safe default response: "choice":"GREEN", "reasoning":"Error... defaulting to GREEN." 22 This ensures that the simulation continues without interruption. 5. Payoff computation. After all nodes have chosen, payoffs are computed using the pairwise payoff matrix. High-level structured pseudocode initialize graph(s) g assign exactly floor(rho * N) random agents create LLM agents for non-random nodes for round_num in 0 .. max_rounds-1: prepare round-specific description for each node i: if random: result = random_choice() else: retry up to 10 times: call LLM attempt JSON parsing if all retries fail: result = default GREEN log result compute payoffs store payoffs in the graph update cooperation statistics end S2 Quantitative analysis S2.1 Definitions ⢠Sources of data: â Empirical: E. â Simulations: S (LLM players and LLM + random players). ⢠Notation: â D âE,S: dataset kind â c: Experimental condition (student, gender, and student + gender). â g: Network structure. â N : Number of network nodes. â r: Index for number of trials or replications in simulation. â i: Agent index. â t: Time step. â a D i (t): action C,D at time t for dataset D, where C and D are cooperate and defect, respectively. Letâs define the indicator variable: y D i (t) = 1 a D i (t) = C â0, 1,(4) thus, y = 1 means cooperation. S2.2 Macro analysis Cooperation rate (macro): For each condition c, topology g, and trial r, compute: X D c,g,r (t) = 1 N N X i=1 y D i (t),(5) 23 which is the fraction of cooperators at time t. Cooperation propensity per agent (micro): P D i = 1 T T X t=1 y D i (t),(6) i.e., how often an agent i cooperates over the full temporal horizon T . S2.2.1 Time evolution analysis Here, we test whether the model can reproduce the evolution of cooperation over time. To accomplish this, we compare simulated and empirical series. First, we compute the trajectory errors to determine whether they are equal at the same iteration. Then, we examine the temporal dependence structure of the empirical data using an autoregressive model. S2.2.2 Comparison of trajectory errors Before using an AR model to examine the temporal dynamics, we assess how closely the simulated series replicates the empirical series in terms of point-by-point accuracy. Two common metrics are used: ⢠Root Mean Squared Error (RMSE): RMSE = v u u t 1 T T X t=1 X E c,g,r (t)â X S c,g,r (t) 2 ,(7) which penalizes larger deviations more heavily. ⢠Mean Absolute Error (MAE): MAE = 1 T T X t=1 X E c,g,r (t)â X S c,g,r (t) ,(8) which measures the average absolute deviation and is less sensitive to outliers. While RMSE and MAE provide a measure of trajectory similarity between the two series, they do not capture the underlying temporal dependence structure. Therefore, we complement this analysis with AR model estimation to assess the persistence properties of the series. S2.2.3 Correlation analysis In addition to trajectory errors, we also compute the Pearson correlation coefficient between the empirical and simulated time series. While RMSE and MAE measure the magnitude of deviations at each time step, they do not capture whether the simulated series follows the same temporal pattern as the empirical data. Pearson correlation quantifies the degree of linear association between the two series, indicating whether increases and decreases in cooperation occur synchronously over time. We calculate Pearson both for the raw correlation series (r) and for the first differences of the series (r â ) to capture correlations in the changes between consecutive time steps. For each correlation, we report the corresponding p-value to evaluate the statistical significance of the observed association, providing a measure of confidence that the observed correlation is not due to random chance. S2.2.4 Autoregressive comparison between empirical and simulated series Before estimating the autoregressive model, we verify whether the empirical and simulated time series are stationary. Stationarity ensures that the mean, variance, and autocovariance structure of the series do not change over time, which is a fundamental assumption for AR model estimation. To do this, we apply the Augmented Dickey-Fuller (ADF) test to both series [24]. In our analysis, we estimate the model using the first-differenced series, âX D c,g,r (t), rather than the level series X D c,g,r (t). The differenced specification captures changes in cooperation over time and satisfies the stationarity requirement, whereas the level series exhibits non-stationary behavior according to the ADF test. 24 To assess whether the simulated time series reproduces the temporal dependence structure of the empirical data, we estimate an autoregressive model of order one (AR(1)) for the differenced series [23]. For a series âX D c,g,r (t), the AR(1) model is given by âX D c,g,r (t) = Îş + Ď âX D c,g,r (tâ 1) + Îľ(t),(9) where Îş is a constant, Ď measures persistence, and Îľ(t) is a white noise error term. We estimate the model separately for the empirical differenced series âX E c,g,r (t) and the simulated series âX S c,g,r (t). The key parameter of interest is the persistence coefficient Ď. To evaluate whether the simulation reproduces the dynamic structure, we use the following analysis metrics: ⢠Difference in persistence coefficients: To assess whether the simulated series reproduces the temporal dependence structure of the empirical data, we compare the estimated AR(1) persistence parameters through the absolute difference Diff Ď =|Ď S â Ď E |.(10) The coefficient Ď measures the strength and direction of temporal dependence. If|Ď S âĎ E | is small, the simulation reproduces the degree of persistence (or anti-persistence) observed in the empirical data. Large differences indicate that the simulated dynamics adjust either too quickly or too slowly relative to the empirical process, even if the overall trajectories appear similar. ⢠Implied mean increment: For a stationary AR(1) process in first differences (|Ď D | < 1), the unconditional mean of the differenced process is Îź D = Îş D 1â Ď D .(11) This quantity represents the average change around which the series increments fluctuate over time. Comparing Îź E and Îź S allows us to verify whether simulated and empirical processes exhibit similar average dynamics. Substantial differences indicate systematic biases in the evolution of cooperation, even if persistence properties are similar. S2.2.5 Null model Since the magnitude of the proposed metrics cannot be directly interpreted, we use a null model to evaluate their relevance. The null model generates a synthetic time series of the same length as the empirical and simulated series. Observations are independently drawn from the empirical marginal dis- tribution. This procedure preserves the unconditional distribution of cooperation levels while eliminating any temporal dependence structure. Let M obs denote the observed value of a metric. To construct a reference distribution under the null hypothesis of no temporal structure, we generate B independent null realizations and compute the corresponding metric valuesM (1) ,...,M (B) . Next, the statistical significance is assessed by comparing M obs to the empirical distribution ofM (b) B b=1 . Here, we use B = 10, 000. The null hypothesis is rejected when M obs lies in the tails of the null distribution. The Monte Carlo two-tailed p-value is defined as p = 1 B B X b=1 1 M (b) â Îź null ⼠M obs â Îź null ,(12) where Îź null = 1 B P B b=1 M (b) is the mean of the null distribution. Here, we adopt this approach because a small p two-tailed indicates that the observed metric is unlikely under the null hypothesis, given deviations in either direction. Note that for Pearson correlation, we adopted the null hypothesis of zero correlation (H 0 : r = 0). We adopt this two-tailed criterion because our goal is to determine whether the observed metric is unusually far from the null expectation, regardless of direction. Under this interpretation, both unusually large and unusually small values of M obs relative to the null distribution are treated as evidence against the null hypothesis. Thus, a small p-value indicates that the observed metric is unlikely under the null model, not only when it exceeds the null expectation, but more generally when it represents an extreme deviation from it. 25 S2.3 Micro analysis S2.3.1 Definitions Switching rate and persistence (micro): Letâs define switching rate: S D i = 1 T â 1 T X t=2 1 a D i (t)̸= a D i (tâ 1) .(13) Persistence is defined as P D i = 1â S D i . Neighborhood cooperation (for conditional cooperation): For each agent i and for t⼠2, we can calculate the previous round cooperation in the agentâs neighborhood N as m D i (tâ 1) = 1 N g tâ1 (i) X jâN g tâ1 (i) y D j (tâ 1),(14) where N g tâ1 (i) denotes the set of neighbors of agent i in the interaction graph g tâ1 . Also for agent i, we can use another variable z D i (tâ 1) = y D i (tâ 1) so that with z D i (tâ 1) and m D i (tâ 1) we can study conditional cooperation. For the conditional-rule analysis, the p-value is computed using the same Monte Carlo criterion defined in Eq. 12. However, each observation is represented by the triplet m D i (tâ 1),z D i (tâ 1),y D i (t) . Under the null hypothesis that y D i (t) is independent of the previous-round context, null realizations are generated by randomly permuting the values of y D i (t) across observations while keeping m D i (tâ 1) and z D i (tâ 1) fixed. This preserves the marginal distribution of current actions while removing any association between y D i (t) and the predictors m D i (tâ 1) and z D i (tâ 1) . S2.3.2 Confidence interval Due to computational constraints, for most tests, only one simulation per condition is performed. In this case, confidence intervals are estimated via bootstrapping across agents rather than simulation runs, capturing uncertainty due to agent heterogeneity. For the confidence interval (CI), we can use a bootstrap method to set pointwise uncertainty. In practice: ⢠Sample N agents/players with replacement (i.e., a player can appear multiple times). ⢠Recompute θ D (t) using the sampled agents. ⢠Repeat B times to obtain bootstrap replicates n θ D,(b) c,S (t) o B b=1 . ⢠Define CI 95% (t) = [ Q 0.025 ,Q 0.975 ] . To quantitatively measure the difference between these trajectories, we use RMSE and MAE to mea- sure the distance between the trajectories θ E c,S andθ S c,S . In addition, we calculate the Pearson correlation coefficient for both for the raw correlation series (r) and for the first differences of the series (r â ) between the trajectoriesθ E,(b) c,S andθ S,(b) c,S to capture correlations in the changes between consecutive time steps. For cases where we run the dynamics R times, we use the Cluster bootstrap to account for observations being grouped by simulation run. Specifically, for each bootstrap replicate, we first resample runs with replacement, then resample agents within each selected run. We then recompute the statistic of interest at each timestep. Repeating this many times yields a bootstrap distribution from which we take the percentiles to have the 95% confidence interval. This procedure preserves both within-run dependence and between-run variability, providing more realistic uncertainty than pooling all agents across runs as independent observations. S2.3.3 Heterogeneity across individuals To evaluate whether simulations reproduce the diversity of behavioral patterns observed in empirical data, we compare the distributions of individual cooperation propensities across agents. Let P E i N i=1 and P S i N i=1 26 denote the cooperation propensities computed for empirical and simulated datasets, respectively. These quantities capture the frequency with which each agent cooperates over the full temporal horizon and therefore encode behavioral heterogeneity at the individual level. We quantify the discrepancy between both distributions using the first-order Wasserstein distance. First, we sort both samples in non-decreasing order: P D (1) ⤠P D (2) â¤Âˇâ¤ P D (N ) , D âE,S. The empirical Wasserstein distance of order one between the two samples is then defined as W 1 (E,S) = 1 N N X k=1 P E (k) â P S (k) .(15) This metric measures the minimal average amount of probability mass that must be transported to transform the simulated distribution into the empirical one. A small value of W 1 (E,S) indicates that the simulation successfully reproduces the level of behavioral heterogeneity observed among real agents, whereas larger values reveal discrepancies such as overly homogeneous or excessively polarized simulated behaviors. S2.3.4 Conditional cooperation rule The conditional cooperation rule is defined as the probability that agent i cooperates at time t given its previous action and the cooperation level in its neighborhood: Pr y D i (t) = 1| m D i (tâ 1),z D i (tâ 1) ,(16) recall that z D i (tâ 1) = y D i (tâ 1). Let b k K k=0 denote the bin boundaries. Define I 1 = [b 0 ,b 1 ], I k = (b kâ1 ,b k ] for k = 2,...,K.(17) Then, for each previous action z â0, 1, define S k,z =(i,t) : m D i (tâ 1)â I k , z D i (tâ 1) = z.(18) The empirical conditional cooperation probability is estimated as b P D k,z = 1 |S k,z | X (i,t)âS k,z y D i (t).(19) Thus, for each bin k and previous action z â0, 1, we obtain four conditional cooperation probabilities: b P E k,0 , b P E k,1 : for experiments, b P S k,0 , b P S k,1 : for LLMs, after previous defection (z = 0) or cooperation (z = 1). Finally, we can do as before and compare these quantities using RMSE and MAE and Pearson correlation. 27 S3 Results S3.1 Visual analysis First, we show the results for Ď = 0. 1020304050 Iteration 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation Empirical: Lattice Empirical: Heterogeneous Simulation: Lattice Simulation: Heterogeneous Fig. S1: Evolution of cooperation over iterations for the STUDENT condition in the setting. Dashed lines show empirical results from the behavioral experiment, while solid lines show simulation results using llama416x17b. The solid lines represent the mean cooperation rate per round (empirical) or per iteration (simulation), and the error bars represent the standard error. 1020304050 Iteration 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation Empirical: Lattice Empirical: Heterogeneous Simulation: Lattice Simulation: Heterogeneous Fig. S2: Evolution of cooperation over iterations for the GENDER condition in the setting. Dashed lines show empirical results from the behavioral experiment, while solid lines show simulation results using llama416x17b. The solid lines represent the mean cooperation rate per round (empirical) or per iteration (simulation), and the error bars represent the standard error. 28 1020304050 Iteration 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation Empirical: Lattice Empirical: Heterogeneous Simulation: Lattice Simulation: Heterogeneous Fig. S3: Evolution of cooperation over iterations for the STUDENTANDGENDER condition in the setting. Dashed lines show empirical results from the behavioral experiment, while solid lines show simulation results using llama416x17b. The solid lines represent the mean cooperation rate per round (empirical) or per iteration (simulation), and the error bars represent the standard error. 1020304050 Iteration 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation Empirical: Lattice Empirical: Heterogeneous Simulation: Lattice Simulation: Heterogeneous Fig. S4: Control: Evolution of cooperation over iterations for the STUDENT condition in the setting. Dashed lines show empirical results from the behavioral experiment, while solid lines show simulation results using llama416x17b. The solid lines represent the mean cooperation rate per round (empirical) or per iteration (simulation), and the error bars represent the standard error. 29 1020304050 Iteration 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation Empirical: Lattice Empirical: Heterogeneous Simulation: Lattice Simulation: Heterogeneous Fig. S5: Control: Evolution of cooperation over iterations for the GENDER condition in the setting. Dashed lines show empirical results from the behavioral experiment, while solid lines show simulation results using llama416x17b. The solid lines represent the mean cooperation rate per round (empirical) or per iteration (simulation), and the error bars represent the standard error. 1020304050 Iteration 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation Empirical: Lattice Empirical: Heterogeneous Simulation: Lattice Simulation: Heterogeneous Fig. S6: Control: Evolution of cooperation over iterations for the STUDENTANDGENDER condition in the setting. Dashed lines show empirical results from the behavioral experiment, while solid lines show simulation results using llama416x17b. The solid lines represent the mean cooperation rate per round (empirical) or per iteration (simulation), and the error bars represent the standard error. 30 0.00.20.40.60.81.0 0 10 20 30 40 50 Number of Players a) Empirical: Lattice 0.00.20.40.60.81.0 0 10 20 30 40 50 60 b) Empirical: Heterogeneous 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 20 40 60 80 100 120 140 Number of Players c) Simulation: Lattice 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 10 20 30 40 50 60 d) Simulation: Heterogeneous Fig. S7: Distribution of per-agent cooperation rates for the STUDENT condition setting. Each histogram displays the fraction of cooperative rounds per player. Bins are uniformly spaced in the interval [0, 1]. Simulated results are generated using llama416x17b. 0.00.20.40.60.81.0 0 10 20 30 40 50 Number of Players a) Empirical: Lattice 0.00.20.40.60.81.0 0 10 20 30 40 50 60 b) Empirical: Heterogeneous 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 20 40 60 80 100 120 140 Number of Players c) Simulation: Lattice 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 10 20 30 40 50 60 d) Simulation: Heterogeneous Fig. S8: Distribution of per-agent cooperation rates for the GENDER condition setting. Each histogram displays the fraction of cooperative rounds per player. Bins are uniformly spaced in the interval [0, 1]. Simulated results are generated using llama416x17b. 31 0.00.20.40.60.81.0 0 10 20 30 40 50 Number of Players a) Empirical: Lattice 0.00.20.40.60.81.0 0 10 20 30 40 50 60 b) Empirical: Heterogeneous 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 20 40 60 80 100 120 140 Number of Players c) Simulation: Lattice 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 10 20 30 40 50 60 d) Simulation: Heterogeneous Fig. S9: Distribution of per-agent cooperation rates for the STUDENTANDGENDER condition setting. Each histogram displays the fraction of cooperative rounds per player. Bins are uniformly spaced in the interval [0, 1]. Simulated results are generated using llama416x17b. 0.00.20.40.60.81.0 0 10 20 30 40 50 60 Number of Players a) Empirical: Lattice 0.00.20.40.60.81.0 0 10 20 30 40 50 b) Empirical: Heterogeneous 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 50 100 150 200 250 Number of Players c) Simulation: Lattice 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 20 40 60 80 100 d) Simulation: Heterogeneous Fig. S10: Control: Distribution of per-agent cooperation rates for the STUDENT condition setting. Each histogram displays the fraction of cooperative rounds per player. Bins are uniformly spaced in the interval [0, 1]. Simulated results are generated using llama416x17b. 32 0.00.20.40.60.81.0 0 10 20 30 40 50 60 Number of Players a) Empirical: Lattice 0.00.20.40.60.81.0 0 10 20 30 40 50 b) Empirical: Heterogeneous 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 50 100 150 200 250 Number of Players c) Simulation: Lattice 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 20 40 60 80 100 d) Simulation: Heterogeneous Fig. S11: Control: Distribution of per-agent cooperation rates for the GENDER condition setting. Each histogram displays the fraction of cooperative rounds per player. Bins are uniformly spaced in the interval [0, 1]. Simulated results are generated using llama416x17b. 0.00.20.40.60.81.0 0 10 20 30 40 50 60 Number of Players a) Empirical: Lattice 0.00.20.40.60.81.0 0 10 20 30 40 50 b) Empirical: Heterogeneous 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 50 100 150 200 250 Number of Players c) Simulation: Lattice 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 20 40 60 80 100 d) Simulation: Heterogeneous Fig. S12: Control: Distribution of per-agent cooperation rates for the STUDENTANDGENDER condition setting. Each histogram displays the fraction of cooperative rounds per player. Bins are uniformly spaced in the interval [0, 1]. Simulated results are generated using llama416x17b. 33 In addition to testing different scenarios, we evaluated our simulation setup with llama416x17b across a range of temperatures. Specifically, temperature values were randomly sampled from the interval [0, 2]. These experiments were conducted exclusively on the lattice network topology (see Figure S13). 1020304050 Iteration 0.00 0.25 0.50 0.75 1.00 Cooperation a) Emp.: Lattice Sim.: Lattice 0.00.51.0 0 20 40 b) Emp.: Lattice 0.00.51.0 c) Sim.: Lattice Cooperative Rounds / Total Rounds Number of Players Fig. S13: Simulation results obtained with llama416x17b on the Lattice network under different tem- peratures sampled uniformly from the interval [0, 2]. To further assess whether the reduced dispersion observed in the simulations is model-dependent, we repeated the analysis using qwen332b. The results shown in Fig. S14 indicate that the tendency of LLM agents to generate less variable cooperation distributions than human participants is not specific to llama4 16x17b. 1020304050 Iteration 0.00 0.25 0.50 0.75 1.00 Cooperation a) Emp.: Lattice Sim.: Lattice Emp.: Heterogeneous Sim.: Heterogeneous 0 50 b) Emp.: Lattice c) Emp.: Heterogeneous 0.00.51.0 0 100 d) Sim.: Lattice 0.00.51.0 e) Sim.: Heterogeneous Cooperative Rounds / Total Rounds Number of Players Fig. S14: Simulation results obtained with qwen332b on the Lattice and Heterogeneous networks. Next, we present the results for other values of Ď. 34 1020304050 Iteration 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation 0 10 50 90 Rand. (%) Fig. S15: Cooperation over iterations for the lattice network, varying Ď, under the STUDENT condition. The simulation results were generated by llama416x17b. The black dashed line shows the empirical data, and the error bars show the standard error. 0.00.20.40.60.81.0 0 10 20 30 40 50 Number of Players a) Empirical 0.00.20.40.60.81.0 0 20 40 60 80 100 b) Hybrid Simulation (80% LLM + 20% Random) 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 20 40 60 80 100 120 140 Number of Players c) Full LLM Simulation 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 20 40 60 80 100 d) LLM Subset (within Hybrid) Fig. S16: Distribution of per-agent cooperation rates for the lattice network with Ď = 0.2 (80% LLM + 20% Random), under the STUDENT condition. Bins are uniformly spaced in the interval [0, 1]. Simulation results were generated using llama416x17b. 35 0.00.20.40.60.81.0 0 10 20 30 40 50 60 Number of Players a) Empirical 0.00.20.40.60.81.0 0 50 100 150 200 b) Hybrid Simulation (80% LLM + 20% Random) 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 50 100 150 200 250 Number of Players c) Full LLM Simulation 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 50 100 150 200 d) LLM Subset (within Hybrid) Fig. S17: Control: Distribution of per-agent cooperation rates for the lattice control network with Ď = 0.2 (80% LLM + 20% Random), under the STUDENT condition. Bins are uniformly spaced in the interval [0, 1]. Simulation results were generated using llama416x17b. 0.00.20.40.60.81.0 0 10 20 30 40 50 60 Number of Players a) Empirical 0.00.20.40.60.81.0 0 10 20 30 40 50 b) Hybrid Simulation (80% LLM + 20% Random) 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 10 20 30 40 50 60 Number of Players c) Full LLM Simulation 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 10 20 30 40 50 d) LLM Subset (within Hybrid) Fig. S18: Distribution of per-agent cooperation rates for the heterogeneous network with Ď = 0.2 (80% LLM + 20% Random), under the STUDENT condition. Bins are uniformly spaced in the interval [0, 1]. Simulation results were generated using llama416x17b. 36 0.00.20.40.60.81.0 0 10 20 30 40 50 Number of Players a) Empirical 0.00.20.40.60.81.0 0 20 40 60 80 b) Hybrid Simulation (80% LLM + 20% Random) 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 20 40 60 80 100 Number of Players c) Full LLM Simulation 0.00.20.40.60.81.0 Cooperative Rounds / Total Rounds 0 20 40 60 80 d) LLM Subset (within Hybrid) Fig. S19: Control: Distribution of per-agent cooperation rates for the heterogeneous control network with Ď = 0.2 (80% LLM + 20% Random), under the STUDENT condition. Bins are uniformly spaced in the interval [0, 1]. Simulation results were generated using llama416x17b. S3.2 Quantitative analysis S3.2.1 Macro analysis We begin our analysis by comparing the three backstory types, as shown in Figure S20 and Table S1. To reduce the number of simulations, we test two scenarios: (i) Hybrid Simulation (80% LLM + 20% Random) and (i) Full LLM Simulation. 37 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation a) Empirical STUDENT GENDER STUDENT AND GENDER 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation b) 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation c) 1020304050 Iteration 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation d) Fig. S20: Evolution of cooperation over iterations for the conditions in the setting. Dashed lines show empirical results from the behavioral experiment, while solid lines show simulation results using llama416x17b. The solid lines represent the mean cooperation rate per round (empirical) or per iter- ation (simulation), and the error bars represent the standard error. The panels depict the following scenarios: (a) lattice, (b) lattice control, (c) heterogeneous, and (d) heterogeneous control. Since the different backstory types do not significantly impact the results, we only consider STUDENT for the following results. First, we compare the different executions for both networks. See Figure S21 and Table S2. 38 ExperimentRMSEMAEDiff Ď Diff Îź r â lattice - GENDER0.104**0.0900.155***0.0010.577***0.367*** lattice - STUDENT0.098**0.083**0.217***0.0020.612***0.276* lattice - STUDENT AND GENDER0.110*0.0950.113***0.0020.495***0.384*** lattice control - GENDER0.089***0.074*0.383***0.003***0.472***0.151 lattice control - STUDENT0.089***0.074*0.464***0.002***0.583***0.214 lattice control - STUDENT AND GENDER0.089***0.075*0.327***0.002***0.506***0.031 heterogeneous - GENDER0.110**0.0980.0170.001*0.606***0.117 heterogeneous - STUDENT0.107**0.0960.0710.0020.622***0.151 heterogeneous - STUDENT AND GENDER0.108**0.0970.038*0.001*0.642***0.224 heterogeneous control - GENDER0.087***0.051*0.2760.005*0.467***0.002 heterogeneous control - STUDENT0.079***0.0510.3590.004*0.440***0.362*** heterogeneous control - STUDENT AND GENDER0.078***0.0500.4200.002*0.367***0.157 Table S1: Trajectory errors, autoregressive comparison, and correlation across different backstory types (for the scenario with no random agents). RMSE and MAE measure point-by-point deviations between empirical and simulated cooperation rates. Diff Ď denotes the absolute difference between the AR(1) persistence coefficients, and Diff Îź denotes the absolute difference between the implied long-run means of the simulated and empirical series. r denotes the Pearson correlation coefficient between empirical and simulated trajectories. For RMSE, MAE, Diff Ď , and Diff Îź , statistical significance is assessed using the Monte Carlo null model described in the text. For r, signif- icance refers to the null hypothesis of zero correlation (H 0 : r = 0). â denotes p < 0.01, â denotes p < 0.05, â denotes p < 0.1, and no symbol indicates that the result is not statistically significant. 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation a) 1020304050 Iteration 0.0 0.2 0.4 0.6 0.8 1.0 Cooperation b) Fig. S21: The average evolution of cooperation over r = 20 runs for both lattice and heterogeneous networks (panels a and b, respectively). The dashed lines show the results of the behavioral experiment, the error bars represent the standard error, and the solid lines show the simulation results. The shaded regions represent the respective standard deviation. 39 NetworkRMSEMAEDiff Ď Diff Îź r â Lattice0.099**0.0880.461***0.004**0.666***0.407*** Heterogeneous0.105**0.096*0.1900.001**0.708***0.179 Table S2: Trajectory errors, autoregressive comparison, and correla- tion across different executions for the scenario with no random agents, on both lattice and heterogeneous networks. RMSE and MAE mea- sure point-by-point deviations between empirical and average simulated cooperation rates. Diff Ď denotes the absolute difference between the AR(1) persistence coefficients, and Diff Îź denotes the absolute differ- ence between the implied long-run means of the simulated and empirical series. r denotes the Pearson correlation coefficient between empirical and simulated trajectories. For RMSE, MAE, Diff Ď , and Diff Îź , statisti- cal significance is assessed using the Monte Carlo null model described in the text. For r, significance refers to the null hypothesis of zero corre- lation (H 0 : r = 0). â denotes p < 0.01, â denotes p < 0.05, â denotes p < 0.1, and no symbol indicates that the result is not statistically sig- nificant. 40 Table S3 reports trajectory errors, AR(1) differences, and correlation metrics for simulations on lattice network, with varying percentages of random agents in the population. Table S4 presents the same metrics across different experimental scenarios (i.e., lattice, heterogeneous, and control setups). Random (%)RMSEMAEDiff Ď Diff Îź r â 00.098**0.083**0.217***0.0020.612***0.276* 100.116***0.099**0.537***0.000**0.665***0.284** 200.109***0.093**0.405***0.001***0.690***0.469*** 300.104***0.092**0.443***0.000***0.691***0.384*** 400.086*0.0740.098***0.0010.675***0.282** 500.057***0.044**0.376***0.0030.638***0.200 600.049**0.037**0.020***0.0030.569***0.236* 700.060***0.048***0.058*0.0030.645***0.332** 800.0910.0790.048*0.005***0.304**-0.097 900.1210.1100.215*0.004**0.457***0.297** Table S3: Trajectory errors, autoregressive comparison, and correlation across different percentages of random agents for the same backstory type on lattice network. RMSE and MAE measure point-by-point deviations between empirical and simulated cooperation rates. Diff Ď denotes the abso- lute difference between the AR(1) persistence coefficients, and Diff Îź denotes the absolute difference between the implied long-run means of the sim- ulated and empirical series. r denotes the Pearson correlation coefficient between empirical and simulated trajectories. For RMSE, MAE, Diff Ď , and Diff Îź , statistical significance is assessed using the Monte Carlo null model described in the text. For r, significance refers to the null hypothesis of zero correlation (H 0 : r = 0). â denotes p < 0.01, â denotes p < 0.05, â denotes p < 0.1, and no symbol indicates that the result is not statistically significant. ExperimentRMSEMAEDiff Ď Diff Îź r â lattice (100% LLM + 0% Rand.)0.098**0.083**0.217***0.0020.612***0.276* lattice (80% LLM + 20% Rand.)0.109***0.093**0.405***0.001***0.690***0.469*** lattice control (100% LLM + 0% Rand.)0.089***0.074*0.464**0.002***0.583***0.214 lattice control (80% LLM + 20% Rand.)0.087***0.046**0.519***0.005***0.509***0.098 heterogeneous (100% LLM + 0% Rand.)0.107**0.0960.0710.0020.622***0.151 heterogeneous (80% LLM + 20% Rand.)0.109***0.101**0.2490.000***0.743***0.343** heterogeneous control (100% LLM + 0% Rand.)0.079***0.0510.3590.004*0.440***0.362*** heterogeneous control (80% LLM + 20% Rand.)0.070**0.0460.3840.003*0.325**0.135 Table S4: Trajectory errors, autoregressive comparison, and correlation across experimental scenarios. RMSE and MAE measure point-by-point deviations between empirical and simulated cooperation rates. Diff Ď denotes the absolute difference between the AR(1) persistence coefficients, and Diff Îź denotes the absolute difference between the implied long-run means of the simulated and empirical series. r denotes the Pearson correlation coefficient between empirical and simulated trajectories. For RMSE, MAE, Diff Ď , and Diff Îź , statistical signifi- cance is assessed using the Monte Carlo null model described in the text. For r, significance refers to the null hypothesis of zero correlation (H 0 : r = 0). â denotes p < 0.01, â denotes p < 0.05, â denotes p < 0.1, and no symbol indicates that the result is not statistically significant. S3.2.2 Micro analysis First, we show the persistence across different percentages of random agents (in the lattice scenario), as shown in Table S5. Next, we present the comparison between the different experiments (see Table S6). 41 ExperimentP D 95% CI Empirical0.640(0.629, 0.652) Lattice (100% LLM + 0% Rand.)0.579(0.573, 0.586) Lattice (90% LLM + 10% Rand.)0.613(0.605, 0.621) Lattice (80% LLM + 20% Rand.)0.633(0.624, 0.642) Lattice (70% LLM + 30% Rand.)0.645(0.634, 0.656) Lattice (60% LLM + 40% Rand.)0.650(0.637, 0.662) Lattice (50% LLM + 50% Rand.)0.610(0.598, 0.621) Lattice (40% LLM + 60% Rand.)0.593(0.582, 0.605) Lattice (30% LLM + 70% Rand.)0.564(0.554, 0.574) Lattice (20% LLM + 80% Rand.)0.540(0.531, 0.549) Lattice (10% LLM + 90% Rand.)0.518(0.511, 0.525) Table S5: Mean cooperation persistence and 95% confidence interval for different percentages of random agents. Confidence intervals are estimated via boot- strapping over agents. ExperimentP D 95% CI lattice (100% LLM + 0% Rand.)0.579(0.573, 0.586) lattice (80% LLM + 20% Rand.)0.633(0.623, 0.642) lattice - Empirical0.640(0.629, 0.651) lattice control (100% LLM + 0% Rand.)0.610(0.604, 0.615) lattice control (80% LLM + 20% Rand.)0.589(0.582, 0.596) lattice control - Empirical0.724(0.711, 0.738) heterogeneous (100% LLM + 0% Rand.)0.604(0.592, 0.617) heterogeneous (80% LLM + 20% Rand.)0.629(0.615, 0.642) heterogeneous - Empirical0.620(0.610, 0.630) heterogeneous control (100% LLM + 0% Rand.)0.578(0.572, 0.584) heterogeneous control (80% LLM + 20% Rand.)0.577(0.570, 0.584) heterogeneous control - Empirical0.693(0.680, 0.707) Table S6: Mean cooperation persistence and 95% confidence interval. Confidence intervals are estimated via bootstrapping over agents. Next, we examine conditional cooperation as a function of neighborhood behavior. Figures S22 and S23 show the average probability of cooperation as a function of the fraction of cooperating neighbors in the previous round, separately after a previous defection and after a previous cooperation. Both tested scenarios (lattice and heterogeneous) comprise 20 dynamics executions. 42 0.0 0.2 0.4 0.6 0.8 1.0 Neighborhood Coop. (After Coop.) a) 1020304050 Timestep 0.0 0.2 0.4 0.6 0.8 1.0 Neighborhood Coop. (After Def.) b) Fig. S22: Average fraction of cooperating neighbors over time for the lattice network. Results are shown separately for agents that cooperated (panel a) or defected (panel b) in the previous round. The continuous blue line represents the simulation, and the black dashed line represents the empirical data. The simulated curve is the mean across all 20 runs, and the error bars represent 95% CIs estimated by Cluster bootstrap across runs. Empirical patterns (dashed) are shown with 95% CIs estimated by bootstrapping over players. 43 0.0 0.2 0.4 0.6 0.8 1.0 Neighborhood Coop. (After Coop.) a) 1020304050 Timestep 0.0 0.2 0.4 0.6 0.8 1.0 Neighborhood Coop. (After Def.) b) Fig. S23: Average fraction of cooperating neighbors over time for the heterogeneous network. Results are shown separately for agents that cooperated (panel a) or defected (panel b) in the previous round. The continuous green line represents the simulation, and the black dashed line represents the empirical data. The simulated curve is the mean across all 20 runs, and the error bars represent 95% CIs estimated by Cluster bootstrap across runs. Empirical patterns (dashed) are shown with 95% CIs estimated by bootstrapping over players. 44 S3.2.3 Heterogeneity across individuals To further analyze the role of stochastic behavior in shaping individual-level heterogeneity, we vary the fraction of random agents Ď in the population. Table S7 reports the resulting values of W 1 (E,S) across different percentages of random agents. The results show a non-monotonic relationship between Ď and heterogeneity mismatch, with intermediate levels of randomness generally yielding lower distances between empirical and simulated distributions. Random (%)W 1 (E, S) 00.101 100.089 200.081 300.084 400.073 500.037 600.039 700.051 800.080 900.119 Table S7: Sensitivity of individual coopera- tionheterogeneityto the fraction of random agents. Wasserstein dis- tance between empirical andsimulateddistri- butionsofindividual cooperationpropen- sities.Lowervalues indicate better repro- duction of behavioral heterogeneity. S3.2.4 Conditional cooperation rule To quantify how well the simulations reproduce the empirical conditional cooperation rule, we compute RMSE and MAE between the empirical and simulated conditional cooperation probabilities. Errors are calculated separately for cases in which the previous action was defection (z = 0) and cooperation (z = 1), as well as aggregated across both cases, see Table S8. 45 Random (%)RMSE z=0 MAE z=0 RMSE z=1 MAE z=1 00.237**0.217***0.397***0.388* 100.208**0.188***0.322***0.310* 200.1800.165***0.255***0.246* 300.1530.144***0.199***0.192* 400.145**0.136***0.117***0.116*** 500.140***0.121***0.096***0.092*** 600.137***0.112***0.060***0.055*** 700.135***0.101***0.049***0.041*** 800.150***0.116***0.042***0.038*** 900.181***0.163***0.035***0.027*** Table S8: Deviation between empirical and simulated con- ditional cooperation rules across different percentages of random agents. Errors are computed separately for cases in which the previous action was defection (z = 0) or cooperation (z = 1). The columns RMSE all and MAE all summarize the overall deviation across both cases. â denotes p < 0.01, â denotes p < 0.05, â denotes p < 0.1, and no symbol indicates that the result is not statistically significant. 46 References [1] Lu, Y., Aleta, A., Du, C., Shi, L., Moreno, Y.: Llms and generative agent-based models for complex systems research. Physics of Life Reviews 51, 283â293 (2024) [2] Chiang, C.-H., Lee, H.-y.: Can large language models be an alternative to human evaluations? In: Proceedings of the 61st Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), p. 15607â15631 (2023) [3] Papachristou, M., Yang, L., Hsu, C.-C.: Leveraging large language models for collective decision- making. Proceedings of the ACM on Human-Computer Interaction 9(7), 1â44 (2025) [4] Flamino, J., Modi, M.S., Szymanski, B.K., Cross, B., Mikolajczyk, C.: Testing the limits of large language models in debating humans. Scientific Reports 15(1), 13852 (2025) [5] Nudo, J., Pandolfo, M.E., Loru, E., Samory, M., Cinelli, M., Quattrociocchi, W.: Generative exag- geration in llm social agents: Consistency, bias, and toxicity. Online Social Networks and Media 51, 100344 (2026) [6] Loru, E., Nudo, J., Di Marco, N., Santirocchi, A., Atzeni, R., Cinelli, M., Cestari, V., Rossi-Arnaud, C., Quattrociocchi, W.: The simulation of judgment in llms. Proceedings of the National Academy of Sciences 122(42), 2518443122 (2025) [7] Hadad, O., Loru, E., Nudo, J., Di Marco, N., Cinelli, M., Quattrociocchi, W.: The statistical signature of llms. arXiv preprint arXiv:2602.18152 (2026) [8] Akata, E., Schulz, L., Coda-Forno, J., Oh, S.J., Bethge, M., Schulz, E.: Playing repeated games with large language models. Nature Human Behaviour 9(7), 1380â1390 (2025) [9] Ashery, A.F., Aiello, L.M., Baronchelli, A.: Emergent social conventions and collective bias in llm populations. Science Advances 11(20), 9368 (2025) [10] Kumar, A., Poungpeth, N., Yang, D., Farrell, E., Lambert, B.L., Groh, M.: When large language models are reliable for judging empathic communication. Nature Machine Intelligence 8(2), 173â185 (2026) [11] Serapio-Garc ĚÄąa, G., Safdari, M., Crepy, C., Sun, L., Fitz, S., Romero, P., Abdulhai, M., Faust, A., Matari Ěc, M.: A psychometric framework for evaluating and shaping personality traits in large language models. Nature Machine Intelligence, 1â15 (2025) [12] Wang, Y., Zhong, W., Li, L., Mi, F., Zeng, X., Huang, W., Shang, L., Jiang, X., Liu, Q.: Aligning large language models with human: A survey. arXiv preprint arXiv:2307.12966 (2023) [13] Slocum, S., Parker-Sartori, A., Hadfield-Menell, D.: Diverse preference learning for capabilities and alignment. In: The Thirteenth International Conference on Learning Representations (2025). https://openreview.net/forum?id=pOq9vDIYev [14] Hamilton, S.: Detecting mode collapse in language models via narration. In: Proceedings of the First Edition of the Workshop on the Scaling Behavior of Large Language Models (SCALE-LLM 2024), p. 65â72 (2024) [15] Reia, S.M., Pfoser, D., et al.: Opportunities and challenges of llms in urban science: Comment onâ llms and generative agent-based models for complex systems researchâ by yikang lu et al. Physics of Life Reviews 53, 305â306 (2025) [16] Ye, J., Cao, L., Chen, D., Ferrara, E.: Stop drawing scientific claims from llm social simulations without robustness audits. arXiv preprint arXiv:2605.18890 (2026) [17] Gracia-L Ěazaro, C., Ferrer, A., Ruiz, G., Taranc Ěon, A., Cuesta, J.A., S Ěanchez, A., Moreno, Y.: Hetero- geneous networks do not promote cooperation when humans play a prisonerâs dilemma. Proceedings of the National Academy of Sciences 109(32), 12922â12926 (2012) 47 [18] Gruji Ěc, J., Fosco, C., Araujo, L., Cuesta, J.A., S Ěanchez, A.: Social experiments in the mesoscale: Humans playing a spatial prisonerâs dilemma. PloS one 5(11), 13749 (2010) [19] Touvron, H., Lavril, T., Izacard, G., Martinet, X., Lachaux, M.-A., Lacroix, T., Rozi`ere, B., Goyal, N., Hambro, E., Azhar, F., et al.: Llama: Open and efficient foundation language models. arXiv preprint arXiv:2302.13971 (2023) [20] Bi, X., Chen, D., Chen, G., Chen, S., Dai, D., Deng, C., Ding, H., Dong, K., Du, Q., Fu, Z., et al.: Deepseek llm: Scaling open-source language models with longtermism. arXiv preprint arXiv:2401.02954 (2024) [21] Yang, A., Li, A., Yang, B., Zhang, B., Hui, B., Zheng, B., Yu, B., Gao, C., Huang, C., Lv, C., et al.: Qwen3 technical report. arXiv preprint arXiv:2505.09388 (2025) [22] Allal, L.B., Lozhkov, A., Bakouch, E., Bl Ěazquez, G.M., Penedo, G., Tunstall, L., Marafioti, A., Kydl ĚÄąËcek, H., Lajar ĚÄąn, A.P., Srivastav, V., et al.: Smollm2: When smol goes bigâdata-centric training of a small language model. arXiv preprint arXiv:2502.02737 (2025) [23] Hamilton, J.D.: Time Series Analysis. Princeton University Press, Princeton, NJ (1994) [24] Dickey, D.A., Fuller, W.A.: Distribution of the estimators for autoregressive time series with a unit root. Journal of the American statistical association 74(366a), 427â431 (1979) 48