This dataset contains R scripts and simulation outputs associated with the article:
Adaptive harvesting of two trophic levels stabilises predator–prey dynamics: Simulations with Eurasian lynx and European roe deer.

The study implements two stochastic stage-structured matrix population models integrating European roe deer (Capreolus capreolus) and Eurasian lynx (Lynx lynx) linked in a coupled predator–prey system. The model evaluates how adaptive harvesting strategies applied to one or both trophic levels influence predator–prey dynamics over a 25-year simulation period.

Simulations were replicated 500 times for each harvest management scenario and combination of initial population densities. The provided .rds files correspond to the output of the simulations used to generate the results and figures presented in the associated publication.

Because demographic stochasticity is included in vital rates, rerunning the simulations may produce slightly different realisations of population trajectories, but will not change the general trend due to the 500 replicates. However, the archived output files allow full reproduction of all published figures and results.

Description of the data and file structure

The dataset contains R scripts defining model parameters and functions, stochastic processes, simulation, and plotting functions, as well as the simulation outputs of the article.

File structure

MPM_load_parameters.R
Defines baseline demographic parameters for roe deer and lynx, including:

• Initial population densities and stage structure
• Reproduction and recruitment parameters
• Stage-specific survival and mortality rates
• Parameters related to predator functional response

MPM_stochastic_parameters.R
Defines stochastic variation applied to species’ vital rates (reproduction and mortality).

MPM_packages_and_functions.R
Loads required R packages and defines all model functions, including:

• Matrix population model construction
• Reproduction and survival processes
• Lynx predation function
• Hunting mortality function
• Red fox predation function
• Density dependence functions
• Functional response function on kill rate
• Data preparation functions (population density summaries, quasi-extinction proportion)
• Plotting functions

MPM_simulations.R
Main simulation script which defines simulation parameters before running such as:

• Number of replicates (500 per scenario)
• Simulation length (25 years)
• Harvest management scenarios
• Combinations of initial population densities

The script runs the stochastic simulations and saves output files as .rds objects in a folder named Results.

MPM_plots.R
Generates figures presented in the manuscript and supplementary material using the saved .rds output files.

Description of simulation output files

Two output files are provided in the Results.zip folder:

• AllSimLyn.rds – Lynx population simulations
• AllSimRoe.rds – Roe deer population simulations

Each .rds file contains nested lists structured as:

Scenario
 └── Initial density
      └── Simulation replicates
           └── Population trajectories

Where:

• "Scenario" represents a unique harvest management strategy.
• "Initial density" represents a specific combination of starting densities for roe deer and lynx populations.
• "Simulation replicates" represent independent stochastic simulations.
• "Population trajectories" represent annual abundance over the 25-year projection period.

In addition, two output files are provided in the Sensitivity.zip folder:

• SimLynRatio.rds – Lynx population simulations
• SimRoeRatio.rds – Roe deer population simulations

Each .rds file contains nested lists structured as:

Scenario
 └── Ratio variation
      └── Simulation replicates
           └── Population trajectories

Where:

• "Scenario" represents a unique harvest management strategy (only two relevant scenario for the sensitivity are simulated here: scenario 7 and 16).
• "Ratio variation" represents a variation in the density threshold for the numerical response (ratio of roe deer/lynx of: 0.3/0.003 or 1.5/0.003)
• "Simulation replicates" represent independent stochastic simulations.
• "Population trajectories" represent annual abundance over the 25-year projection period.

Population abundance is expressed as total number of individuals per species per year. No missing data codes are used; all trajectories are fully simulated over the projection period for each of the species stages.

The .rds files contain the complete simulation results used to generate manuscript figures and supplementary materials.

Sharing/Access information

This dataset is publicly available via the Dryad Digital Repository.

Links to other publicly accessible locations of the data:

• GitHub repository: https://github.com/CecileAEC/MPM_roedeer-lynx

Data was not derived from external archived datasets. All simulation outputs were generated using the model code provided in this archive.

Code and software version

All simulations were conducted using:

R version 4.2.0 (2022-04-22 ucrt) – “Vigorous Calisthenics”
Platform: x86_64-w64-mingw32/x64 (64-bit)

All required R packages are loaded within MPM_packages_and_functions.R.

Workflow

To reproduce analyses:

1. Create the Results folder to store the output of the simulation.
2. Run MPM_simulations.R to generate simulation output files (AllSimLyn.rds and AllSimRoe.rds).
3. Run relevant sections of MPM_plots.R to reproduce figures presented in the manuscript and supplementary material.

Simulation parameters (number of replicates, simulation time span, harvest scenarios, and initial densities) can be modified at the beginning of MPM_simulations.R for exploratory or reduced test runs.

The .rds files in .zip folders correspond exactly to the simulations used in the published article.
OPTIONAL: To produce the exact figures of the manuscript, unzip the Results.zip and Sensitivity.zip folders to access the .rds of the manuscript simulations outputs.