xjtulyc/agent-based-social

Agent-Based Social Modeling with Mesa

One-line summary: Build and analyze agent-based social models using Mesa: Schelling segregation, opinion dynamics, social contagion (SIR), and network diffusion with parameter sweeps.

---

When to Use This Skill

• When simulating emergent social phenomena from individual-level rules
• When modeling opinion formation, polarization, or information spread
• When implementing Schelling segregation or residential sorting models
• When running SIR/SIS contagion models on social networks
• When conducting parameter sweeps to find critical thresholds
• When visualizing spatial or network dynamics of social agents

Trigger keywords: agent-based model, ABM, Mesa, Schelling segregation, opinion dynamics, social contagion, information diffusion, SIR model, emergence, simulation, complex systems, social simulation

---

Background & Key Concepts

Agent-Based Modeling

ABM simulates individual agents following local rules to study emergent macro-level behavior:

• Agents: Autonomous decision-makers with state and behavior
• Environment: Grid, network, or continuous space
• Schedule: Defines order of agent activation (random, sequential, staged)
• Data collection: Track population-level statistics over time

Schelling Segregation

Agents move if fewer than $T$% of neighbors share their type. Even with low $T$ (30%), high global segregation emerges — demonstrating how individual preferences aggregate into unintended outcomes.

Opinion Dynamics (Bounded Confidence)

Agents update opinions only if neighbors' opinions are within a "confidence bound" $\epsilon$:

$$ x_i(t+1) = x_i(t) + \mu \sum_{j: |x_j - x_i| < \epsilon} (x_j(t) - x_i(t)) $$

This leads to polarization into isolated opinion clusters.

---

Environment Setup

Install Dependencies

pip install mesa>=2.3 numpy>=1.24 pandas>=2.0 networkx>=3.2 matplotlib>=3.7

Verify Installation

import mesa
print(f"Mesa version: {mesa.__version__}")

# Quick test
from mesa import Agent, Model
from mesa.time import RandomActivation

class TestAgent(Agent):
    def step(self): pass

class TestModel(Model):
    def __init__(self):
        super().__init__()
        self.schedule = RandomActivation(self)
        for i in range(10):
            self.schedule.add(TestAgent(i, self))
    def step(self):
        self.schedule.step()

m = TestModel()
m.step()
print("Mesa OK — 10 agents stepped")

---

Core Workflow

Step 1: Schelling Segregation Model

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mesa import Agent, Model
from mesa.space import MultiGrid
from mesa.time import RandomActivation
from mesa.datacollection import DataCollector

# ------------------------------------------------------------------ #
# Schelling (1971) segregation model
# ------------------------------------------------------------------ #

class SchellingAgent(Agent):
    """Residential agent with type (0 or 1) and similarity preference."""

    def __init__(self, unique_id, model, agent_type, homophily):
        super().__init__(unique_id, model)
        self.type = agent_type
        self.homophily = homophily  # Minimum fraction of same-type neighbors
        self.is_happy = False

    def step(self):
        """Check happiness; if unhappy, move to random empty cell."""
        neighbors = self.model.grid.get_neighbors(
            self.pos, moore=True, include_center=False
        )
        same_type = sum(1 for n in neighbors if n.type == self.type)
        total = len(neighbors)

        if total == 0:
            self.is_happy = True
            return

        self.is_happy = (same_type / total) >= self.homophily

        if not self.is_happy:
            # Move to random empty cell
            empty_cells = list(self.model.grid.empties)
            if empty_cells:
                new_pos = self.random.choice(empty_cells)
                self.model.grid.move_agent(self, new_pos)


