Case Study: Predator-Prey Agent-Based Model - Modelling and Simulation Cookbook

This case study develops a spatial predator–prey model on a two-dimensional lattice. The ecosystem contains two types of agents: Prey, which move and reproduce, and predators, which move, die, consume neighbouring prey and reproduce following successful predation.

Topics:

Agent-Based Modelling

Objectives:

Which difficulties are involved in agent-based models with spatial structure?
Why is it important to update agents simultaneously in models with discrete time steps?
Why is it tricky to remove agents from an agent-based model?

from abc import abstractmethod
from typing import Any, Self

import holoviews as hv
import matplotlib.pyplot as plt
import numpy as np

from modelling_and_simulation_cookbook import (
    frames_to_json_data,
    show_grid_abm_animation,
)

hv.extension("bokeh")

Modelling Purpose¶

In principle, the model in this case study could be motivated analogous to the model in Simulation Circle. I.e. we want to help a farmer decising for an intervention into an existing predator-prey system, however, this time we found that consideration of spatial effects is crucial for the model’s validity. Therefore a microscopic modelling strategy was chosen in favour of the macroscopic System Dynamics approach.

However, the agent-based model(s) derived/used in this case-study are primarily used to teach about features and problems of spatially-discrete agent-based models. Therefore the focus lies more on clarity and computational simplicity than on validity compared with the real system. We will discuss the biological implications of the model(s) in the final section of this notebook.

Conceptual Modelling¶

We define an agent-based model (ABM) on the following modelling assumptions:

Each agent in the model represents either prey or a predator.
All agents inhabit a commin environment in form of a finite rectangular lattice. Each cell can contain no more than one agent.
Agents interact through a spatial neighborhood.
The model is updated in discrete time steps and consists of agent movement, where agents update their positions, and population dynamics, where agents die, reproduce and interact with each other.
Movement is random within the agent’s neighourhood and is not influenced by e.g. food availability, predator avoidance or habitat quality.
Predator have a fixed probability of dying at the end of each timestep.
A predator may consume up to one neighbouring prey agent each timestep with a fixed probability.
Following successful predation, the prey is removed and a new predator is placed in its former location.
Each prey agent may reproduce into an empty neighbouring cell with a fixed probability.
All probabilities remain constant throughout a simulation.

We discuss conceptual details later on together with the implementation.

Model Parameters - Conceptual¶

The model has five parameters: the height $m\in \mathbb{N}$ and width $n\in \mathbb{N}$ of the grid, the initial numbers of prey $N_{prey}\in \mathbb{N}$ and predators $N_{pred}\in \mathbb{N}$ with $m\cdot n<N_{prey}+N_{pred}$ , and three probabilities for predator death $p_d\in [0,1]$ , prey birth $p_b\in [0,1]$ and predation success $p_h\in [0,1]$ .

Model Parameters - Implemented¶

For the implemented version, the actual conceptual model parameters often fall together with other experiment specific parameters. In this case, we must also define the seed for generation of the random numbers. Note that running the notebook with the same software environment and seed produces the same sequence of movement and interaction events. Changing the seed does not modify the model itself, but it generates a different realisation of the same stochastic process.

# Spatial and temporal configuration
GRID_SIZE = (15, 15)
N_PREY_INITIAL = 100
N_PREDATORS_INITIAL = 20

# Transition probabilities per time step
PREDATOR_DEATH_PROB = 0.05
PREDATION_PROB = 0.07
PREY_REPRODUCTION_PROB = 0.20

# Random seed for reproducibility
SEED = 4321

Agents - Conceptual¶

Each model agent represents either a prey or a predator animal. Since the model does not require any specific agent states apart from each position on the grid, this is, conceptually, its only state.

Agents - Implemented¶

In the implemented model, we would extend this by a second attribute for technical reasons:

agent_id: a unique identifier used to track removal, reproduction and conflict resolution.
position: a (row, column) tuple specifying the occupied lattice cell.

Adding a unique agent id helps to identify agents in the implemented model. However, in particular in population-dynamic agent-based models it is not always simple to assign such ids (see later). To differentiate optimally between predator and prey animals, we may use class ineritance to derive a Prey and a Predator class from the same Agent class. Addig a copy method is often a good strategy to help making things reproducible. Also some serialization method like the defined to_dict method is useful to create something object-independent to save.

class Agent:
    def __init__(self, agent_id: int, position: tuple[int, int]):
        """ "
        Consructor of an agent
        :param agent_id: unique id
        :param position: coordinate of the agent on the grid
        """
        self.agent_id = agent_id
        self.position = position

    @abstractmethod
    def copy(self) -> Self:
        """
        Returns a copy of this agent
        """

    def to_dict(self) -> dict[str, Any]:
        """
        Returns a dictionary representation of the agent
        """


