Symbolic regression for empirically realistic population dynamic time series

This study evaluates the effectiveness of symbolic regression in recovering population dynamic models from realistic field-based time series, revealing that while high sampling density enables equation recovery, the method's success is limited by process noise and current evaluation workflows often fail to consistently identify the true model among candidates.

The Gist

Imagine you are a detective trying to solve a mystery: How does a specific population of giant kelp grow and shrink over time?

In the past, scientists would guess the rules of the game based on their intuition (like guessing a recipe by tasting the soup). But today, we have a powerful new tool called Symbolic Regression. Think of this tool as a super-smart, robotic chef that looks at a pile of data (the soup) and tries to reverse-engineer the exact recipe (the mathematical equation) that created it.

This paper asks a very practical question: Does this robotic chef work when the data is messy, like real life, or does it only work in the perfect, sterile kitchen of a computer simulation?

Here is the breakdown of their investigation, explained simply:

1. The Setup: The "Kelp Factory"

The researchers didn't just look at real kelp; they built a digital "Kelp Factory." They created a perfect, known recipe for how kelp grows (a complex equation involving time delays, like how long it takes a baby kelp to grow up).

• The Goal: Feed the data from this factory into the robotic chef and see if it can figure out the original recipe.
• The Twist: They didn't just feed it perfect data. They messed it up to mimic real-world problems:
  • Low Sampling Density: Instead of taking a photo of the kelp every second, they took a photo only once every few days (or even once a week).
  • Process Noise: They added "chaos" to the system, like random storms or temperature spikes that make the kelp grow unpredictably.
  • Fake Clues: They added extra variables that had nothing to do with the kelp (like the number of seagulls) to see if the robot would get distracted.

2. The Investigation: The "Four Detectives"

Once the robotic chef generated a list of possible recipes (equations), the researchers had to pick the right one. They tested four different ways to choose the winner:

1. The Visual Detective: Looking at a graph and picking the simplest recipe that fits well.
2. The Logarithmic Detective: A slightly different way of looking at the graph.
3. The Scorekeeper: A computer algorithm that automatically picks the best balance between simplicity and accuracy.
4. The Statistician: Using a strict mathematical rule (BIC) to penalize complex recipes.

3. The Findings: What Worked and What Didn't

The "Too Few Photos" Problem (Sampling Density)
This was the biggest deal-breaker.

• The Analogy: Imagine trying to guess the plot of a movie by watching only 5 random frames. You might guess the genre, but you won't know the story.
• The Result: If the researchers took fewer than 10 to 25 photos per cycle of the kelp's growth, the robotic chef failed completely. It couldn't find the recipe.
• The Good News: Once they took 50 or more photos per cycle, the chef started getting it right. It could find the true recipe, even with the "chaos" (noise) added in.

The "Chaos" Surprise (Process Noise)

• The Analogy: Usually, we think noise is bad. But here, the "random storms" actually helped!
• The Result: Surprisingly, adding a little bit of chaos made the data easier to understand. It forced the kelp to explore different growth states, giving the robotic chef more clues to work with. It's like shaking a box of puzzle pieces to help them fall into place.

The "Fake Clue" Trap (Spurious Variables)

• The Result: When the data was high-quality (lots of photos), the robot ignored the fake clues (seagulls) and focused on the real ones. But when the data was sparse, the robot got confused and started blaming the seagulls for the kelp's growth.

The "Selection" Problem
This is the most critical finding.

• The Result: Even when the robotic chef did find the perfect, true recipe, the "Detectives" (the selection workflows) often missed it. They picked a slightly different, simpler-looking recipe instead.
• The Analogy: It's like the chef cooks the perfect dish, but the judge picks a slightly different dish because it looks prettier on the plate, even though it tastes worse. The true answer was there, but the tools to pick the winner weren't good enough.

4. The Bottom Line

Symbolic regression is a powerful tool, but it has strict requirements:

1. You need a lot of data: You can't just check your population once a year. You need to check it frequently (at least 25–50 times per growth cycle) to get a clear picture.
2. A little chaos is okay: Random environmental changes might actually help you understand the system better than a perfectly calm one.
3. We need better judges: The algorithm is great at cooking (finding the equation), but we need better ways to taste-test (select the equation). Currently, the tools we use to pick the "best" equation often miss the true one.

In short: If you want to use this technology to understand nature, make sure you have high-quality, frequent data, and be very careful about how you choose the final answer. The robot can do the math, but humans still need to be smart about how they interpret the results.

Technical Summary

1. Problem Statement

While Symbolic Regression (SR) has shown promise in reverse-engineering mechanistic population models (e.g., Logistic growth, Lotka-Volterra) from idealized, densely sampled simulated or laboratory data, its applicability to empirical field data remains unproven. Field-based time series often suffer from:

• Low sampling density: Far fewer data points per population cycle compared to previous benchmarks.
• Process noise: Inherent stochasticity in population dynamics, distinct from observation noise.
• Asymmetry: Real-world cycles often exhibit rapid growth and slow decay, unlike the symmetric cycles used in prior studies.
• Spurious variables: In field contexts, researchers often do not know the true causal drivers a priori and must test extraneous variables.
• Selection ambiguity: Even if the correct equation is generated, there is a lack of robust criteria to select it from a suite of competing candidates along the Pareto frontier (the trade-off between complexity and fit).

The authors aim to quantify SR's utility under these realistic constraints and evaluate different workflows for model selection.

2. Methodology

Generative Model

The study utilized a modified Bence & Nisbet (1989) delay-differential equation for giant kelp (Macrocystis pyrifera) populations.