class SchellingModel(Model):
    """Schelling segregation model on a grid."""

    def __init__(self, width=30, height=30, density=0.85, pct_type1=0.50, homophily=0.30):
        super().__init__()
        self.schedule = RandomActivation(self)
        self.grid = MultiGrid(width, height, torus=True)
        self.happiness = 0.0

        # Place agents
        agent_id = 0
        for x in range(width):
            for y in range(height):
                if self.random.random() < density:
                    a_type = 0 if self.random.random() < pct_type1 else 1
                    agent = SchellingAgent(agent_id, self, a_type, homophily)
                    self.grid.place_agent(agent, (x, y))
                    self.schedule.add(agent)
                    agent_id += 1

        self.datacollector = DataCollector(
            model_reporters={
                "Happiness":     lambda m: self._happiness(),
                "Segregation":   lambda m: self._segregation_index(),
            }
        )

    def _happiness(self):
        agents = self.schedule.agents
        if not agents:
            return 0.0
        return sum(1 for a in agents if a.is_happy) / len(agents)

    def _segregation_index(self):
        """Dissimilarity index: proportion of one type that would need to move."""
        agents = self.schedule.agents
        type0 = [a for a in agents if a.type == 0]
        type1 = [a for a in agents if a.type == 1]
        N0, N1 = len(type0), len(type1)
        if N0 == 0 or N1 == 0:
            return 0.0
        D = 0.0
        for x in range(self.grid.width):
            for y in range(self.grid.height):
                cell = self.grid.get_cell_list_contents([(x, y)])
                n0 = sum(1 for a in cell if a.type == 0)
                n1 = sum(1 for a in cell if a.type == 1)
                D += abs(n0/N0 - n1/N1)
        return D / 2

    def step(self):
        self.datacollector.collect(self)
        self.schedule.step()


# ---- Run simulation -------------------------------------------- #
model = SchellingModel(width=30, height=30, density=0.85, pct_type1=0.50, homophily=0.35)
for _ in range(50):
    model.step()

data = model.datacollector.get_model_vars_dataframe()

print(f"Final happiness:   {data['Happiness'].iloc[-1]:.3f}")
print(f"Final segregation: {data['Segregation'].iloc[-1]:.3f}")

# ---- Visualization --------------------------------------------- #
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Final grid state
grid_array = np.zeros((model.grid.width, model.grid.height))
for agent in model.schedule.agents:
    x, y = agent.pos
    grid_array[x, y] = agent.type + 1  # 0=empty, 1=type0, 2=type1

axes[0].imshow(grid_array.T, cmap='RdBu', vmin=0, vmax=2, aspect='equal')
axes[0].set_title(f"Final Spatial Distribution\n(homophily={0.35})")
axes[0].axis('off')

axes[1].plot(data['Happiness'] * 100, 'g-', linewidth=2)
axes[1].set_xlabel("Time step"); axes[1].set_ylabel("Happiness (%)")
axes[1].set_title("Agent Happiness Over Time"); axes[1].grid(True, alpha=0.3)

axes[2].plot(data['Segregation'], 'r-', linewidth=2)
axes[2].set_xlabel("Time step"); axes[2].set_ylabel("Dissimilarity Index")
axes[2].set_title("Segregation Index Over Time"); axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("schelling_segregation.png", dpi=150)
plt.show()

Step 2: Opinion Dynamics (Bounded Confidence)

import numpy as np
import matplotlib.pyplot as plt
from mesa import Agent, Model
from mesa.time import SimultaneousActivation
from mesa.datacollection import DataCollector

# ------------------------------------------------------------------ #
# Deffuant-Weisbuch bounded confidence model
# ------------------------------------------------------------------ #

class OpinionAgent(Agent):
    """Agent with a continuous opinion in [0, 1]."""

    def __init__(self, unique_id, model, initial_opinion):
        super().__init__(unique_id, model)
        self.opinion = initial_opinion
        self.new_opinion = initial_opinion

    def step(self):
        """Interact with a random neighbor within confidence bound."""
        # Select random agent to interact with
        partner = self.random.choice(self.model.schedule.agents)
        if partner is self:
            return
        diff = abs(self.opinion - partner.opinion)
        if diff < self.model.epsilon:
            mu = self.model.mu  # Convergence parameter
            self.new_opinion = self.opinion + mu * (partner.opinion - self.opinion)

    def advance(self):
        """Apply the update (simultaneous activation)."""
        self.opinion = np.clip(self.new_opinion, 0, 1)