class Predator(Agent):
    def copy(self):
        return Predator(self.agent_id, self.position)

    def to_dict(self):
        return {
            "type": "prey",
            "row": int(self.position[0]),
            "col": int(self.position[1]),
        }


class Prey(Agent):
    def copy(self):
        return Prey(self.agent_id, self.position)

    def to_dict(self):
        return {
            "type": "predator",
            "row": int(self.position[0]),
            "col": int(self.position[1]),
        }

Agent Initialisation - Conceptual¶

To initialise the agent population we need to create the defined number of prey and predator agents and assign random unique grid-points to them.

Agent Initialisation - Implemented¶

As seen above, the initial numbers of prey and predator agents are defined by N_PREY_INITIAL and N_PREDATORS_INITIAL and the size of the grid is defined by GRID_SIZE. Since positions must be unique, it is not possible to simply draw a random position for each agent, because this would not prevent different agents to get the same assigned gridpoint. In our implementation the task is solved by rng.choice without replacement together with the divmod function to map the unique 1D position onto a 2D grid. Using shuffling would be another viable strategy.

def initialize_agents(
    grid_size: tuple[int, int], n_prey: int, n_predators: int, rng: np.random.Generator
) -> list[Agent]:
    """
    Create the initial prey and predator populations at unique random positions.

    :param grid_size: Number of rows and columns in the lattice.
    :param n_prey: Initial number of prey agents.
    :param n_predators: Initial number of predator agents.
    :param rng: Random number generator
    :returns: Initial population with no more than one agent per cell.
    """
    n_rows, n_cols = grid_size
    total_agents = n_prey + n_predators
    total_cells = n_rows * n_cols

    if n_prey < 0 or n_predators < 0:
        raise ValueError("Population sizes must be non-negative.")

    if total_agents > total_cells:
        raise ValueError("The number of agents cannot exceed the number of grid cells.")

    selected_cells = rng.choice(
        total_cells, size=total_agents, replace=False
    )  # draw without replacement

    positions = [
        divmod(int(cell_index), n_cols) for cell_index in selected_cells
    ]  # create a 2d index from a scalar one

    agents = []  # create the agent list

    for agent_id, position in enumerate(positions[:n_prey]):
        agents.append(Prey(agent_id, position))

    for offset, position in enumerate(positions[n_prey:]):
        agents.append(Predator(n_prey + offset, position))

    return agents

We furthermore write some utility functions to help with the analysis and visualisation of the agent-population and test our initialisation routine.

def count_agents(agents: list[Agent]) -> tuple[int, int]:
    """
    Return the current numbers of prey and predator agents.
    :param agents: list of agents
    :returns: number of prey and predators
    """
    n_prey = len([agent for agent in agents if isinstance(agent, Prey)])
    n_predators = len([agent for agent in agents if isinstance(agent, Predator)])

    return n_prey, n_predators

def check_unique_positions(agents: list[Agent]) -> bool:
    """
    Checks whether all agents occupy unique positions.
    :param agents: list of agents
    :returns: true, if no gridpoint is occupied twice
    """
    positions = [agent.position for agent in agents]
    return len(positions) == len(set(positions))

# Colour-blind-friendly colours from the Okabe–Ito palette
PREY_COLOR = "#0072B2"  # blue
PREDATOR_COLOR = "#D55E00"  # vermillion


def plot_agents(
    agents: list[Agent],
    grid_size: tuple[int, int],
    title: str,
    show_counts: bool = True,
    show_legend: bool = True,
    wiggle: bool = False,
):
    """
    Creates a figure with current spatial positions of prey and predator agents.

    :param agents: list of agents
    :param grid_size: size of the grid as integer pair
    :param title: title for the image
    :param show_counts: if true, the routine will add the counts for predator and prey
        to the title
    :param show_legend: if true, the routine will display a legend to the right of the
        plot
    :param wiggle: if true, the method will add a random noise to the positions of the
        agents (so that it would be visible if agents do not have a unique position)
    """
    n_rows, n_cols = grid_size

    if show_counts:
        n_prey, n_predators = count_agents(agents)
        title = f"{title}\nPrey: {n_prey} | Predators: {n_predators}"
    fig = plt.gcf()
    ax = plt.gca()
    ax.set(
        xlim=(0, n_cols),
        ylim=(0, n_rows),
        aspect="equal",
        title=title,
        xticks=np.arange(n_cols + 1),
        xticklabels=["" for x in np.arange(n_cols + 1)],
        yticks=np.arange(n_rows + 1),
        yticklabels=["" for x in np.arange(n_rows + 1)],
    )
    ax.grid(True, color="black", linewidth=0.35, alpha=0.2)

    # determine how large to make the markers
    plt.gcf().canvas.draw()
    p0 = ax.transData.transform((0, 0))
    p1 = ax.transData.transform((0.5, 0))
    dx_pixels = p1[0] - p0[0]
    marker_size = (dx_pixels * 72 / fig.dpi) ** 2

    data = {
        "Prey": {
            "positions": np.array(
                [a.position for a in agents if isinstance(a, Prey)], dtype=np.float64
            ),
            "color": PREY_COLOR,
            "marker": "o",
        },
        "Predators": {
            "positions": np.array(
                [a.position for a in agents if isinstance(a, Predator)],
                dtype=np.float64,
            ),
            "color": PREDATOR_COLOR,
            "marker": "^",
        },
    }
    if wiggle:
        rng = np.random.default_rng(SEED)
        # make sure to create an own random number generator for this
        # to avoid any interactions with the simulation
        data["Prey"]["positions"] += 0.3 * (
            rng.random([len(data["Prey"]["positions"]), 2]) - 0.5
        )
        data["Predators"]["positions"] += 0.3 * (
            rng.random([len(data["Predators"]["positions"]), 2]) - 0.5
        )

    for label, group in data.items():
        positions = group["positions"]

        if len(positions) > 0:
            rows, cols = positions[:, 0], positions[:, 1]

            # Centre each marker within its cell and display row 0 at the top.
            ax.scatter(
                cols + 0.5,
                n_rows - rows - 0.5,
                s=marker_size,
                color=group["color"],
                marker=group["marker"],
                edgecolor="black",
                linewidth=0.4,
                label=label,
                zorder=3,
            )
    if show_legend:
        ax.legend(loc="center left", bbox_to_anchor=(1.03, 0.5), frameon=True)

