This project impelments a spatial predator-prey cellular automaton (CA) to investigate the Hydra effect, a counterintuitive ecological phenomenon where increased prey mortality paradoxically leads to higher prey population densities. The model explores the rise of emergent population dynamics through spatial structure and local interactions that differe fundemntally from well-mixed (mean-field) predictions.

The codebase uses Numba JIT compilation for computationally intensive kernels and is tuned for high performance computing environments. It supporrs parameter sweeps, finite size scaling analysis, and self-organized criticality exploration.

In classical Lotka-Volterra dynamics, increasing prey mortality always reduces the prey population. However, theoretical work has identified conditions where the opposite effect can be observed. The Hydra effect emerges in spatially structured systems where

1. Local interactions prevent predators from immediately exploiting all available prey.

2. Spatial refugia from naturally through predator-prey dynamics

3. Increased prey mortality can disrupt predator populations disproportionally.

This study uses a cellular automaton framework toi study how spatial strcuture generates the Hydra effecr and whether the system exhibits signatures of self-organized criticality at the transition point.

SOC refers to systems that naturally evolve toward a critical state without external tuning. At criticality, such systems exhibit:

- Power-law cluster size distributions: $P(s) \sim s^{-\tau}$

- Scale-free correlations: No characteristic length scale dominates

- Critical slowing down: Perturbations decay slowly near the critical point

- Fractal spatial patterns: Self-similar structure across scales

In the predator-prey context, SOC would manifest as the system self-tuning toward a critical prey mortality rate where population dynamics become scale-free and the Hydra effect is maximized.

The model uses a 2D lattice with periodic boundary conditions. Each cell occupies one of three states:

| State | Value | Description |

|-------|-------|-------------|

| Empty | 0 | Unoccupied cell |

| Prey | 1 | Prey organism |

| Predator | 2 | Predator organism |

The model uses asynchronous updates: cells are processed in random order each timestep, with state changes taking effect immediately. This prevents artificial synchronization artifacts common in parallel update schemes.

For each prey cell, in order:

- Death: With probability `prey_death`, the prey dies and the cell becomes empty

- Reproduction: If alive and a randomly selected neighbor is empty, with probability `prey_birth`, a new prey is placed in that neighbor cell

For each predator cell, in order:

- Death: With probability `predator_death`, the predator dies (starvation)

- Hunting/Reproduction: If alive and a randomly selected neighbor contains prey, with probability `predator_birth`, the predator consumes the prey and reproduces into that cell

The model supports both neighborhood types:

- Moore neighborhood: (default): 8 adjacent cells (including diagonals)

- von Neumann neighborhood: 4 adjacent cells (cardinal directions only)

Moore neighborhoods are used throughout the experiments as they provide more realistic movement and interaction patterns.

The model implements two distinct neighbor selection strategies that qualitatively affect dynamics.

In the standard mode, each organism selects a single random neighbor for interaction:

neighbor_idx = np.random.randint(0, n_neighbors)

if grid[neighbor] == EMPTY and random() < prey_birth:

grid[neighbor] = PREY

if grid[neighbor] == PREY and random() < predator_birth:

grid[neighbor] = PREDATOR

This creates a blind interaction model where organisms are not aware of their surroundings.

The directed mode implements "intelligent" neighbor selection:

- Prey: Scan all neighbors for empty cells, then randomly select one empty cell for reproduction

- Predators: Scan all neighbors for prey, then randomly select one prey cell to hunt

empty_neighbors = [n for n in neighbors if grid[n] == EMPTY]

if empty_neighbors and random() < prey_birth:

target = random.choice(empty_neighbors)

grid[target] = PREY

prey_neighbors = [n for n in neighbors if grid[n] == PREY]

if prey_neighbors and random() < predator_birth:

target = random.choice(prey_neighbors)

grid[target] = PREDATOR

This increases the effective reproduction and predation rates without changing the nominal probability parameters.

The `CA` base class provides:

- Grid initialization with specified densities

- Neighborhood definitions (Moore/von Neumann)

- Per-cell parameter evolution framework

- Validation and run loop infrastructure

The `PP` class adds:

- Predator-prey specific parameters and defaults

- Synchronous/asynchronous update dispatch

- Integration with Numba-optimized kernels

from models.CA import PP