class OpinionModel(Model):
    """Bounded confidence opinion dynamics."""

    def __init__(self, n_agents=500, epsilon=0.3, mu=0.5):
        super().__init__()
        self.epsilon = epsilon
        self.mu = mu
        self.schedule = SimultaneousActivation(self)

        for i in range(n_agents):
            op = self.random.uniform(0, 1)
            agent = OpinionAgent(i, self, op)
            self.schedule.add(agent)

        self.datacollector = DataCollector(
            agent_reporters={"Opinion": "opinion"},
            model_reporters={
                "Mean":    lambda m: np.mean([a.opinion for a in m.schedule.agents]),
                "Std":     lambda m: np.std([a.opinion for a in m.schedule.agents]),
                "N_clusters": lambda m: self._count_clusters(m),
            }
        )

    def _count_clusters(self, m, tol=0.05):
        """Count opinion clusters (groups within tol of each other)."""
        opinions = sorted([a.opinion for a in m.schedule.agents])
        if not opinions:
            return 0
        clusters, cluster_start = 1, opinions[0]
        for op in opinions[1:]:
            if op - cluster_start > tol:
                clusters += 1
                cluster_start = op
        return clusters

    def step(self):
        self.datacollector.collect(self)
        self.schedule.step()


# ---- Compare two epsilon values -------------------------------- #
fig, axes = plt.subplots(2, 2, figsize=(13, 9))

for col, eps in enumerate([0.15, 0.45]):
    model = OpinionModel(n_agents=300, epsilon=eps, mu=0.5)
    snapshots = {0: None, 50: None, 200: None}

    for step in range(201):
        model.step()
        if step in snapshots:
            snapshots[step] = [a.opinion for a in model.schedule.agents]

    # Opinion distribution at different times
    for t, opinions in snapshots.items():
        axes[0][col].hist(opinions, bins=40, range=(0, 1),
                          alpha=0.6, density=True, label=f"t={t}")
    axes[0][col].set_xlabel("Opinion"); axes[0][col].set_ylabel("Density")
    axes[0][col].set_title(f"Opinion Distribution (ε={eps})")
    axes[0][col].legend(fontsize=9); axes[0][col].grid(True, alpha=0.3)

    # Cluster count over time
    cluster_data = model.datacollector.get_model_vars_dataframe()['N_clusters']
    axes[1][col].plot(cluster_data, color='purple', linewidth=1.5)
    axes[1][col].set_xlabel("Time step"); axes[1][col].set_ylabel("Number of clusters")
    axes[1][col].set_title(f"Cluster Formation (ε={eps})"); axes[1][col].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("opinion_dynamics.png", dpi=150)
plt.show()

Step 3: SIR Contagion on Social Network

import numpy as np
import pandas as pd
import networkx as nx
import matplotlib.pyplot as plt
from mesa import Agent, Model
from mesa.time import RandomActivation
from mesa.datacollection import DataCollector

# ------------------------------------------------------------------ #
# SIR contagion model on a small-world network
# ------------------------------------------------------------------ #

class SIRAgent(Agent):
    """S→I→R agent on a contact network."""

    def __init__(self, unique_id, model, state='S'):
        super().__init__(unique_id, model)
        self.state = state
        self.days_infected = 0

    def step(self):
        if self.state == 'I':
            # Infect susceptible neighbors
            neighbors_ids = list(self.model.network.neighbors(self.unique_id))
            for nid in neighbors_ids:
                neighbor = self.model.agents_dict.get(nid)
                if neighbor and neighbor.state == 'S':
                    if self.random.random() < self.model.beta:
                        neighbor.state = 'I'

            # Recover
            self.days_infected += 1
            if self.days_infected >= self.model.gamma_inv:
                self.state = 'R'