# let's test out initialisation strategy
rng = np.random.default_rng(SEED)
agents = initialize_agents(GRID_SIZE, N_PREY_INITIAL, N_PREDATORS_INITIAL, rng)

print(f"Initial number of agents: {len(agents)}")
print(f"Initial nmber of prey and predators: {count_agents(agents)}")
print("All agents occupy unique positions:", check_unique_positions(agents))

plt.figure(figsize=(4, 4))
plot_agents(agents, GRID_SIZE, title="Initial Spatial Distribution")
plt.show()

Update Rules - Conceptual¶

As mentioned, we update the model in discrete time-step, whereas each time-step should be divided into two sequential phases:

Movement phase Agents (may) move to an empty neighbouring cell
Population-dynamics phase Predators may die or consume neighbouring prey, while prey may reproduce into an empty neighbouring cell.

Hereby, phase 2 is started as soon as phase 1 is finished for all agents.

Neighbourhood - Conceptual¶

As seen, the term neighbourhood plays a central role for the model update. In out model we will use the so called Moore neighbourhood, defined by the eight cells that share either an edge or a corner with the focal cell Moore (1962). At the edge of the lattice, torodorial boundary conditions are applied in both spatial dimensions. That means, a position extending beyond one boundary re-enters through the opposite boundary. On a two-dimensional rectangular lattice, connecting the left and right boundaries and the upper and lower boundaries produces a toroidal topology (donut) which gives these boundary conditions their name Chopard & Droz (1998),Wilensky (1999). This choice ensures that boundary cells have the same number of potential neighbours as interior cells. It therefore avoids the systematic reduction in interaction opportunities that would occur at open edges, although the toroidal geometry remains a modelling simplification rather than a literal ecological representation.

Neighbourhood - Implemented¶

Typically, the Moore neighourhood is implemented using the eight relative offsets $\{(i,j):i\in \{-1,0,1\},j\in \{-1,0,1\},i\neq 0\vee j\neq 0\}$ . Note that the target cell itself is not counted to its own neighbourhood (therefore $i\neq 0\vee j\neq 0$ ). Implementation of torodorial neighbouhood is particularly easy, because indices can be wrapped using the remainder operator (%) with the number of rows and columns. Note that this computation is marginally more complicated in programming languages in which indexing starts with 1 and not with 0

Periodic boundary conditions are applied in both spatial dimensions. The row and column indices are wrapped using the remainder operator (%), so a position extending beyond one boundary re-enters through the opposite boundary. On a two-dimensional rectangular lattice, connecting the left and right boundaries and the upper and lower boundaries produces a toroidal topology Chopard & Droz (1998),Wilensky (1999). This choice ensures that boundary cells have the same number of potential neighbours as interior cells. It therefore avoids the systematic reduction in interaction opportunities that would occur at open edges, although the toroidal geometry remains a modelling simplification rather than a literal ecological representation.

def get_neighbour_positions(
    position: tuple[int, int], grid_size: tuple[int, int]
) -> set[tuple[int, int]]:
    """
    Return a list of indices in the Moore neighbourhood around a position.
    Periodic boundary conditions are applied.
    :param position: target position
    :param grid_size: size of the grid
    :returns: set with positions (8 elements)
    """
    row, col = position
    n_rows, n_cols = grid_size
    assert n_rows > 2 and n_cols > 2, ValueError(
        "Grid must have at least three cols/rows"
    )

    neighbour_positions = set()

    for d_row in [-1, 0, 1]:
        for d_col in [-1, 0, 1]:
            if d_row == 0 and d_col == 0:
                continue
            neighbour = ((row + d_row) % n_rows, (col + d_col) % n_cols)
            neighbour_positions.add(neighbour)
    return neighbour_positions