• Equation: $\frac{1}{A(t)}\frac{dA(t)}{dt} = se^{-\alpha\tau}A(t)[1 - a_A A(t-\tau)]_+ - m$
• Features: Includes a time-delayed variable ($A(t-\tau)$), non-linear operators ($[\cdot]_+$), and exponential terms.
• Noise: The model incorporated multiplicative process noise in juvenile settlement and adult mortality ($s \to s e^{\epsilon_1}$, $m \to m e^{\epsilon_2}$), a factor often neglected in prior SR benchmarks.

Case Studies & Data Generation

Six case studies were generated to test specific variables:

1. Deterministic: Symmetric vs. Asymmetric cycles.
2. Stochastic: Asymmetric cycles with Low vs. High process noise.
3. Preprocessing: Both Discrete-time (log-ratio) and Continuous-time (cubic spline derivative) approaches were used to estimate per capita growth rates.
4. Sampling Densities: Data was downsampled to 5, 10, 25, 50, and 100 points per cycle.
5. Predictors: The algorithm was given the true drivers ($A(t)$, $A(t-2)$) plus spurious, autocorrelated variables ($A(t-1)$, $A(t-3)$).

Symbolic Regression Implementation

• Tool: PySR (v1.3.1), a Python library with a Julia backend.
• Protocol: 100 independent runs per case study/density combination (3,000 total searches).
• Fitness Metric: Mean Squared Error (MSE).
• Search Space: 96 semi-independent populations with migration to explore equation space.

Evaluation Workflows

Four workflows were tested to select the "best" equation from the results:

1. Workflow 1 (Subjective): Visual inspection of MSE vs. Complexity; selecting the simplest equation with the largest additive MSE drop.
2. Workflow 2 (Subjective): Visual inspection using $\ln(\text{MSE})$; selecting the simplest equation with the largest multiplicative MSE drop.
3. Workflow 3 (Objective): Using PySR's built-in score (negative discrete log loss change per unit complexity).
4. Workflow 4 (Objective): Bayesian Information Criterion (BIC).

Success Metrics:

• Variable Recovery: Did the selected equation contain only the correct variables ($A(t)$, $A(t-2)$)?
• Equation Recovery: Did the selected equation match the Bence-Nisbet model in functional form and parameters?
• Presence: Did the true model appear anywhere in the top 10 equations (regardless of selection)?

3. Key Results

Sampling Density is the Primary Driver

• Critical Threshold: SR failed to recover the underlying equation at sampling densities below 10 points per cycle.
• High Density Requirement: Even at densities of 25–50 points per cycle, the true model was frequently generated by the algorithm but not consistently selected by the evaluation workflows.
• Performance: Success rates for equation recovery were generally low (<20%) except in high-density (50+ points/cycle) scenarios with high process noise using subjective workflows.

The Role of Process Noise

• Contrary to the intuition that noise hinders recovery, process noise increased success rates.
• Stochastic perturbations expanded the system's exploration of state space, making the data more informative and helping distinguish between dynamically equivalent models.
• High noise levels performed as well as or better than low/noise scenarios.

Equation Generation vs. Selection

• Generation: The algorithm successfully evolved the correct Bence-Nisbet equation in nearly all case studies with sufficient sampling (25+ points/cycle).
• Selection Failure: The true model often ranked first in MSE but was not selected by the workflows.
  • The true model frequently appeared "behind" the Pareto frontier or was overshadowed by more complex equations that were algebraically equivalent but had extraneous constants.
  • Objective workflows (3 and 4) generally performed worse than subjective visual inspections (1 and 2) in identifying the true model.

Variable Identification

• At high sampling densities, the algorithm correctly identified the true variables ($A(t)$ and $A(t-2)$) even when the full equation structure wasn't selected.
• At low densities, spurious variables ($A(t-1)$, $A(t-3)$) were frequently included due to autocorrelation, leading to incorrect variable selection.

Asymmetry and Preprocessing

• Cycle Asymmetry: Had negligible impact on performance; asymmetric cycles performed comparably to symmetric ones.
• Preprocessing: No consistent advantage was found between discrete-time and continuous-time approaches for estimating growth rates.

4. Key Contributions

1. Realistic Benchmarking: This is the first study to rigorously test SR against field-relevant conditions, specifically low sampling density, process noise, and spurious variables, moving beyond idealized laboratory simulations.
2. Decoupling Generation and Selection: The study highlights a critical distinction: SR is often capable of generating the correct equation, but current selection workflows fail to identify it. The bottleneck is not the search algorithm but the post-hoc criteria.
3. Process Noise Benefit: Demonstrated that process noise can be beneficial for equation discovery by increasing data informativeness, challenging the view that noise is purely detrimental.
4. Sampling Density Guidelines: Established that 25–50 points per cycle is likely a minimum threshold for reliable SR application in population dynamics, far exceeding the "hundreds of points" used in previous literature.

5. Significance and Implications

• For Ecologists: The paper provides a "guide" for when SR is appropriate. It suggests that for typical field data (often <10 points/cycle), SR is unlikely to succeed without significant methodological adjustments.
• Methodological Shift: The authors argue that the field must move beyond relying solely on the Pareto frontier and goodness-of-fit metrics (MSE, BIC). Future workflows need structural identifiability criteria or dynamical diagnostics (e.g., checking if the model reproduces key statistical properties of the time series) to select the correct model among dynamically equivalent candidates.
• Future Directions: The study suggests that strategies like irregular sampling, compressive sampling, or state-space formulations (which avoid derivative estimation errors) may be necessary to make SR viable for real-world ecological time series.

In conclusion, while symbolic regression holds the potential to reverse-engineer complex ecological mechanisms, its current application is limited by data sparsity and a lack of robust model selection criteria. Success requires high-density data and a shift in how researchers evaluate the output of these algorithms.