https://doi.org/10.5061/dryad.c2fqz61k7

Description of the data and file structure

This dataset contains the C-language simulation code and Mathematica scripts used to generate and visualize the results presented in the forthcoming paper:

Tomizuka, H. & Tachiki, Y. (2025). Eco-evolutionary metapopulation dynamics of Batesian mimicry: Conditions for mimics without models. Journal of Theoretical Biology. (https://doi.org/10.1016/j.jtbi.2025.112299)

Short Summary of the Study
Batesian mimicry is adopted by palatable prey species (mimics) to avoid predator attack by resembling unpalatable species (models). Despite numerous studies on this phenomenon, several aspects of its evolution remain unclear. One of the most interesting questions is whether mimics can inhabit allopatrically from their models. Classical theory suggests that mimics should only inhabit sympatrically with their models because predators need to learn the model phenotype. However, several studies have reported mimics that live outside the distribution range of their models (mimics without models; MWM). In this study, we constructed and analyzed an eco-evolutionary dynamics model of Batesian mimicry incorporating two populations under different environmental conditions to clarify the conditions necessary for MWM to occur. We identified that MWM occurred when specific conditions are satisfied by the following six parameters: (1) carrying capacities; (2) migration rate of mimics; (3) migration rate of predators; (4) evolution rate of mimic phenotype; (5) toxicity of the model and nutritional level of the mimic; and (6) population density of predators. Our findings differed from previous predictions that considered population and evolutionary dynamics independently. Overall, this study asserts that the interaction between ecological and evolutionary processes makes the conditions for MWM more stringent.

Overview of Folders/Files
The repository is organized into seven folders, each corresponding to a figure or supplementary figure in the study. The files within each folder are as follows:

File–Figure Crosswalk

Note: Panel labels (A–H) match the article’s figures.

---

Brief descriptions of each figure panel

Conditions under which mimics without models occur (Fig.3)

Axes for all panels: vertical = carrying capacity of model-species at site 1 (KD,1); horizontal = carrying capacity at site 2 (KD,2).

• (A) Final density of model-species at site 1 under local learning (q = 1) with mimic migration (mM = 0.027).
  Data: Fig3a.csv (from Fig3a.c)
• (B) Final mean phenotype of mimic-species at site 1 (zM,1) under the same condition as (A); darker region denotes evolution toward crypsis; lighter/white toward mimetic trait.
  (source file same as (A))
• (C) Final model density at site 1 under local learning, no mimic migration (mM = 0).
  Data: Fig3b.csv (from Fig3b.c)
• (D) Final zM,1 under the same condition as (C).
  (source file same as (C))
• (E) Final model density at site 1 under global learning (q = 0.5) with mimic migration (mM = 0.027).
  Data: Fig3c.csv (from Fig3c.c)
• (F) Final zM,1 under the same condition as (E).
  (source file same as (E))
• (G) Final model density at site 1 under global learning, no mimic migration.
  Data: Fig3d.csv (from Fig3d.c)
• (H) Final zM,1 under the same condition as (G).
  (source file same as (G))
  Visualization: Fig3.nb

The effect of the migration rate on the outcome (Fig.4)

• (A) Phase diagram of final model density at site 1 across model migration (mD) and mimic migration (mM); orange indicates sustained oscillations (limit cycles).
  Data: Fig4ab.csv (from Fig4ab.c)
• (B) Phase diagram of final zM,1 across (mD, mM).
  Data: Fig4ab.csv (from Fig4ab.c)
• (C) Time series at mD = 0.0005, mM = 0.0015: model density at site 1 stays near KD,1 while the mimic at site 1 evolves toward crypsis.
  Data: Fig4c.csv (from Fig4c.c)
• (D) Time series at mD = 0.0005, mM = 0.01: densities and z1 show sustained oscillations.
  Data: Fig4d.csv (from Fig4d.c)
• (E) Time series at mD = 0.0005, mM = 0.027 : model goes locally extinct at site 1 while site 2 retains coexistence; gene flow maintains mimetic trait at site 1.
  Data: Fig4e.csv (from Fig4e.c)
• (F) Time series at mD = 0.004, mM = 0.0015: initial oscillations decay to equilibrium with coexistence.
  Data: Fig4f.csv (from Fig4f.c)
  Visualization: Fig4.nb

The effect of predator memory on the outcome (Fig.5)

Bifurcation diagrams versus q; solid = site 1, dotted = site 2; red/blue = model/mimic densities; green = zM; gray = attack probability on model.

• (A) mM = 0:
  Data: Fig5a.csv (from Fig5a.c)
• (B) mM = 0.01:
  Data: Fig5b.csv (from Fig5b.c)
• (C) mM = 0.03:
  Data: Fig5c.csv (from Fig5c.c)
  Visualization: Fig5.nb

The effect of other parameters on the outcome (Fig.6)

Bifurcation diagrams (solid = site 1, dotted = site 2; colors as above).

• (A) Evolution rate G:
  Data: Fig6a.csv (from Fig6a.c)
• (B) Cost-to-benefit ratio c/b:
  Data: Fig6b.csv (from Fig6b.c)
• (C) Predator density P:
  Data: Fig6c.csv (from Fig6c.c)
  Visualization: Fig6.nb

Simulation of predator learning process (Fig.S2)

Simulates how predators memory retention time and migration rate affect to predation behavior.
Data: FigS2.csv (from FigS2.c); visualization: FigS2.nb

Parameter sensitivity of the predation threshold (Fig.S3)

Sensitivity of the attack threshold x* to key quantities (e.g., relative abundances, phenotype means, c/b).
Data: FigS3.csv (from FigS3.c); visualization: FigS3.nb

Outcomes considering only population dynamics (Fig.S4)

Ecological consequences of mimicry complex. (A & B) The vertical and horizontal axes represent the migration rates of model- and mimic-species, respectively. (A) Population density of model-species at site 1. (B) Phenotype of mimic-species at site 1.
Data: FigS4.csv (from FigS4.c); visualization: FigS4.nb

The effect of the migration rate of mimics on the outcome (Fig.S5)

Dependence on migration rate of mimic-species (mM). Solid lines correspond to site 1, and dotted lines correspond to site 2. (A) Equilibrium population density of model- and mimic-species. (B) Equilibrium mean phenotype of mimic-species. (C) Predation threshold z_star decided by predators.
Data: FigS5.csv (from FigS5.c); visualization: FigS5.nb

The outcome when the model went extinct (Fig.S6)

Bifurcation diagrams of our system to vary the local learning rate (q). Solid lines correspond to site 1, and dotted lines correspond to site 2. (A) Equilibrium population density of model- and mimic-species. (B) Equilibrium the phenotype of mimic-species.
Data: FigS6.csv (from FigS6.c); visualization: FigS6.nb

Each folder contains:

・[.c]: C-language source code for evolutionary simulations. The main script for each experiment is provided here.

・[.csv]: Rdata generated from simulations.

・[.nb,  .pdf]: Mathematica (.nb) scripts for generating figures. A corresponding PDF file is also included for reference.

Column order (for CSV files without headers)

Global conventions

t = time;
D1, D2 = population densities of the model at sites 1 and 2;
M1, M2 = population densities of the mimic at sites 1 and 2;
z1, z2 = mean phenotype of the mimic at sites 1 and 2;
KD1, KD2 / KM1, KM2 = carrying capacities of model / mimic at sites 1 and 2;
mD, mM = per-capita migration rates between sites (model / mimic);
q = weight on local information in predator decision-making;
G = evolution rate of the mimic phenotype;
P = predator density;
c_over_b = cost-to-benefit ratio;
z1_star, z2_star = predation thresholds at sites 1 and 2 (appear only in some files);
#_max, #_min = the maximum/minimum of variable # within a simulation;
Learning_time = predator memory retention time;
migration_rate = predator migration rate between sites.

Per-file column order

• Fig3a.csv : KD1, KD2, D1, M1, z1, D2, M2, z2

• Fig3b.csv : KD1, KD2, D1, M1, z1, D2, M2, z2

• Fig3c.csv : KD1, KD2, D1, M1, z1, D2, M2, z2

• Fig3d.csv : KD1, KD2, D1, M1, z1, D2, M2, z2

• Fig4ab.csv : mD, mM, D1, M1, z1, D2, M2, z2

• Fig4c.csv : t, D1, M1, z1, D2, M2, z2, z1_star, z2_star

• Fig4d.csv : t, D1, M1, z1, D2, M2, z2, z1_star, z2_star

• Fig4e.csv : t, D1, M1, z1, D2, M2, z2, z1_star, z2_star

• Fig4f.csv : t, D1, M1, z1, D2, M2, z2, z1_star, z2_star

• Fig5a.csv : q, D1_max, M1_max, z1_max, D2_max, M2_max, z2_max, D1_min, M1_min, z1_min, D2_min, M2_min, z2_min, RD1_max, RD1_min, RD2_max, RD2_min

• Fig5b.csv : q, D1_max, M1_max, z1_max, D2_max, M2_max, z2_max, D1_min, M1_min, z1_min, D2_min, M2_min, z2_min, RD1_max, RD1_min, RD2_max, RD2_min

• Fig5c.csv : q, D1_max, M1_max, z1_max, D2_max, M2_max, z2_max, D1_min, M1_min, z1_min, D2_min, M2_min, z2_min, RD1_max, RD1_min, RD2_max, RD2_min

• Fig6a.csv : G, D1_max, M1_max, z1_max, D2_max, M2_max, z2_max, D1_min, M1_min, z1_min, D2_min, M2_min, z2_min, RD1_max, RD1_min, RD2_max, RD2_min

• Fig6b.csv : c_over_b, D1_max, M1_max, z1_max, D2_max, M2_max, z2_max, D1_min, M1_min, z1_min, D2_min, M2_min, z2_min, RD1_max, RD1_min, RD2_max, RD2_min

• Fig6c.csv : P, D1_max, M1_max, z1_max, D2_max, M2_max, z2_max, D1_min, M1_min, z1_min, D2_min, M2_min, z2_min, RD1_max, RD1_min, RD2_max, RD2_min

• FigS2.csv : Learning_time, migration_rate, q

• FigS3.csv : z1_star, z2_star, q, D1, M1, z1, D2, M2, z2, c_over_b

• FigS4.csv : mD, mM, D1, M1, z1, D2, M2, z2

• FigS5.csv : mM, D1_max, M1_max, z1_max, D2_max, M2_max, z2_max, D1_min, M1_min, z1_min, D2_min, M2_min, z2_min, RD1_max, RD1_min, RD2_max, RD2_min

• FigS6.csv : q, D1_max, M1_max, z1_max, D2_max, M2_max, z2_max, D1_min, M1_min, z1_min, D2_min, M2_min, z2_min, RD1_max, RD1_min, RD2_max, RD2_min

Description of Variables

The meanings of each column name are shown below.

Sharing/Access Information

Data was derived from the following sources: The dataset is generated entirely through computational simulations and does not contain empirical data.

Code/Software

Computational Environment

Operating System: Windows 11 Pro

Processor: 12th Gen Intel Core i9-12900H (2.90 GHz)

RAM: 64 GB

Software Used

C Compiler: GCC (TDM-GCC version 10.3.0)

Data Visualization: Wolfram Mathematica (version 13.1.0.0)