Movement Phase - Conceptual¶

In the movement phase, every agent should randomly move to an empty neighbouring cell.

Movement Phase - Implementation¶

Although the concept sounds well defined, surprisingly, it is not unproblematic. To show this we create several versions of the implementation.

Movement Phase - Implementation Version 1¶

In this first version, we simply update each agent’s position by randomly selecting an empty place in its neighbourhood. Note that this requires us to search for an empty site. The naive implementation would be

for agent in agents:
  pos = agent.position
  neighbour_positions = get_neighbor_positions(pos)
  for agent2 in agents:
      if agent2.position in neighbour_positions:
          neighbour_positions.pop(agent2.position)

However, since this creates a $\mathcal{O}(N^2)$ computational effort, we first create a mapping between positions and agents.

def create_position_agent_mapping(agents: list[Agent]) -> dict[tuple[int, int], Agent]:
    """
    Creates a mapping between agents and positions
    Only works properly if every agent at most occupies one site.
    :param agents: list of agents
    :returns: dictionary mapping position tuples onto agents
    """
    mapping = dict()
    for agent in agents:
        mapping[agent.position] = agent
    return mapping


def get_empty_neigbour_positions(
    position: tuple[int, int],
    grid_size: tuple[int, int],
    position_agent_mapping: dict[tuple[int, int], Agent],
) -> list[tuple[int, int]]:
    """
    Returns all neighbour positions, which are currently not occupied by an agent
    :param position: target position
    :param grid_size: size of the grid
    :param position_agent_mapping: dictionary mapping position tuples onto agents
    :returns: set with positions (maximum 8 elements)
    """
    nbs = get_neighbour_positions(position, grid_size)
    empty = list()
    for nb in nbs:
        if nb not in position_agent_mapping:
            empty.append(nb)
    return empty


def movement_phase_v1(
    agents: list[Agent], grid_size: tuple[int, int], rng: np.random.Generator
) -> None:
    """
    Move agents randomly by randomly selecting a empty position in the Moore
    neighbourhood
    :param agents: list of agents
    :param grid_site: size of the grid
    :param rng: random number generator
    """
    position_agent_mapping = create_position_agent_mapping(agents)

    for agent in agents:
        possible_positions = get_empty_neigbour_positions(
            agent.position, grid_size, position_agent_mapping
        )
        if len(possible_positions) > 0:
            i = rng.integers(0, len(possible_positions))
            new_position = possible_positions[i]
            agent.position = new_position
        else:
            pass  # agent must stay where it is


plt.figure(figsize=(12, 5))
rng = np.random.default_rng(SEED)
agents = initialize_agents(GRID_SIZE, N_PREY_INITIAL, N_PREDATORS_INITIAL, rng)
plt.subplot(1, 3, 1)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 0, unique positions: {check_unique_positions(agents)}",
    wiggle=True,
    show_legend=False,
)
plt.subplot(1, 3, 2)
movement_phase_v1(agents, GRID_SIZE, rng)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 1, unique positions: {check_unique_positions(agents)}",
    wiggle=True,
    show_legend=False,
)
plt.subplot(1, 3, 3)
movement_phase_v1(agents, GRID_SIZE, rng)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 2, unique positions: {check_unique_positions(agents)}",
    wiggle=True,
)
plt.show()

Just by looking at the pictures, one discovers very obvious problems with the implementation: as all agents move independently of one another, it is possible that two or more of them may, by chance, head for the same location. This violates the idea that every site can only be occupied by a single agent.

Movement Phase - Implementation Version 2¶

We try to resolve this problem by updating the position_agent_mapping after every single agent movement.

def movement_phase_v2(
    agents: list[Agent], grid_size: tuple[int, int], rng: np.random.Generator
) -> None:
    """
    Move agents randomly by randomly selecting a empty position in the Moore
    neighbourhood.
    Empty spaces are updated everytime an agent moves
    :param agents: list of agents
    :param grid_site: size of the grid
    :param rng: random number generator
    """
    position_agent_mapping = create_position_agent_mapping(agents)

    for agent in agents:
        possible_positions = get_empty_neigbour_positions(
            agent.position, grid_size, position_agent_mapping
        )
        if len(possible_positions) > 0:
            i = rng.integers(0, len(possible_positions))
            new_position = possible_positions[i]
            position_agent_mapping.pop(agent.position)  # remove the original position
            agent.position = new_position
            position_agent_mapping[agent.position] = (
                agent  # add the agent to the new position
            )
        else:
            pass  # agent must stay where it is