model = PP(

rows=100,

cols=100,

densities=(0.30, 0.15), # (prey_fraction, predator_fraction)

params={

"prey_birth": 0.2,

"prey_death": 0.05,

"predator_birth": 0.8,

"predator_death": 0.05,

},

seed=42,

synchronous=False, # Always False for this study

directed_hunting=False,

)

for step in range(1000):

model.update()

prey_count = np.sum(model.grid == 1)

pred_count = np.sum(model.grid == 2)

The model supports per-cell parameter evolution, where offspring inherit (with mutation) their parent's parameter values:

model.evolve(

"prey_death",

sd=0.01, # Mutation standard deviation

min_val=0.0, # Minimum allowed value

max_val=0.20, # Maximum allowed value

)

When evolution is enabled, each prey cell maintains its own `prey_death` value. Upon reproduction, offspring inherit the parent's value plus Gaussian noise. This allows investigation of whether the population self-organizes toward a critical mortality rate.

The computational bottleneck is the update kernel, which must process every occupied cell each timestep. For a 1000×1000 grid with 50% occupancy, this means ~500,000 cell updates per step.

Key optimizations:

- **JIT Compilation**: Core kernels use `@njit(cache=True)` for ahead-of-time compilation

2. **Pre-allocated Buffers**: The `PPKernel` class maintains reusable arrays to avoid allocation overhead

3. **Efficient Shuffling**: Fisher-Yates shuffle implemented in Numba for random cell ordering

4. **Cell Lists for PCF**: Pair correlation functions use spatial hashing for O(N) instead of O(N²) complexity

class PPKernel:

"""Wrapper for predator-prey kernel with pre-allocated buffers."""

def __init__(self, rows, cols, neighborhood="moore", directed_hunting=False):

self.rows = rows

self.cols = cols

self.directed_hunting = directed_hunting

self._occupied_buffer = np.empty((rows * cols, 2), dtype=np.int32)

# Neighbor offset arrays

if neighborhood == "moore":

self._dr = np.array([-1, -1, -1, 0, 0, 1, 1, 1], dtype=np.int32)

self._dc = np.array([-1, 0, 1, -1, 1, -1, 0, 1], dtype=np.int32)

else:

self._dr = np.array([-1, 1, 0, 0], dtype=np.int32)

self._dc = np.array([0, 0, -1, 1], dtype=np.int32)

Clusters are contiguous groups of same-species cells (using Moore connectivity). The implementation uses stack-based flood fill with periodic boundary conditions.

from models.numba_optimized import get_cluster_stats_fast

stats = get_cluster_stats_fast(grid, species=1) # Prey clusters

print(f"Number of clusters: {stats['n_clusters']}")

print(f"Largest cluster: {stats['largest']} cells")

print(f"Largest fraction: {stats['largest_fraction']:.3f}")

Metrics collected:

- Cluster size distribution: Power-law indicates criticality

- Largest cluster fraction: Order parameter for percolation transition

The `Config` dataclass centralizes all experimental parameters:

from config import PHASE1_CONFIG, Config

cfg = PHASE1_CONFIG

cfg = Config(

grid_size=500,

n_replicates=20,

warmup_steps=1000,

measurement_steps=2000,

collect_pcf=True,

pcf_sample_rate=0.5,

)

print(f"Grid: {cfg.grid_size}x{cfg.grid_size}")

print(f"Estimate: {cfg.estimate_runtime(n_cores=32)}")

Each phase has a predefined configuration (`PHASE1_CONFIG` through `PHASE5_CONFIG`) with appropriate defaults for that analysis.

Results are saved in JSONL format (one JSON object per line) for efficient streaming and parallel processing:

{"prey_birth": 0.2, "prey_death": 0.05, "seed": 12345, "final_prey": 28500, ...}

{"prey_birth": 0.2, "prey_death": 0.06, "seed": 12346, "final_prey": 27200, ...}

Each result dictionary contains:

- Input parameters

- Final population counts

- Cluster statistics

- Evolution statistics (if enabled)

- Time series (if enabled)

Metadata files (`phase{N}_metadata.json`) accompany each results file with configuration details and runtime information.

**Required:**

- Python 3.8+

- NumPy

- Numba (for JIT compilation)

- tqdm (progress bars)

- joblib (parallelization)

**Optional:**

- matplotlib (visualization)

- scipy (additional analysis)