DOI: 10.5061/dryad.c866t1gkr
title: Multifractal analysis of microscopic fluctuation and quasi-equilibrium landscape of evolutionary games with environmental feedback on complex networks

description: The evolution of societies within civilizations is characterized by the
 emergence of cooperative traits through continuous adaptation within their
 unique complex networks. Cooperative behaviors can sustain local
 persistence for extended durations before succumbing to total defection,
 motivating this study to investigate evolutionary games with environmental
 feedback across diverse network topologies, integrating replicator
 dynamics with resource-dependent payoffs and agent-based simulations. By
 exploring and quantifying macroscopic (potential landscape) and
 microscopic (multifractal fluctuation) perspectives, the research
 elucidates the influence of network structure and population size on the
 transition from oscillatory to quasi-equilibrium cooperation regimes. This
 study introduces the innovative application of multifractal detrended
 fluctuation analysis (MF-DFA) to evolutionary game dynamics, revealing
 that intermittent, scale-invariant bursts in cooperation are not random
 noise but carry predictive information about metastable states. Employing
 MF-DFA quantifies these bursts of cooperative behavior, demonstrating how
 network size affects cooperation’s stability and potential landscape.
 Findings reveal that larger networks exhibit deeper and narrower potential
 wells, enhancing metastability and prolonging the endurance of cooperative
 interactions. These results advance understanding of cooperative evolution
 within finite, structured populations, providing valuable insights into
 mechanisms sustaining cooperation in real-world adaptive systems and
 suggesting a novel method to detect early-warning signals of cooperation
 collapse or stability.