plt.figure(figsize=(12, 5))
rng = np.random.default_rng(SEED)
agents = initialize_agents(GRID_SIZE, N_PREY_INITIAL, N_PREDATORS_INITIAL, rng)
plt.subplot(1, 3, 1)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 0, unique positions: {check_unique_positions(agents)}",
    show_legend=False,
)
plt.subplot(1, 3, 2)
movement_phase_v2(agents, GRID_SIZE, rng)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 1, unique positions: {check_unique_positions(agents)}",
    show_legend=False,
)
plt.subplot(1, 3, 3)
movement_phase_v2(agents, GRID_SIZE, rng)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 2, unique positions: {check_unique_positions(agents)}",
)
plt.show()

This solution now ensures that the location assignment is unique, but it gives rise to two further problems:

First, the order in which the agents are called leads to a bias in the model. The agent with index 0 can always carry out its movement first, and the one with index N always last. Although this does not yet result in any obvious advantage or disadvantage for individual agents in the current model (which involves only movement and no further interaction), it can quickly lead to such a situation when combined with other rules.
The second problem is that the implementation violates a fundamental law of discrete-time agent-based modelling: simultaneous update. The update for agent 1 takes into account the already updated state of agent 0, but also the states of agents 2 to N, which have not yet been updated.

Movement Phase - Implementation Version 3¶

An obvious solution to the first problem is to iterate through the agents in a randomised order. In many cases this approach is often used in ABMs to rule out systematic bias. However, strictly speaking, this is not a formally correct time-discrete implementation, because the second problem remains.

A correct solution that resolves both problems is actually not that easy to implement. One approach involves a two-stage iteration: the first iteration is carried out analogous to Version 1, but the randomly determined new positions are not accepted immediately. Instead, they are collected as requests. In the second iteration, the requests are resolved by randomly selecting one agent from among those requesting the same field. This agent is allowed to move whereas the other ones must remain where they are.

def movement_phase_v3(
    agents: list[Agent], grid_size: tuple[int, int], rng: np.random.Generator
) -> None:
    """
    Move agents randomly by randomly selecting a empty position in the Moore
    neighbourhood
    Target positions are first requested and afterwards evaluated.
    :param agents: list of agents
    :param grid_site: size of the grid
    :param rng: random number generator
    """
    position_agent_mapping = create_position_agent_mapping(agents)

    requests = dict()
    for agent in agents:
        possible_positions = get_empty_neigbour_positions(
            agent.position, grid_size, position_agent_mapping
        )
        if len(possible_positions) > 0:
            i = rng.integers(0, len(possible_positions))
            requested_position = possible_positions[i]
            if requested_position not in requests:
                requests[requested_position] = list()
            requests[requested_position].append(
                agent
            )  # the agent henceforth request the position
        else:
            pass  # agent must stay where it is

    for position, requestees in requests.items():
        if len(requestees) == 1:
            # if there is only one agent who wants to move there, let it move
            requestees[0].position = position
        else:
            # otherwise, pick one at random
            i = rng.integers(0, len(requestees))
            requestees[i].position = position


plt.figure(figsize=(12, 5))
rng = np.random.default_rng(SEED)
agents = initialize_agents(GRID_SIZE, N_PREY_INITIAL, N_PREDATORS_INITIAL, rng)
plt.subplot(1, 3, 1)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 0, unique positions: {check_unique_positions(agents)}",
    show_legend=False,
)
plt.subplot(1, 3, 2)
movement_phase_v3(agents, GRID_SIZE, rng)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 1, unique positions: {check_unique_positions(agents)}",
    show_legend=False,
)
plt.subplot(1, 3, 3)
movement_phase_v3(agents, GRID_SIZE, rng)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 2, unique positions: {check_unique_positions(agents)}",
)
plt.show()

This version finally preserves uniqueness, does not introduce bias and updated all states at the same time. Abviously, this strategy is more complicated, but quite usual in time and spatially discrete ABMs. One example is evacuation simulation Tanimoto et al. (2010). Here modellers even make a virtue of necessity and consider e.g. heterogeneity when it comes to add deliberate bias in selecting one of the multiple requesting agents.

In any case, it is clear that this implementation involves significantly more complex processes than those simply described as “random movement” in the conceptual model. In other words, the implementation cannot be deduced from the conceptualisation. The former must therefore be carried out with greater care.

Movement Phase - Conceptual (Revisited)¶

In the movement phase, every agent randomly moves to an empty site in its Moore neighbourhood. In case two or more agents aim to move to the same site, one of them is randomly chosen to move whereas the other ones must stay where they are.

Population Dynamics Phase - Conceptual¶

After movement, local interactions are applied according to agent type.

Predator agents:

The agent dies with a probability $p_d$ .
If there is, at least, one prey-agent in any of the cells in the agent’s neighbourhood, one of them is randomly selected and consumed with probability $p_h$ .
In case 2. happens, a new predator agent is created and located on the former place of the consumed prey agent.