class SIRModel(Model):
    """SIR contagion on Watts-Strogatz small-world network."""

    def __init__(self, n_agents=500, k=6, p_rewire=0.1,
                 beta=0.05, gamma_inv=7, seed_frac=0.01):
        super().__init__()
        self.beta     = beta       # Transmission probability per contact per day
        self.gamma_inv = gamma_inv  # Mean infectious period (days)

        # Create Watts-Strogatz network
        self.network = nx.watts_strogatz_graph(n_agents, k, p_rewire, seed=42)
        self.schedule = RandomActivation(self)
        self.agents_dict = {}

        for node in self.network.nodes():
            state = 'I' if self.random.random() < seed_frac else 'S'
            agent = SIRAgent(node, self, state)
            self.schedule.add(agent)
            self.agents_dict[node] = agent

        self.datacollector = DataCollector(
            model_reporters={
                "S": lambda m: sum(1 for a in m.schedule.agents if a.state=='S'),
                "I": lambda m: sum(1 for a in m.schedule.agents if a.state=='I'),
                "R": lambda m: sum(1 for a in m.schedule.agents if a.state=='R'),
            }
        )

    def step(self):
        self.datacollector.collect(self)
        self.schedule.step()
        # Stop when epidemic ends
        return sum(1 for a in self.schedule.agents if a.state=='I') > 0

    @property
    def R0(self):
        """Basic reproduction number: R₀ = β × <k> × γ_inv"""
        avg_degree = np.mean([d for _, d in self.network.degree()])
        return self.beta * avg_degree * self.gamma_inv


# ---- Simulation ------------------------------------------------- #
n_agents = 500
model = SIRModel(n_agents=n_agents, k=6, p_rewire=0.1,
                 beta=0.04, gamma_inv=7, seed_frac=0.02)

print(f"R₀ = {model.R0:.2f}  (epidemic if R₀>1)")

for day in range(100):
    active = model.step()
    if not active and day > 10:
        print(f"Epidemic ended at day {day}")
        break

data = model.datacollector.get_model_vars_dataframe()
final = data.iloc[-1]
print(f"\nFinal state: S={final['S']:.0f} ({final['S']/n_agents*100:.1f}%), "
      f"I={final['I']:.0f}, R={final['R']:.0f} ({final['R']/n_agents*100:.1f}%)")

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(data['S']/n_agents*100, 'b-', label='Susceptible', linewidth=2)
ax.plot(data['I']/n_agents*100, 'r-', label='Infectious',  linewidth=2)
ax.plot(data['R']/n_agents*100, 'g-', label='Recovered',   linewidth=2)
ax.set_xlabel("Day"); ax.set_ylabel("Population fraction (%)")
ax.set_title(f"SIR Epidemic on Watts-Strogatz Network  (R₀={model.R0:.2f})")
ax.legend(); ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("sir_epidemic.png", dpi=150)
plt.show()

---

Advanced Usage

Parameter Sweep with BatchRunner

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# ------------------------------------------------------------------ #
# Sweep homophily threshold in Schelling model
# ------------------------------------------------------------------ #

homophily_values = np.arange(0.1, 0.9, 0.1)
n_runs = 5
results = []

for homophily in homophily_values:
    for run in range(n_runs):
        model = SchellingModel(width=20, height=20, density=0.80,
                               pct_type1=0.50, homophily=float(homophily))
        for _ in range(30):
            model.step()
        data = model.datacollector.get_model_vars_dataframe()
        results.append({
            'homophily': homophily,
            'run': run,
            'segregation': data['Segregation'].iloc[-1],
            'happiness':   data['Happiness'].iloc[-1],
        })

df = pd.DataFrame(results)
summary = df.groupby('homophily')[['segregation','happiness']].agg(['mean','std'])

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
for ax, col, title in zip(axes, ['segregation','happiness'], ['Segregation','Happiness']):
    means = summary[col]['mean']
    stds  = summary[col]['std']
    ax.plot(homophily_values, means, 'bo-', linewidth=2, markersize=7)
    ax.fill_between(homophily_values, means-stds, means+stds, alpha=0.2, color='blue')
    ax.set_xlabel("Homophily threshold"); ax.set_ylabel(title)
    ax.set_title(f"{title} vs. Homophily (Schelling)"); ax.grid(True, alpha=0.3)
