Dataset DOI: 10.5061/dryad.9w0vt4bsk

Related publication

This dataset supports the following article:

Tomizuka, H., and Y. Tachiki. 2026. Maintenance and breakdown: how population dynamics
determine the evolutionary fate of Batesian mimicry. The American Naturalist.
https://doi.org/10.1086/743066

Description of the data and file structure

This repository contains all code and simulation output required to reproduce Figures 2–5 of the associated article.

The article develops a mathematical model of a Batesian mimicry complex consisting of an unpalatable model-species (which predators learn to avoid) and a palatable mimic-species (which gains protection by resembling the model-species). Predators decide whether to attack a prey item using a signal-detection rule. The mimic-species evolves a one-dimensional continuous appearance trait z; the model-species has a fixed trait z_D = 5.0. The dataset comprises (i) deterministic/analytical results computed in Mathematica and (ii) stochastic agent-based simulations written in C, whose outputs are stored as CSV files.

Two ecological regimes recur throughout the dataset and are referred to by name:

• Type I — the model-extinction equilibrium is locally stable (i.e. r_D - S_D - alpha_D * P < 0). The model-species cannot increase from low density when predators have not learned to avoid it.
• Type II — the model-extinction equilibrium is unstable (i.e. r_D - S_D - alpha_D * P > 0). The model-species can increase from low density.

Directory layout

Fig2A/                          Mathematica only; no tabular data
    fig2A.nb, fig2a.m, params_set1.m ... params_set6.m, util.m
Fig2BC/                         Mathematica only; no tabular data
    fig2BC.nb
Fig3-1/                         Mathematica only; no tabular data
    fig3-1.nb, params_set1.m, params_set2.m, util.m
Fig3-2/
    Fig3-2.nb
    C_code/
        Agent_Based_Simulation.c
        Dynamics.csv, Evolution.csv
        include/                (place MT.h here; see "Code/software")
Fig4/
    fig4.nb
    Type1_Simulation/
        Agent_Based_Simulation.c
        Dynamics.csv, Evolution.csv
        include/
    Type2_Simulation/
        Agent_Based_Simulation.c
        Dynamics.csv, Evolution.csv
        include/
    Fitness_Calculation/
        Fitness_Calculation.c
        Fitness.csv
        include/
Fig5/
    Fig5.nb
    C_code/
        Type1_Simulation/
            Agent_Based_Simulation_EnvChange.c
            Dynamics.csv, Evolution.csv
            include/
        Type2_Simulation/
            Agent_Based_Simulation_EnvChange.c
            Dynamics.csv, Evolution.csv
            include/

Fig2A/, Fig2BC/ and Fig3-1/ contain no tabular data files: those figures are produced entirely by numerical analysis inside the Mathematica notebooks, using the parameter values supplied in the accompanying .m files. All tabular data in this dataset are the CSV files listed above.

Figures reproduced by each folder

Units and conventions used in all figures. Time is expressed in dimensionless model time units (all rates are per unit time). Population sizes are counts of individuals. The phenotype z is a dimensionless, one-dimensional index of overall appearance (colour and pattern): z = z_D = 5.0 corresponds to perfect resemblance to the model-species, and z = 0 corresponds to the non-mimetic phenotype, which carries the lowest physiological cost.

Abbreviations used in the figures.

Fig2A/ — Figure 2A

Bifurcation diagrams of the ecological subsystem. For each fixed mimic phenotype z_M (x-axis, dimensionless, 0–5), the panels show the equilibrium ratio of model-species to mimic-species population size, D*/M* (y-axis, dimensionless). Solid curves are stable equilibria; dashed curves are unstable equilibria.

In the Type I rows a stable model-extinction equilibrium (D*/M* = 0) coexists with a stable coexistence equilibrium, so the ecological outcome depends on the initial densities. In the Type II rows the model-extinction equilibrium is unstable and coexistence is the only stable outcome.

Fig2BC/ — Figure 2B and 2C

Figure 2B. Phase planes for three of the equilibria labelled in Figure 2A. The x-axis is the density of the model-species and the y-axis the density of the mimic-species (both in individuals). The red curve is the model-species nullcline (dD/dt = 0) and the blue curve the mimic-species nullcline (dM/dt = 0); grey arrows show the direction of the ecological dynamics.