Prey agents:

When at least one neighbouring cell is empty, one of them is selected at random. With probaility $p_b$ a new prey agent is created an put at this place.

Population Dynamics Phase - Implementation¶

The defined rules open several places for interpretation, in particular referring to when agents are removed from the model and when new agents are added. For example, if a predator eats a certain prey animal that, in terms of update order, comes after the predator, it might already be removed before its update. If it was the other way round, the prey could potentially spawn offspring before it is eaten. However, if one adheres to the golden rule of simultaneous updates, many of these questions answer themselves: neither agents newly created nor agents removed during this time step may influence the update of agent states.

Furthermore, one encounters a common problem when implementing population-dynamic agent-based models: updating the agent list. In particular, removing agents from a list is a computationally expensive task. For each of these agents, (a) the list must be searched for the agent and (b) the list must be reconstructed without that agent - both $\mathcal{O}(N)$ operations. There are at least three common ways to avoid this problems:

One idea is to use sets instead of lists, in which both finding and removing objects have a time complexity of $\mathcal{O}(1)$ . A drawback of this is that the iteration order is undefined, which might lead to problems with reproducibility when the iteration order is relevant.
Another idea is to simply rebuilding the agent list anew at each step. This is done, e.g. in Bicher et al. (2018).
Finally, it is also possible to simply leave all agents in the list, but only flag those who died. This is clearly problematic when memory is a problem, but is advantageous when it comes to unique ids: the index in the list can be used as agent-id.

def get_prey_neigbours(
    position: tuple[int, int],
    grid_size: tuple[int, int],
    position_agent_mapping: dict[tuple[int, int], Agent],
) -> list[Prey]:
    """
    Returns all Prey agents o neighbouring cells
    :param position: target position
    :param grid_size: size of the grid
    :param position_agent_mapping: dictionary mapping position tuples onto agents
    :returns: set with Prey agents (maximum 8 elements)
    """
    nbs = get_neighbour_positions(position, grid_size)
    prey = list()
    for nb in nbs:
        if nb in position_agent_mapping and isinstance(
            position_agent_mapping[nb], Prey
        ):
            prey.append(position_agent_mapping[nb])
    return prey


def population_dynamics_phase_v1(
    agents: list[Agent], grid_size: tuple[int, int], rng: np.random.Generator
) -> None:
    """
    Evaluates the population dynamics phase for the agents.
    I.e. predators may die, prey may produce offspring and predators eat prey

    :param agents: list of agents
    :param grid_site: size of the grid
    :param rng: random number generator
    """
    position_agent_mapping = create_position_agent_mapping(agents)
    next_agent_id = (
        max([agent.agent_id for agent in agents]) + 1
    )  # this is sloppy but sufficient for this purpose

    dying_agents = set()
    new_agents = [x for x in agents]  # agents to have in the next timestep
    for agent in agents:
        if isinstance(agent, Predator):
            if rng.random() < PREDATOR_DEATH_PROB:
                dying_agents.add(agent)
            prey_neighbors = get_prey_neigbours(
                agent.position, grid_size, position_agent_mapping
            )
            if len(prey_neighbors) > 0 and rng.random() < PREDATION_PROB:
                prey = rng.choice(prey_neighbors)
                dying_agents.add(prey)
                new_agents.append(Predator(next_agent_id, prey.position))
                next_agent_id += 1

        elif isinstance(agent, Prey):
            if rng.random() < PREY_REPRODUCTION_PROB:
                possible_positions = get_empty_neigbour_positions(
                    agent.position, grid_size, position_agent_mapping
                )
                if len(possible_positions) > 0:
                    i = rng.integers(0, len(possible_positions))
                    new_agents.append(Prey(next_agent_id, possible_positions[i]))
                    next_agent_id += 1
    agents.clear()  # clear the original list
    # (this is necessary since we operate on the same list...)
    agents.extend(
        [x for x in new_agents if x not in dying_agents]
    )  # add all new agents which have not died


plt.figure(figsize=(12, 5))
rng = np.random.default_rng(SEED)
agents = initialize_agents(GRID_SIZE, N_PREY_INITIAL, N_PREDATORS_INITIAL, rng)
plt.subplot(1, 3, 1)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 0, unique positions: {check_unique_positions(agents)}",
    show_legend=False,
    wiggle=True,
)
plt.subplot(1, 3, 2)
population_dynamics_phase_v1(agents, GRID_SIZE, rng)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 1, unique positions: {check_unique_positions(agents)}",
    show_legend=False,
    wiggle=True,
)
plt.subplot(1, 3, 3)
population_dynamics_phase_v1(agents, GRID_SIZE, rng)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 2, unique positions: {check_unique_positions(agents)}",
    wiggle=True,
)
plt.show()