description: # Multifractal analysis of microscopic fluctuation and quasi-equilibrium
 landscape of evolutionary games with environmental feedback on complex
 networks Dataset DOI: [10.5061/dryad.c866t1gkr](10.5061/dryad.c866t1gkr)
 ## Description of the data and file structure The code implements an
 agent-based simulation framework grounded in evolutionary game theory to
 study the coevolution of cooperative behavior and environmental feedback
 on complex networks. The simulations model populations of agents
 distributed across various network topologies, including Toroidal lattice,
 Erdős-Rényi, Watts-Strogatz, and Barabási-Albert networks. Agents engage
 in pairwise evolutionary games with dynamic, resource-dependent payoff
 matrices that evolve based on the state of a shared environmental
 resource. Strategy updates follow the Fermi pairwise comparison rule,
 capturing bounded rationality and local interactions. The code generates
 time series data of cooperative fractions, environmental resource levels,
 spatial correlations, and strategic diversity across extended simulation
 runs. These datasets underpin the multifractal detrended fluctuation
 analysis (MF-DFA) and potential landscape calculations, providing
 microscopic and macroscopic insights into the emergence and stability of
 quasi-equilibrium cooperation regimes under varying population sizes and
 network structures. This computational framework supports reproducible
 exploration of eco-evolutionary feedback mechanisms in structured
 populations, facilitating analysis of intermittent cooperative bursts,
 metastability, and early-warning signals for cooperation collapse or
 persistence. ### Files and variables #### File:
 For_figure_1-2_metrices_for_different_network.py **Description:** Main
 script implementing agent-based evolutionary game simulations on various
 network topologies with environmental feedback. Captures cooperator
 fraction time series and related metrics. **The output files that will be
 generated upon running this code:** **`Lattice_metrices.csv`**\ Contains
 simulation results for the feedback-evolving game on a Toroidal lattice
 network topology. Generated by the lattice simulation code section.
 **`ER_metrices.csv`**\ Contains simulation results for the
 feedback-evolving game on an Erdős-Rényi random network topology.
 Generated by the ER simulation code section. **`BA_metrics.csv`**\
 Contains simulation results for the feedback-evolving game on a
 Barabási-Albert scale-free network topology. Generated by the BA
 simulation code section. **`WS_metrices.csv`**\ Contains simulation
 results for the feedback-evolving game on a Watts-Strogatz small-world
 network topology. Generated by the WS simulation code section. Each CSV
 file contains rows of time series data collected over multiple simulation
 episodes. **Variables and Abbreviations** Each simulation saves the
 following columns into the CSV files (all variables are dimensionless
 unless otherwise noted): | **Variable Name** | **Description** | **Units**
 | **Comments** | | :----------------------- |
 :--------------------------------------------------------------------------------------- | :--------------------------------- | :-------------------------------------------------------------------------------------------------------------------- | | time | Simulation time at each recorded time step | Arbitrary units (simulation steps) | Evenly spaced over the simulation duration (e.g., from 0 to 5000) | | real\_eig1 | Real part of the first eigenvalue of the Jacobian matrix at the current state | Dimensionless | Stability analysis of game dynamics | | real\_eig2 | Real part of the second eigenvalue of the Jacobian matrix | Dimensionless | 
 
 | | imag\_eig1 | Imaginary part of the first eigenvalue | Dimensionless |
 
 
 | | imag\_eig2 | Imaginary part of the second eigenvalue | Dimensionless
 | 
 
 | | cooperator\_fraction (x) | Fraction of agents currently adopting the
 cooperative strategy ("C") | Between 0 and 1 | Calculated as the
 number of cooperators divided by the total population size | |
 random\_expectation | Expected fraction of edges connecting agents of the
 same strategy assuming random mixing | Between 0 and 1 | Calculated
 as x2+(1−x)2, where x is the global cooperator fraction | |
 resource\_level (m) | The current environmental resource abundance level
 is influenced by agents’ behavior | Between 0 and 1 | Models replete (m=1)
 to depleted (m=0) resource states | | spatial\_correlation | Fraction of
 edges connecting agents with identical strategies | Between 0 and 1 |
 Indicates clustering of cooperators or defectors on the network | |
 episode | Identifier for the simulation run (episode number) | Integer |
 Multiple episodes are simulated for averaging. Episodes are numbered
 starting from 1; averaged data is marked with -1 | **Simulation Parameters
 (Constants)** R0, S0, T0, P0: Payoff matrix parameters in depleted
 environment R0=5.0, S0=1.0, T0=3.0, P0=0.0 R1, S1, T1, P1: Payoff matrix
 parameters in repleted environment R1=3.0, S1=0.0, T1=5.0, P1=1.0 theta:
 Cooperation enhancement parameter for environmental feedback (value: 2.0)
 epsilon: Speed of environmental change (value: 0.1) beta: Selection
 strength parameter related to Fermi update (value used: 10) ***Initial
 values for cooperator_fraction (x) =0.5 and resource_level (m) =0.4 are
 taken for all the scripts and simulations.** **Missing Values** The CSV
 files do not contain missing values; all variables are recorded at each
 timestep for every episode. Averaged rows across episodes are appended at
 the end of each file, with the episode field set to -1 to distinguish
 averaged values from individual runs. #### File: For_figure_3_snapshot.py
 **Description:**  **The output files that will be generated upon running
 this code:** **snapshot_ep{episode}_t{time_step}.png**\ A series of
 snapshot image files depicting the spatial distribution of agents'
 strategies (Cooperator or Defector) on a Toroidal lattice network at
 discrete time points during simulation episodes. {episode}: Simulation
 episode number (e.g., 0 for the first run). {time_step}: Simulation time
 step (in seconds). Snapshots are saved every 200 time steps up to 10,000.
 These images visually represent evolving cooperative (blue nodes) and
 defective (dark red nodes) clusters in a population of 10,000 agents,
 illustrating the dynamics of strategy changes and environmental feedback
 over time. **Simulation Constants** | **Variable** | **Description** |
 **Units** | **Comments** | | :------------- |
 :----------------------------------------------------------- |
 :---------------------------- | :----------------------------------- | |
 R0, S0, T0, P0 | Payoff matrix parameters for the depleted environment |
 Dimensionless | R0=5.0, S0=1.0, T0=3.0, P0=0.0 | | R1, S1, T1, P1 | Payoff
 matrix parameters for the replete environment | Dimensionless |
 R1=3.0, S1=0.0, T1=5.0, P1=1.0 | | theta | Cooperation enhancement
 parameter for environmental feedback | Dimensionless | Value: 2.0 | |
 epsilon | Speed of environmental change | Dimensionless (per time step) |
 Value: 0.1 | | beta | Selection intensity parameter (inverse of noise) |
 Dimensionless | Used as 10 to define the noise level | **Key Variables
 During Simulation** | Variable | Data Type | Description | Units |
 Comments | | :------------------------- | :-------------- |
 :---------------------------------------------------------------------------------------- | :------------- | :-------------------------------------------- | | m | Float | Environmental resource level | Dimensionless  | Updated every step; influences payoffs | | A\_evo | 2×2 numpy array | Current payoff matrix constructed as a weighted sum of depleted and repleted environments | Dimensionless | A(m)=(1−m) A0+mA1 | | Dg\_prime, Dr\_prime | Float | Gamble-Intending Dilemma (GID) and Risk-Averting Dilemma (RAD) strengths | Dimensionless | Computed from the current payoff matrix | | Agents | List of objects | Collection of agent instances each containing: strategy ("C" or "D"), payoff, neighbors | N/A | Individual decision-making and updates | | strategy (per agent) | String | Current strategy: "C" for Cooperator, "D" for Defector | N/A | Used to color nodes in snapshots | | point (per agent) | Float | Accumulated payoff from interactions with neighbors | Dimensionless | Used in Fermi update rule | | next\_strategy (per agent) | String | Strategy to adopt in next time step after probabilistic update | N/A | Updated every step based on payoff comparison | | neighbors\_id (per agent) | List\[int] | IDs of neighboring agents on the Toroidal lattice | N/A | Defined by network topology | | t | numpy array | Array of simulation time values (time steps) | Seconds | Range: 0 to 10,000 steps, step size 1 | | episode | Integer | Simulation run identifier | N/A | Typically, 0 or multiple independent runs | **Functions with Important Roles** A(m): Returns the dynamic payoff matrix based on the current resource level m. calculate_dilemma_strengths(A_evo): Computes dilemma strengths (Dg_prime, Dr_prime) from payoff matrix. pw_fermi: Implements the pairwise Fermi update rule for strategy adoption reflecting bounded rationality. snapshot_colored: Saves spatial snapshot images showing agent strategies on the lattice. **Network and Simulation Setup** **Network type:** Toroidal lattice (2D grid with periodic boundary conditions). **Population size:** 10,000 or 625 agents placed on lattice nodes. **Neighborhood:** Each agent interacts with its 4 von Neumann neighbors (up, down, left, right). **Initial Condition:** Approximately 50% agents randomly initialized as cooperators ("C"). **Simulation Duration:** 10,000 time steps with snapshots recorded every 200 steps. **Episodes:** Typically 1 or more independent simulation runs for statistical robustness. Missing Values and Data Integrity No missing values are present in the simulation data. All variables are updated and recorded at every simulation step. Snapshot images always depict the full lattice state at the respective time points without missing nodes. 
 
 #### File: For_figure_4_**Panel A – Quasi-stable Phase Analysis.py**
 **Description:** **The output files that will be generated upon running
 this code: cooperation_fractions.csv**\ This CSV file contains processed
 data representing the evolution of cooperative fractions over simulation
 time across multiple population sizes. Each row corresponds to a recorded
 simulation time step. Columns include: * time (simulation time points,
 units are arbitrary simulation steps), *
 cooperation_fraction_ corresponding to the averaged
 cooperative fraction over 100 independent episodes for each population
 size. For example, cooperation_fraction_100 for a population size of 100
 individuals. **Variables and Abbreviations** | **Variable Name** |
 **Description** | **Units** | **Comments** | | :--------------------------
 | :----------------------------------------------------- |
 :--------------------------- |
 :----------------------------------------------------------------------- |
 | time | Discrete simulation time points | Arbitrary (simulation steps) |
 Linearly spaced from 0 to 1,000,000 over 100,000 intervals. | |
 cooperation\_fraction\_\ | Average fraction of agents cooperating
 at given time t | Dimensionless  | Computed as the mean over 100
 simulation episodes for population size N. | **Key Parameters and
 Constants (used internally in the code, not directly saved in a file):** |
 **Parameter** | **Description** | **Units** | **Typical Value / Comments**
 | | :----------------- |
 :------------------------------------------------------------- |
 :------------ | :---------------------------------- | | R0, S0, T0, P0 |
 Payoffs for depleted environment game matrix components | Dimensionless |
 5.0, 1.0, 3.0, 0.0 | | R1, S1, T1, P1 | Payoffs for repleted environment
 game matrix components | Dimensionless | 3.0, 0.0, 5.0, 1.0 | | theta |
 Cooperation enhancement factor for environmental feedback | Dimensionless
 | 2.0 | | epsilon | Environmental change rate | Dimensionless | 0.1 | |
 beta | Selection intensity parameter for Fermi update (inverse noise) |
 Dimensionless | 10 | | population\_sizes | Population sizes simulated |
 Integer | 100, 225, 400, 625, 900, 1225, 1600 | | episodes\_per\_pop |
 Number of independent simulation episodes per population | Integer | 100 |
 **Internal Simulation Variables (relevant for understanding the data
 generation)** | **Variable** | **Description** | **Units** | **Comments**
 | | :----------------------------- |
 :------------------------------------------------------------------------
 | :------------- |
 :--------------------------------------------------------------------------------------- | | m | Environmental resource level (dynamic during simulation) | Dimensionless  | Updated each time step based on cooperation level and feedback dynamics. | | A\_evo | Dynamic payoff matrix depending on m | Dimensionless | Weighted average of payoff matrices for depleted and repleted environments (2×2 matrix). | | Dg\_prime, Dr\_prime | Gamble-Intending Dilemma and Risk-Averting Dilemma strengths respectively | Dimensionless | Calculated from A\_evo, used to update agent payoffs during strategy evolution. | | Agent.strategy | Strategy adopted by each agent ("C" for Cooperation, "D" for Defection) | String | Updated via pairwise Fermi process at each time step. | | Agent.point | Accumulated payoff of an agent from interactions with neighbors | Dimensionless | Used to calculate adoption probabilities for strategies. | | episode\_cooperator\_fractions | Time series data of cooperator fractions per episode | Dimensionless  | Recorded at every simulation time step for averaging and further statistical analysis. | Data Format and Missing Values * The CSV file uses comma  as the delimiter with headers indicating variable names. * No missing values are expected or represented in the dataset; each column contains complete time series data for the respective population size. * All recorded values are numeric and represent model outputs with no placeholders such as NA or empty cells. Additional Notes * The simulation runs 100 independent episodes per population size to ensure statistical robustness. * The cooperation_fraction columns represent averaged values over these episodes, smoothing stochastic fluctuations inherent in agent-based simulations. * Time steps cover a long interval (0 to 1,000,000 simulation steps) with a fine resolution (100,000 steps), enabling detailed temporal analysis of quasi-stable cooperation phases. * The code uses multiprocessing (multiprocessing library) for parallel execution across episodes, improving computational efficiency.   #### File: For_figure_4_panel_B-Entropy.py **Description:**  **The output files that will be generated upon running this code**: **mean_entropy.csv** This CSV file contains the average Shannon entropy values calculated over time from multiple independent simulation episodes for various population sizes. Rows correspond to discrete time points in the simulation. Columns include: ·        Time: The simulation time steps at which entropy was recorded (units in simulation steps). ·        Entropy_: The mean entropy value (averaged over 100 episodes) for the respective population size. For example, Entropy_100 gives the average entropy for a population size 100 agents. **Variables and Abbreviations** | **Variable Name** | **Description** | **Units** | **Comments** | | :---------------- | :-------------------------------------------------------------------------------------------- | :--------------------------- | :--------------------------------------------------------------------------------------------------------------------------------------------------------- | | Time | Simulation time points at which entropy is measured | Arbitrary (simulation steps) | Recorded every save\_every timesteps (default every 10 steps) during the simulation | | Entropy\_\ | Average Shannon entropy of the population strategy distribution at time for population size N | Bits (Shannon entropy) | Calculated as H = − 𝑝𝑐𝑜𝑜𝑝×log2 (𝑝𝑐𝑜𝑜𝑝)− 𝑝𝑑𝑒𝑓 ×log2 (𝑝𝑑𝑒𝑓) where 𝑝𝑐𝑜𝑜𝑝 and 𝑝𝑑𝑒𝑓 are the proportions of cooperators and defectors | **Key Parameters and Constants (used internally in the code but not saved to file):** | **Parameter** | **Description** | **Units** | **Typical Value / Comments** | | :---------------- | :------------------------------------------------------------ | :------------ | :---------------------------------- | | R0, S0, T0, P0 | Payoff matrix entries for the depleted environment game | Dimensionless | 5.0, 1.0, 3.0, 0.0 | | R1, S1, T1, P1 | Payoff matrix entries for repeated environment game | Dimensionless | 3.0, 0.0, 5.0, 1.0 | | theta | Environmental cooperation feedback factor | Dimensionless | 2.0 | | epsilon | Speed of environmental change | Dimensionless | 0.1 | | beta | Inverse temperature parameter controlling selection intensity | Dimensionless | 10 | | population\_sizes | List of simulated population sizes | Integer | 100, 225, 400, 625, 900, 1225, 1600 | | num\_episodes | Number of independent simulation repeats per population size | Integer | 100 | | total\_timesteps | Total number of simulation steps per episode | Integer | 1,000,000 | | save\_every | Interval of simulation steps between entropy recordings | Integer | 10 | **Internal Simulation Variables (for understanding data generation):** | **Variable** | **Description** | **Units** | **Comments** | | :------------- | :----------------------------------------------------------------------- | :----------------- | :------------------------------------------------------------------------------- | | agents | List of Agent objects representing individuals | N/A | Each agent holds strategy, payoff, and neighbor information | | Agent.strategy | Current strategy of the agent ("C" for Cooperator, "D" for Defector) | Categorical string | Evolving via pairwise Fermi update | | m | Environmental resource level | Dimensionless | Dynamically updated during simulation based on cooperation level | | A\_evo | Environment-dependent payoff matrix | 2×2 NumPy array | Combined payoff matrix weighted by resource level | | Dg, Dr | Gamble-Intending and Risk-Averting Dilemma strengths derived from A\_evo | Dimensionless | Used to compute payoffs in agent interactions | | entropy\_list | List of Shannon entropy values recorded during simulation | Bits | Calculated from the agent population's strategy distribution at saved time steps | **Data Format and Missing Values** The CSV file uses comma, as the delimiter with appropriate headers for variables. No missing values are in the dataset; every time point records entropy values for all population sizes simulated. Shannon entropy is defined as 0 if the population is fully homogeneous (all cooperators or all defectors). **Additional Notes** The entropy values provide a measure of diversity and disorder in the strategic composition of populations during the evolutionary simulation. Larger population sizes tend to exhibit prolonged periods of higher entropy, indicating sustained coexistence of cooperative and defective strategies, aligning with the quasi-stable phase analysis. The simulations employ multiprocessing (multiprocessing.Pool) to parallelize independent episodes, improving computational efficiency. #### File: For_figure_5-Phase_space_lattice.py **Description:**  **The output files that will be generated upon running this code**: **potential_landscape_data.csv**\ Contains the empirical potential landscape data representing the system’s effective potential, derived from the stationary distribution of the cooperator fraction across different population sizes on a toroidal lattice network. **Columns:** ·        Population Size (Integer): Number of agents simulated (e.g., 1225, 2025, 3025, 4225, 5625). ·        Center (Float): Value of the cooperator fraction x at which the potential is evaluated. Range approximately between 0 and 1. ·        Potential (Float): Normalized effective potential V(x)=−lnP(x), where P(x) is the estimated probability density from kernel density estimation (KDE). **well_and_barrier_data.csv**\ Summarizes extracted features from the potential landscape for each population size, quantitatively describing system stability. **Columns:** ·        Population Size (Integer): Number of agents as above. ·        Well Depth (Float): Depth of the potential well, calculated as the difference between maximum and minimum potential values within the landscape, measuring the stability strength of the quasi-equilibrium state. ·        Barrier Height (Float): Maximum potential height within the landscape, indicating the energy barrier to transition away from metastable states. ·        Well Position (Float): Cooperator fraction x*x* at the potential well minimum, representing the most probable stable state of cooperation. **potential_landscape_vs_population.png**\ Visualization of the potential landscape V(x) vs. cooperator fraction x for different population sizes. Illustrates how the shape and depth of potential wells evolve as population size increases, revealing stronger confinement of cooperation in larger populations. **well_and_barrier_vs_population.png**\ Plot showing how well depth and barrier height metrics scale with population size, reflecting increased stability and resistance to perturbations in larger networks. **Variables and Abbreviations** | **Variable Name** | **Description** | **Units / Type** | **Comments** | | :---------------- | :--------------------------------------------------------------------------- | :--------------- | :----------------------------------------------------------- | | Population Size | The number of agents in the lattice network simulation | Integer | E.g., 1225, 2025, 3025, 4225, 5625 | | Center | Points in the cooperator fraction x domain where KDE is evaluated | Dimensionless  | Represents fraction of cooperators in the population | | Potential | Effective potential V(x)=−lnP(x) estimated from KDE | Dimensionless | Normalized so the minimum potential is zero | | Well Depth | Depth of the potential well (max potential value minus minimum) | Dimensionless | Larger values indicate stronger metastability | | Barrier Height | Maximum value of the potential landscape | Dimensionless | Represents an energy barrier to leave the quasi-stable state | | Well Position | Cooperator fraction x at which the potential function attains its minimum | Dimensionless  | Most probable state of the population | | timesteps | Number of discrete time points simulated within each episode | Integer | Default: 100,000 | | num\_episodes | Number of independent simulation repeats for averaging | Integer | Default: 100 | | theta | Cooperation enhancement parameter influencing environmental feedback | Dimensionless | Constant: 2.0 | | epsilon | Rate constant for environmental resource change | Dimensionless | Constant: 0.1 | | beta | Selection intensity parameter (inverse temperature) in the Fermi update rule | Dimensionless | Used with a value of 10 in simulations | **Description and Methodology** The potential landscape data results from aggregating long time series data of the cooperator fraction x over many independent simulation episodes for each population size. A kernel density estimation (KDE) with a fixed bandwidth (default: 0.1) smooths the empirical distribution, from which the effective potential V(x)=−lnP(x) is computed and normalized. The landscape informs about the system’s macroscopic stability: deep, narrow wells indicate strong quasi-equilibrium states where the system tends to remain for extended durations. Well depth and barrier height metrics quantify the resilience and stability of cooperation within the network. The code simulates agents on a 2D toroidal lattice with periodic boundary conditions, where each agent updates its strategy based on local interactions and environment-dependent payoffs evolving over time. **Data Format and Missing Values** CSV files use comma separators with headers as described above. The datasets contain no missing values; each row corresponds to a complete record of computed variables for a given population size and cooperator fraction. Numeric values are always present with no blank cells or "NA" entries. #### File: For_figure_6_-final_mutifractal_analysis.py **Description:** Main script implementing agent-based evolutionary game simulations on various network topologies and performing multifractal detrended fluctuation analysis (MF-DFA). Captures cooperator fraction time series and related metrics. **The output files that will be generated upon running this code:** | **Filename / Folder** | **Type** | **Description** | | :--------------------------------------- | :----------- | :-------------------------------------------------------------------------------------------------------------------------------------------- | | multifractal\_spectrum.csv (for lattice) | CSV | Stores the singularity spectrum (α, f(α)) from MF-DFA of cooperator fraction time series for different networks. | | tau\_q\_vs\_q.csv (for lattice) | CSV | Time-dependent mass exponents (tau(q)) vs generalized q values from MF-DFA for cooperative fraction time series. | | h\_2\_values\_Lattice.csv  | CSV | Time evolution of Hurst exponent h(2) values representing long-range correlations in time series for various networks. | | simulation\_results.csv (for lattice) | CSV | Summary of simulation outputs: average cooperative fraction, resource level, and dilemma strengths (Dg, Dr) over time. | | h\_q\_vs\_q\_original\_tXXXX.csv  | CSV | Intermediate MF-DFA results storing generalized Hurst exponents h(q) for original, shuffled, and phase-randomized series at given time steps. | | plots/ | *(optional)* | Folder for generated figures and plots (if saved). Not present by default in code, but suggested for organization. | **Variables and Definitions** | **Variable Name** | **Type** | **Description** | **Units** | | :---------------- | :------- | :------------------------------------------------------------------------ | :-------- | | R0, S0, T0, P0 | float | Payoff matrix parameters for a depleted environment | Unitless | | R1, S1, T1, P1 | float | Payoff matrix parameters for a replete environment | Unitless | | theta | float | Feedback parameter controlling the cooperator's impact on the environment | Unitless | | epsilon | float | Scaling factor for environmental resource dynamics | Unitless | | beta | float | Inverse temperature parameter controlling Fermi update selection strength | Unitless | **Key Classes and Their Variables** | **Class Name** | **Attribute / Method** | **Description** | **Type / Unit** | | :------------- | :----------------------- | :-------------------------------------------------------------------------------------------------- | :---------------------- | | Agent | id | Unique agent identifier | Integer | | 
 
 | point | Payoff accumulated by the agent | Float (unitless) | | 
 
 | strategy | Current strategy: "C" for
 cooperator, "D" for defector | String | | 
 
 | next\_strategy | Strategy to be adopted in the next update step |
 String | | 
 
 | neighbors\_id | List of neighboring agent IDs | List\[int] | | Society
 | size | Number of agents in the population | Integer | | 
 
 | network\_type | Network topology type
 ('lattice', 'ER', 'BA', 'WS') |
 String | | 
 
 | topology | NetworkX graph representing the agent interaction network |
 networkx.Graph | | 
 
 | agents | List of Agent instances | List\[Agent] | | 
 
 | count\_fraction() | Returns the fraction of cooperating agents
 (strategy == "C") | Float (0 to 1) | | Decision | kappa | Noise
 parameter in Fermi update rule (kappa = 1/beta) | Float (unitless) | | 
 
 | Dg, Dr | Gamble-Intending Dilemma strength (Dg) and Risk-Averting
 Dilemma strength (Dr), dynamically updated | Float (unitless) | | 
 
 | count\_payoff(agents) | Computes payoff for each agent based on current
 strategies | Updates point attribute | | 
 
 | pw\_fermi(agents) | Applies Fermi rule to update next\_strategy of
 agents | Updates next\_strategy | | 
 
 | update\_strategy(agents) | Updates strategy from next\_strategy after
 payoff calculation | Updates strategy |   **Key Functions** | **Function
 Name** | **Arguments** | **Returns / Effects** | **Description** | |
 :--------------------------------------------- | :-----------------------
 | :------------------------------------------------------------ |
 :----------------------------------------------------------------- | |
 A(m) | m (float) | 2×2 numpy array | Computes current payoff matrix based
 on resource level m | | calculate\_dilemma\_strengths(A\_evo) | 2×2 matrix
 | Tuple (Dg\_prime, Dr\_prime) | Calculates dilemma strengths based on
 current payoff matrix | | choose\_init\_c(num\_agent, initial\_fraction) |
 int, float | List of agent IDs selected as cooperators | Randomly selects
 initial cooperators | | init\_strategy(agents, init\_c) | List, List |
 None (modifies agent strategies) | Initializes agent strategies
 to "C" or "D" | | shuffle\_time\_series(time\_series)
 | numpy array | Shuffled numpy array | Randomizes time series to destroy
 temporal correlations | | phase\_randomized\_surrogate(time\_series) |
 numpy array | Surrogate time series with randomized phase (FFT domain) |
 Preserves linear properties; destroys nonlinear temporal structure | |
 mdfa(time\_series, q\_values, scales, order) | numpy array, arrays/ints |
 Tuple with MF-DFA outputs – F\_q, h\_q, t\_q, alpha, f\_alpha | Performs
 Multifractal Detrended Fluctuation Analysis | **Simulation Parameters and
 Output Variables** | **Variable Name** | **Type** | **Description** |
 **Units** | | :----------------------------------------------- |
 :------------- |
 :----------------------------------------------------------------------- |
 :------------------------ | | num\_agent | int | Number of agents
 (population size), typically 10,000 | Agents count | | num\_episodes | int
 | Number of simulation runs | Count | | initial\_fraction | float |
 Initial fraction of cooperators (0 ≤ initial\_fraction ≤ 1) | Unitless | |
 t | numpy array | Time points for simulation, typically 0 to 5000 discrete
 steps | Time unit (arbitrary) | | q\_values | numpy array | Generalized
 moment orders q for MF-DFA, e.g., from -5 to 5 | Unitless | | scales |
 numpy array | Segment/window sizes for detrending in MF-DFA | Discrete
 length (samples) | | episode\_cooperator\_fractions | list of floats |
 Time series of cooperation fraction per episode | Fraction (0 to 1) | |
 episode\_resource\_levels | list of floats | Resource level time series
 per episode | Unitless (range 0 to 1) | | tau\_q\_values | list of arrays
 | Time series of mass exponent τ(q) values computed every 100 time steps |
 Unitless | | h\_2\_values | list of floats | Time series of Hurst exponent
 h(2) (long-range correlation measure) | Unitless | |
 alpha\_values, f\_alpha\_values | list of arrays | Singularity spectrum α
 and fractal dimension f(α) at the last time point | Unitless | |
 average\_cooperator\_fraction\_over\_time | numpy array | Average
 cooperative fraction across episodes at each time | Fraction | |
 average\_resource\_levels\_over\_time | numpy array | Average resource
 level across episodes at each time | Unitless | |
 average\_Dg\_over\_time, average\_Dr\_over\_time | numpy array | Average
 dilemma strengths (GID and RAD) across episodes at each time | Unitless |
 **Abbreviations & Units Summary** | **Abbreviation** | **Meaning** | |
 :--------------- |
 :----------------------------------------------------------------------- |
 | C | Cooperator (strategy) | | D | Defector (strategy) | | PD |
 Prisoner’s Dilemma (game type) | | MF-DFA | Multifractal Detrended
 Fluctuation Analysis (time series analysis) | | GID (Dg) |
 Gamble-Intending Dilemma strength (temptation to defect) | | RAD (Dr) |
 Risk-Averting Dilemma strength (fear of being exploited) | | FFT | Fast
 Fourier Transform (for phase randomized surrogate generation) | | q |
 Moment order parameter in MF-DFA | | α and f(α) | Lipschitz-Holder
 exponent and fractal dimension in multifractal spectrum | | h(q) |
 Generalized Hurst exponent (scaling exponent) | | τ(q) | Holder
 exponent/Mass exponent in multifractal analysis | **Missing Values** The
 code and outputs do not explicitly handle missing values. Numerical time
 series and MF-DFA outputs are fully computed; incomplete runs might lead
 to: Blank cells or missing rows in CSV files if outputs are truncated
 during writing. "NA" or "NaN" values are not used in
 current code. Users are advised to verify file integrity and completeness
 accordingly. #### File:
 Calculating_avg_Spectrum_strength_For_Network_Lattice.py
 **Description:** Main script for evolutionary game simulations,
 environment-dependent payoffs, and MF-DFA on a toroidal lattice network.
 Includes 100 simulation episodes executed with multifractal spectrum
 strength analysis to calculate the average spectrum strength. This script
 is just an extended version of the previous script **(For figure 6
 -final_mutifractal_analysis.py)** from multifractal analysis. **The output
 files that will be generated upon running this code:** | **Filename /
 Folder** | **Type** | **Description** | |
 :------------------------------------------------- | :--------- |
 :--------------------------------------------------------------------------------------------------------------------------------------------------- | | mf\_Lattice\_spectrum\_strength\_all\_episodes.csv | CSV | Contains the calculated multifractal spectrum strength (Δα) for each of the 100 simulation episodes for the lattice network. | | multifractal\_Lattice\_spectrum.csv | CSV | Contains multifractal spectrum (α, f(α)) for the last simulation episode in the lattice network case. | | tau\_q\_vs\_q.csv | CSV | Holder exponent / Mass exponent τ(q) vs q values computed at regular intervals (every 100 time steps) aggregated from all episodes. | | h\_2\_values\_Lattice.csv | CSV | Hurst exponent h(2) values representing long-range correlations computed at every 100 time steps across episodes. | | simulation\_Lattice\_results.csv | CSV | Summary of averaged simulation results over episodes: time, average cooperator fraction, resource level, and dilemma strengths Dg and Dr. | | multifractal\_Lattce\_spectrum.png | Image plot | Visualization of the multifractal spectrum f(α) vs α for the last simulation episode in lattice topology. | | surrogation\_Lattice.png | Image plot | Plot comparing generalized Hurst exponents h(q) between original and surrogate (shuffled/phase randomized) data for testing multifractality origins. | **Key Variables and Data Descriptions** | Variable Name | Type | Description | Units / Notes | | :----------------------------- | :----------------- | :------------------------------------------------------------------------------------------------------------- | :------------------------------------ | | num\_agent | int | Number of agents in population (10,000 for lattice simulations) | Count of agents | | num\_episodes | int | Total number of independent simulation runs (100 episodes) | Count | | initial\_fraction | float | Initial proportion of cooperating agents (0.5) | Fraction (range 0 to 1) | | t | np.ndarray | Time points for simulation (linearly spaced over 5000 steps) | Abstract time units | | q\_values | np.ndarray | Array of generalized moment orders q for MF-DFA analysis (-5 to 5 in 21 increments) | Unitless | | scales | np.ndarray | Array of segment sizes (scales) used in MF-DFA averaging (logarithmically spaced) | Discrete time/sample lengths | | order | int | Polynomial detrending order in MF-DFA (linear detrending = 1) | Integer | | episode\_cooperator\_fractions | list of float | Time series of cooperator fractions for each episode | Fraction (0 to 1) | | episode\_resource\_levels | list of float | Time series of resource levels for each episode | Dimensionless resource level (bound) | | alpha\_values | list of np.ndarray | Lipschitz-Holder exponent arrays defining singularity strengths for multifractal spectra at the last time step | Unitless | | f\_alpha\_values | list of np.ndarray | Corresponding multifractal fractal dimensions f(α) at the last time step | Unitless | | tau\_q\_values | list of np.ndarray | Holder exponent / Mass exponents τ(q) are computed every 100 steps for each episode | Unitless (scaling exponents) | | h\_2\_values | list of floats | Hurst exponent values at q=2 indicate long-range correlation strength | Unitless | | mf\_strengths | list of float | Multifractal spectrum strengths Δα, calculated as max(α) - min(α) per episode across f(α)&gt;0 region | Unitless (measure of multifractality) | **Classes and Methods (Extensions Used)** The previously described classes (Agent, Society, Decision) and core methods remain unchanged. **Extended code introduces:** * Parallelized execution using Python multiprocessing.Pool to run 100 independent simulation episodes in parallel (function run_episode). * Within each run_episode, MF-DFA is computed at every 100 time steps on the cooperator fraction time series; results for h(q), τ(q), α, and f(α) are collected. * Calculation of multifractal spectrum strength Δα=max(α)−min(α), filtering f(α)&gt;0 as valid region. * Averaged statistics computed over all episodes for cooperator fraction, resource levels, and dilemma strengths (Dg, Dr). **File Format Details and Data Columns** | **File Name** | **Columns / Variables** | **Description / Notes** | | :------------------------------------------------- | :---------------------------------------------------------------------- | :---------------------------------------------- | | mf\_Lattice\_spectrum\_strength\_all\_episodes.csv | episode (int), mf\_spectrum\_strength (float) | Spectrum strength Δα per episode | | multifractal\_Lattice\_spectrum.csv | alpha (float array), f\_alpha (float array) | Detailed multifractal spectrum for last episode | | tau\_q\_vs\_q.csv | Time\_Step (int), q\_0...q\_20 (float), tau\_q\_0...tau\_q\_20 (float) | MF-DFA mass exponents τ(q) per q and time step | | h\_2\_values\_Lattice.csv | Time\_Step (int), h\_2 (float) | Hurst exponent at q=2 over time | | simulation\_Lattice\_results.csv | Time, Avg\_Cooperator\_Fraction, Avg\_Resource\_Level, Avg\_Dg, Avg\_Dr | Ensemble average simulation outcomes | **Abbreviations** The previously described abbreviations in the main multifractal analysis script **Handling Missing Values** * The simulation and MF-DFA computations do not explicitly handle missing values. * Output CSV files are fully populated for each simulation run. * Missing or incomplete data may occur only if a run is interrupted or an error stops the generation of specific output files. * In such cases, missing cells in CSVs would appear as blank cells. * No usage of "NA" or "NaN" strings for missing values is present in the current outputs. **Notes on Usage and Reproducibility** * Parallel execution utilizes all available CPU cores (up to 48 in this code) to speed up multiple simulation episodes. * MF-DFA is applied consistently at fixed intervals (every 100 steps). * Surrogate time series analyses (shuffling and phase randomization) serve to validate the multifractal origin but are not parallelized inside this code. * Data produced facilitates examination of how multifractal properties evolve during feedback-evolving cooperation dynamics under environmental feedback on lattice networks. * Similar structured analysis is performed analogously for other network types (ER, BA, WS) as described previously. #### File: Calculating_avg_Spectrum_strength_For_Network_BA.py **Description:** Same as File: Calculating_avg_Spectrum_strength_For_Network_Lattice.py, but for the BA network. #### File: Calculating_avg_Spectrum_strength_For_Network_ER.py **Description:** Same as File: Calculating_avg_Spectrum_strength_For_Network_Lattice.py, but for the ER network. #### File: Calculating_avg_Spectrum_strength_For_Network_WS.py **Description:** Same as File: Calculating_avg_Spectrum_strength_For_Network_Lattice.py, but for the WS network. ## Code/software Software Environment The analysis and simulations for the evolutionary games with environmental feedback, along with multifractal detrended fluctuation analysis (MF-DFA), are implemented entirely in Python. All data files included in this submission are generated by Python scripts performing agent-based simulations on various network topologies, followed by detailed multifractal analysis of cooperative fraction time series. The code was developed and tested using Python versions 3.8 and 3.10, but compatibility with Python 3.7 and newer is expected. **Required Python Packages** The following open-source Python packages are essential for running the simulation and analysis scripts: numpy (≥ 1.18): Numerical computations and array operations. pandas (≥ 1.0): Data organization, manipulation, and CSV input/output. networkx (≥ 2.5): Construction and management of complex networks including lattices, Erdős-Rényi, Watts-Strogatz, and Barabási-Albert graphs. matplotlib (≥ 3.1): Generation of all plots and visualizations. scipy (≥ 1.4), specifically scipy.signal: Kernel density estimation (KDE), detrending for MF-DFA, eigenvalue computation, and related operations. tqdm (≥ 4.62): Progress bars for long-running simulations (optional for user feedback). random (built-in Python module): Stochastic elements in agent strategy updates and initialization. multiprocessing (built-in Python module): Enables parallel execution of multiple independent simulation episodes to drastically improve computational efficiency (optional but strongly recommended for multi-core systems). functools (built-in): Used for advanced function handling such as partial function application in parallel runs. **Software Versions Used for Development and Testing** Software / Package Python 3.8, 3.10 numpy 1.21.0 pandas  1.3.0 network 2.6 matplotlib 3.4.2 scipy 1.7.1 tqdm 4.62.3 Note: Compatibility with newer versions of these packages is expected but not explicitly verified. **Workflow and Script Overview** **Simulation Scripts:** The primary Python scripts simulate evolutionary game dynamics on various network structures including Toroidal lattices, Erdős-Rényi (ER), Watts-Strogatz (WS), and Barabási-Albert (BA) networks, with typical population sizes ranging up to 10,000 agents. Each simulation models coupled dynamics of agent strategy evolution and environmental resource feedback via dynamically changing payoff matrices. Central classes such as Society, Agent, and Decision manage network generation, agent attributes, and stochastic strategy update rules (using the Fermi pairwise comparison method). **MF-DFA (Multifractal Detrended Fluctuation Analysis):** The code (for figure 6 associated Codes) implements the full MF-DFA procedure, including time series integration, segmentation, local detrending, calculation of fluctuation functions, extraction of scaling exponents (generalized Hurst exponents), and multifractal spectra. Surrogate data methods (time series shuffling and phase randomization) are included to assess and validate the origin of multifractality. **Parallel Execution:** To manage computational demands, simulation episodes (typically 100 or more per population size) are executed in parallel using Python’s multiprocessing.Pool. This ensures efficient use of available CPU cores and produces statistically robust results. **Data Output:** Results such as cooperative fractions, environmental resource levels, dilemma strengths, Shannon entropy, and multifractal metrics (e.g., singularity spectra, mass exponents) are saved in CSV format, enabling easy import to other analysis environments. Additionally, publication-quality PNG figures illustrating temporal dynamics, potential landscapes, multifractal spectra, and surrogate comparisons are generated and saved automatically or otherwise the Origin software is used for getting the publication quality of the figures. **How to Run the Code** **Set up Python environment** Install Python 3.7 or higher along with all required packages. A typical setup via pip is: pip install numpy, pandas, network, matplotlib, scipy, tqdm etc. Alternatively, use a package/environment manager such as conda for easier dependency handling. **Prepare scripts and input data** Ensure all agent-based simulation scripts, MF-DFA analysis scripts, and any surrogate testing routines are located within your working directory. **Run simulations and analyses** Execute the main simulation scripts (e.g., for Figures 1–6) to generate time series data, statistical summaries, and multifractal outputs. Scripts support customization of parameters like network topology, population size, number of episodes, and simulation length. **Optional parallel execution** If using multiprocessing features, confirm sufficient CPU cores are available. Run the scripts with parallelism enabled by default (adjustable in the code). **Review results** Examine generated CSV data files such as: ·        cooperation_fractions.csv (cooperative fraction time series) ·        mean_entropy.csv (Shannon entropy series) ·        potential_landscape_data.csv (empirical potential landscapes) ·        well_and_barrier_data.csv (summary landscape features) ·        MF-DFA outputs (mf_Lattice_spectrum_strength_all_episodes.csv, tau_q_vs_q.csv, h_2_values_Lattice.csv, etc.) ·        And inspect PNG figures for visualization. **Summary** ·        All code is written in standard Python and uses only free, open-source libraries. ·        The modular architecture divides the tasks into simulation, data processing, and multifractal analysis components with clearly defined functions/classes. ·        Parallel computation is incorporated to handle computationally intensive ensemble runs efficiently. ·        Output file formats (CSV and PNG) are widely supported, facilitating downstream analysis or figure generation in other software. ·        This setup ensures reproducibility and extensibility for further research or adaptation to related evolutionary game problems. ## Access information Other publicly accessible locations of the data: * N/A Data was derived from the following sources: * N/A