Figure 2C. Parameter regions producing each regime. The x-axis is the predation coefficient on the model-species (alpha_D) and the y-axis is the intrinsic growth rate of the model-species (r_D). Columns vary the background mortality rate of the model-species (S_D = 0.0, 0.5, 1.0) and rows vary the predator density (P = 0.1, 1.0, 5.0). Grey = Type I; white = Type II; blue = r_D - S_D < 0, i.e. the model-species is inviable irrespective of mimicry. The straight boundary is r_D - S_D - alpha_D * P = 0.

Fig3-1/ — Figure 3, bifurcation diagrams and PIPs

Adaptive-dynamics analysis of how the model-species carrying capacity shapes the evolution of the mimic phenotype, for Type I (alpha_D = 0.6) and Type II (alpha_D = 0.25).

Bifurcation diagrams. Evolutionarily singular strategies z_M* (y-axis, dimensionless) as a function of the model-species carrying capacity K_D (x-axis, individuals).
Grey solid = CSS; pink solid = evolutionary branching point (disruptive selection, which can generate polymorphism); grey dashed = evolutionary repeller.

PIPs. For selected values of K_D, the x-axis is the resident (wild-type) mimic phenotype z_M and the y-axis the mutant phenotype z'. Dark regions indicate W(z'; z_M) > 0, i.e. a rare mutant can invade; white regions indicate W(z'; z_M) < 0. Blue regions mark resident phenotypes for which no stable coexistence equilibrium exists, so the ecological system settles at the model-extinction equilibrium; there, invasion fitness is evaluated at that equilibrium.

Fig3-2/ — Figure 3, individual-based simulation panel

A single stochastic agent-based simulation in the Type I regime with K_D = 5000. The x-axis is simulation time and the y-axis the mimic phenotype z.

Starting from a non-mimetic population (z = 0), mimicry initially increases because more mimetic mutants are attacked less often. The improved mimicry, however, increases predation on the model-species and drives it extinct (around t ≈ 250,000). Once the model-species is gone, mimicry confers no protection while still incurring a physiological cost, so selection reverses and the phenotype returns to z ≈ 0. The panel therefore shows a transient gain of mimicry followed by its loss under constant external conditions.

Data files: Dynamics.csv (population sizes) and Evolution.csv (phenotypes).

Fig4/ — Figure 4

Whether evolutionary branching produces a lasting polymorphism differs between the two regimes.

• Panel A (Type1_Simulation/) — Type I (alpha_D = 0.6, K_D = 12,000, K_M = 50,000). Mimic phenotype trajectories from an agent-based simulation. The population briefly splits into a more mimetic and a less mimetic lineage, but the polymorphism collapses. Plotted from Type1_Simulation/Evolution.csv.
• Panel C (Type1_Simulation/) — population dynamics of the same run. The y-axis is population size (individuals); the red line is the model-species and the blue line the mimic-species. The model-species goes extinct, which removes the benefit of mimicry and causes the collapse seen in panel A. Plotted from Type1_Simulation/Dynamics.csv.
• Panel B (Type2_Simulation/) — Type II (alpha_D = 0.25, K_D = 25,000, K_M = 50,000). Here branching produces a persistent dimorphism: one lineage converges on accurate mimicry (z ≈ z_D = 5.0) and the other on the non-mimetic phenotype (z ≈ 0). Plotted from Type2_Simulation/Evolution.csv.
• Panel D (Fitness_Calculation/) — the mechanism maintaining that dimorphism in Type II. The y-axis is the per-capita growth rate of each morph; the x-axis is the ratio of mimetic-morph abundance to model-species abundance (dimensionless). Each point is one recorded time step of a simulation in which the mimic population is held as exactly two fixed morphs (mutation switched off): a mimetic morph with z = z_D = 5.0 and a non-mimetic morph with z = 0. When the mimetic morph is rare relative to the model-species, model-like prey are mostly genuinely unpalatable, deception works, and the mimetic morph has the higher fitness. When the mimetic morph becomes common, model-like prey increasingly include palatable mimics, predator avoidance weakens, and the cheaper non-mimetic morph is favoured. This reversal is the negative frequency-dependent selection that maintains the polymorphism. Plotted from Fitness_Calculation/Fitness.csv.

Fig5/ — Figure 5

The response of mimicry evolution to a temporary environmental perturbation. The x-axis is simulation time and the y-axis the mimic phenotype z.

In both panels the carrying capacity of the model-species is changed on the following schedule:

The carrying capacity of the mimic-species is held at K_M = 10,000 throughout.

• Panel A (Type1_Simulation/) — Type I (alpha_D = 0.6). The temporary reduction in K_D pushes the system to the model-extinction equilibrium. The mimic phenotype evolves to the non-mimetic state (z ≈ 0) and does not recover when the environment is restored: mimicry is lost irreversibly.
• Panel B (Type2_Simulation/) — Type II (alpha_D = 0.25). The mimetic/non-mimetic dimorphism persists throughout the perturbation and the system returns to its pre-perturbation state once conditions are restored.

Both panels are plotted from the corresponding Evolution.csv.

Files and variables

Notation: variable names in the C code vs. symbols in the article

Because this README is plain Markdown, subscripts are written with an underscore and Greek letters are spelled out. For example, K_D here is the symbol printed as K with subscript D in the article, and alpha_D, sigma_N and delta_max correspond to the Greek symbols in the article. An asterisk marks an equilibrium or singular value (for example D*, M*, z*, z_M*).

The mimic-species mortality rate is S(z) = (delta_max - delta_min) * abs(z)^{n / z_D}n + delta_min, i.e. the physiological cost of mimicry rises with resemblance to the model-species.

File: Dynamics.csv

Time series of population sizes. Present in Fig3-2/C_code/, Fig4/Type1_Simulation/, Fig4/Type2_Simulation/, Fig5/C_code/Type1_Simulation/ and Fig5/C_code/Type2_Simulation/. One row per recorded time point. 4 columns.

File: Evolution.csv

The phenotype distribution of the mimic-species over time. Present in the same five folders as Dynamics.csv. 3 columns.

Important. At each recorded time point, one row is written for every trait bin that contains at least one individual; empty bins are omitted. Consequently many consecutive rows share the same value of time, and the number of rows per time point varies with the phenotypic spread of the population. Summing column 3 over all rows with the same time reproduces M in Dynamics.csv for that time point.

File: Fitness.csv

Present only in Fig4/Fitness_Calculation/. In this simulation the mutation rate is set to zero and the mimic population consists of exactly two fixed morphs: a non-mimetic morph (z = 0) and a mimetic morph (z = z_D = 5.0). The file records their abundances and realised fitnesses as the populations oscillate. 7 columns.

Fitness (columns 6 and 7) is computed for a morph with phenotype z as

w(z) = r_M(1 - M/K_M) - S(z) - alpha_M R_M(z) P

i.e. the per-capita birth rate minus the physiological (background) mortality rate minus the predation rate. A positive value means the morph is increasing at that instant.

Code/software

Computational environment

• C compiler: GCC (TDM-GCC 10.3.0) on Windows 11 (64-bit)
• Mathematica: Wolfram Mathematica 13.1.0.0
• Workflow: the stochastic simulations are run in C and write CSV files; the CSV files are then read, analysed and plotted in the Mathematica notebooks. The Mathematica notebooks for Figures 2 and 3-1 perform their numerical analysis internally and require no CSV input.

External dependency: Mersenne Twister

The C programs draw random numbers with the Mersenne Twister (MT19937) generator of Matsumoto and Nishimura, which is not redistributed here. Each include/ directory is provided empty; before compiling, place a single-header implementation of the generator in it, named MT.h. The programs use three functions:

• void init_genrand(unsigned long s) — seed the generator
• double genrand_real1(void) — uniform random number on [0, 1]
• double genrand_real2(void) — uniform random number on [0, 1)

Reproducibility of the stochastic simulations

The random number generator is seeded from the system clock, so re-running a simulation produces a different realisation and the CSV files are not reproduced byte-for-byte; the qualitative outcomes are robust across realisations. The CSV files deposited here are the exact realisations used to produce the figures in the article.

Access information

All data and code in this repository were generated by the authors. No external or third-party datasets were used.

Contact

Haruto Tomizuka
Department of Biological Sciences, Tokyo Metropolitan University, Tokyo 192-0397, Japan
Email: haruto.tomizuka.ecol@gmail.com
ORCID: 0009-0000-2767-226X