It is not well visible in the plots, but the method has a rare but clear flaw in terms of uniqueness of positions. The problem lies in the creation of new prey agents. Currently, there is no guard which prevents two or more prey agents to create offsping into the same site. Note that we do not run into this problem for the creation of predator agents, because the site with the eaten prey is alway vacant.

Population Dynamics Phase - Implementation Version 2¶

From movement phase, version 3, we already know how to circumvent this problem, namely, by creating and resolving requests. In this case it is even easier to resolve, since it does not matter which prey agent eventually produced the offspring.

def population_dynamics_phase_v2(
    agents: list[Agent], grid_size: tuple[int, int], rng: np.random.Generator
) -> None:
    """
    Evaluates the population dynamics phase for the agents.
    I.e. predators may die, prey may produce offspring and predators eat prey.
    Creation of prey animals is

    :param agents: list of agents
    :param grid_site: size of the grid
    :param rng: random number generator
    """
    position_agent_mapping = create_position_agent_mapping(agents)
    next_agent_id = (
        max([agent.agent_id for agent in agents]) + 1
    )  # this is sloppy but sufficient for this purpose

    dying_agents = set()
    new_agents = [x for x in agents]  # agents to have in the next timestep
    requests_prey_birth_positions = set()
    for agent in agents:
        if isinstance(agent, Predator):
            if rng.random() < PREDATOR_DEATH_PROB:
                dying_agents.add(agent)
            prey_neighbors = get_prey_neigbours(
                agent.position, grid_size, position_agent_mapping
            )
            if len(prey_neighbors) > 0 and rng.random() < PREDATION_PROB:
                prey = rng.choice(prey_neighbors)
                dying_agents.add(prey)
                new_agents.append(Predator(next_agent_id, prey.position))
                next_agent_id += 1

        elif isinstance(agent, Prey):
            if rng.random() < PREY_REPRODUCTION_PROB:
                possible_positions = get_empty_neigbour_positions(
                    agent.position, grid_size, position_agent_mapping
                )
                if len(possible_positions) > 0:
                    i = rng.integers(0, len(possible_positions))
                    requests_prey_birth_positions.add(possible_positions[i])
    for pos in requests_prey_birth_positions:
        new_agents.append(Prey(next_agent_id, pos))
        next_agent_id += 1
    agents.clear()  # clear the original list
    # this is necessary since we operate on the same list...
    agents.extend(
        [x for x in new_agents if x not in dying_agents]
    )  # add all new agents which have not died


plt.figure(figsize=(12, 5))
rng = np.random.default_rng(SEED)
agents = initialize_agents(GRID_SIZE, N_PREY_INITIAL, N_PREDATORS_INITIAL, rng)
plt.subplot(1, 3, 1)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 0, unique positions: {check_unique_positions(agents)}",
    show_legend=False,
    wiggle=True,
)
plt.subplot(1, 3, 2)
population_dynamics_phase_v2(agents, GRID_SIZE, rng)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 1, unique positions: {check_unique_positions(agents)}",
    show_legend=False,
    wiggle=True,
)
plt.subplot(1, 3, 3)
population_dynamics_phase_v2(agents, GRID_SIZE, rng)
plot_agents(
    agents,
    GRID_SIZE,
    title=f"Step 2, unique positions: {check_unique_positions(agents)}",
    wiggle=True,
)
plt.show()

Full Model - Implementation¶

We can finally put all things together and create a full iterative model.

def run_model(steps: int, with_movement: True = True) -> list[list[Agent]]:
    rng = np.random.default_rng(SEED)
    agents = initialize_agents(GRID_SIZE, N_PREY_INITIAL, N_PREDATORS_INITIAL, rng)
    agentss = [[x.copy() for x in agents]]
    for _i in range(steps):
        if len(agents) == 0:  # agents have died out...
            agentss.append([])
        else:
            if with_movement:
                movement_phase_v3(agents, GRID_SIZE, rng)
            population_dynamics_phase_v2(agents, GRID_SIZE, rng)
            agentss.append([x.copy() for x in agents])
    return agentss


steps = 27
agentss = run_model(steps)
plt.figure(figsize=(12, 12))
for i in range(9):
    plt.subplot(3, 3, i + 1)
    k = int(i * steps / 9)
    plot_agents(agentss[k], GRID_SIZE, title=f"Step {k}", show_legend=i % 3 == 2)
plt.tight_layout()
plt.show()

Model Analysis¶

The fully implemented model is now ready for deeper analysis. One way to look at the model is given by the aggregated curves of total predator and prey animals.

def compute_totals(agentss: list[list[Agent]]) -> tuple[list[int], list[int]]:
    """
    Computes timeseries of total number of prey and predator animals
    :param agentss: agent-list for each simulation time-step
    :returns: lists or total prey and predator counts
    """
    totals = [count_agents(x) for x in agentss]
    return [x[0] for x in totals], [x[1] for x in totals]