plt.tight_layout(); plt.savefig("schelling_sweep.png", dpi=150); plt.show()

---

Troubleshooting

Error: AttributeError: 'NoneType' object has no attribute 'step'

Cause: Agent scheduler not initialized.

Fix:

class MyModel(Model):
    def __init__(self):
        super().__init__()
        self.schedule = RandomActivation(self)  # Must initialize before adding agents

Error: ImportError: cannot import name 'MultiGrid'

Cause: Mesa 2.x moved some classes.

Fix:

# Mesa >= 2.0
from mesa.space import MultiGrid        # Still available
from mesa.time import RandomActivation  # Still available

Slow simulation with large grids

# Use numpy vectorized operations instead of agent loops for large N
# Or reduce grid resolution: width=height=20 runs much faster than 100×100

Version Compatibility

| Package | Tested versions | Notes | |:--------|:----------------|:------| | mesa | 2.3, 2.4 | API stable since 2.0; BatchRunner removed in 2.x — use loops | | networkx | 3.2 | No known issues |

---

External Resources

Official Documentation

• Mesa documentation
• Mesa examples repository

Key Papers

• Schelling, T.C. (1971). Dynamic models of segregation. Journal of Mathematical Sociology.
• Deffuant, G. et al. (2000). Mixing beliefs among interacting agents. Advances in Complex Systems.

---

Examples

Example 1: Wealth Distribution (Sugarscape)

import numpy as np
import matplotlib.pyplot as plt

# Simplified Sugarscape: agents accumulate wealth, consume, and die
np.random.seed(42)
n = 500
wealth = np.random.exponential(10, n)  # Initial wealth
metabolism = np.random.randint(1, 5, n)  # Daily consumption

for t in range(100):
    # Sugar grows (simplified as random income)
    income = np.random.poisson(3, n)
    wealth += income - metabolism
    # Dead agents replaced with wealth=random
    dead = wealth <= 0
    wealth[dead] = np.random.exponential(5, dead.sum())

# Compute Gini
w_sorted = np.sort(wealth)
n = len(w_sorted)
G = (2 * np.sum(np.arange(1,n+1) * w_sorted) / (n * w_sorted.sum())) - (n+1)/n
print(f"Sugarscape Gini coefficient: {G:.4f}")

fig, ax = plt.subplots(figsize=(8, 4))
ax.hist(wealth, bins=50, color='gold', edgecolor='brown', linewidth=0.5)
ax.set_xlabel("Wealth"); ax.set_title(f"Sugarscape Wealth Distribution (Gini={G:.3f})")
ax.grid(alpha=0.3); plt.tight_layout(); plt.savefig("sugarscape.png", dpi=150); plt.show()

Example 2: Network Influence Cascade

import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

G = nx.barabasi_albert_graph(500, 3, seed=42)
# Simple threshold diffusion: adopt if >30% of neighbors adopted
adopted = set(list(G.nodes())[:5])  # Seed with 5 adopters
cascade_size = [len(adopted)]

for _ in range(50):
    new_adopters = set()
    for node in G.nodes():
        if node in adopted:
            continue
        neighbors = list(G.neighbors(node))
        if not neighbors:
            continue
        frac_adopted = sum(1 for n in neighbors if n in adopted) / len(neighbors)
        if frac_adopted > 0.30:
            new_adopters.add(node)
    adopted.update(new_adopters)
    cascade_size.append(len(adopted))
    if not new_adopters:
        break

print(f"Final cascade size: {len(adopted)}/{G.number_of_nodes()} ({len(adopted)/G.number_of_nodes()*100:.1f}%)")
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(cascade_size, 'b-o', markersize=4, linewidth=1.5)
ax.set_xlabel("Time step"); ax.set_ylabel("Adopters"); ax.set_title("Information Cascade on BA Network")
ax.grid(alpha=0.3); plt.tight_layout(); plt.savefig("cascade.png", dpi=150); plt.show()

---

Last updated: 2026-03-17 | Maintainer: @xjtulyc Issues: GitHub Issues