def plot_population_curves(agentss: list[list[Agent]]) -> hv.Overlay:
    """
    Plots timeseries of total number of prey and predator animals
    :param agentss: agent-list for each simulation time-step
    :returns: lists or total prey and predator counts
    """
    prey_totals, predator_totals = compute_totals(agentss)
    time = np.arange(len(prey_totals))
    c1 = hv.Curve((time, prey_totals), label="Prey").opts(color=PREY_COLOR)
    c2 = hv.Curve((time, predator_totals), label="Predators").opts(color=PREDATOR_COLOR)

    return (c1 * c2).opts(show_legend=True, xlabel="timestep", ylabel="agents")


steps = 400
agentss_base = run_model(steps)
display(
    plot_population_curves(agentss_base).opts(
        responsive=True, height=400, max_width=800, title="Model with base parameters"
    )
)
plt.figure(figsize=(4, 4))
plot_agents(
    agentss_base[-1], GRID_SIZE, title="Final state for model with base parameters"
)
plt.show()
agentss_nomove = run_model(steps, False)
display(
    plot_population_curves(agentss_nomove).opts(
        responsive=True, height=400, max_width=800, title="Model without movement"
    )
)
plt.figure(figsize=(4, 4))
plot_agents(
    agentss_nomove[-1], GRID_SIZE, title="Final state for model without movement"
)
plt.show()

The model shows reasonable predator-prey behaviour, noticeable by the orange curves reacting upon the development of the blue curves. After an initial transient phase, both species reach a somewhat stable equilibrium. The lower plot shows the model without the movement phase. Here, the equilibrium state is much more sensitive to stochastic fluctuations due to the higher impact of locality.

Note that these values belong to this particular stochastic realisation and should not be interpreted as a (deterministic) forecast. A robust quantitative analysis would require repeated simulations across multiple random seeds and a summary of the resulting distributions, which is out of scope of this notebook.

The structural difference between the two versions (with and without movement) can also be seen on the visualisations of the final states of the simulation. Here the model without movement shows a much clearer clustering. However, the difference is even better visisble with an animated plot of the grid. Note that the corresponding visualisation routine was developed using Python’s Javascript bindings for plotting.

agent_styles = {
    "prey": {"label": "Prey", "colour": PREY_COLOR, "shape": "circle"},
    "predator": {"label": "Predators", "colour": PREDATOR_COLOR, "shape": "triangle"},
}
print("Base parameters")
animation_data = frames_to_json_data(agentss_base)

show_grid_abm_animation(
    frames=animation_data,
    grid_size=GRID_SIZE,
    agent_styles=agent_styles,
    cell_size=20,
    interval=250,
    show_population_plot=True,
)

print("Without movement")
animation_data = frames_to_json_data(agentss_nomove)

show_grid_abm_animation(
    frames=animation_data,
    grid_size=GRID_SIZE,
    agent_styles=agent_styles,
    cell_size=20,
    interval=250,
    show_population_plot=True,
)

Vaidation?¶

The model demonstrates how population-level predator–prey dynamics can emerge from simple individual rules. No equation directly determines the number of prey or predators at the next time step. Instead, the aggregate trajectory is produced by repeated movement, mortality, predation and reproduction events.

Although the model is inspired by the predator-prey system, the modelling assumptions are very extreme and far from what can be observed in reality. For example, the assumption that a prey animal that has been eaten is instantly replaced by a new predator is a highly questionable compared to what happens in reality. Nevertheless, it is interesting to note that the results do indeed correspond with the observations from the macroscopic predator-prey model as, e.g. presented in Simulation Circle. This correlation is not by chance and can e.g. be explained using so-called mean-field analysis.

Anyhow, the model should be interpreted as a conceptual ABM rather than as a validated ecological prediction and this case-study should be interpreted in terms of methods for implementation and conceptualisation of ABMs in general.

References¶

Moore, E. F. (1962). Machine Models of Self-Reproduction. In Mathematical Problems in the Biological Sciences (Vol. 14, pp. 17–34). American Mathematical Society.
Chopard, B., & Droz, M. (1998). Cellular Automata Modeling of Physical Systems. Cambridge University Press. 10.1017/CBO9780511549755
Wilensky, U. (1999). NetLogo. Center for Connected Learning and Computer-Based Modeling, Northwestern University. https://ccl.northwestern.edu/netlogo/
Tanimoto, J., Hagishima, A., & Tanaka, Y. (2010). Study of bottleneck effect at an emergency evacuation exit using cellular automata model, mean field approximation analysis, and game theory. Physica A: Statistical Mechanics and Its Applications, 389(24), 5611–5618. https://doi.org/10.1016/j.physa.2010.08.032
Bicher, M., Urach, C., & Popper, N. (2018). Gepoc ABM: a generic agent-based population model for Austria. 2018 Winter Simulation Conference (WSC), 2656–2667. 10.1109/WSC.2018.8632170