Discovering stochastic dynamical equations from biological time series data

TL;DR

The paper tackles the inverse problem of inferring governing stochastic differential equations from biological time-series data. It combines drift/diffusion estimation from jump moments with sparse regression to learn concise, symbolic models for $\mathbf f$ and $\mathbf G$, using cross-validated sparsity to prevent overfitting. Diagnostic tools for noise structure, model self-consistency through simulated data, and validation on left-out data bolster reliability, and the authors provide PyDaDDy, an open-source package, to implement the workflow. Demonstrations on synthetic SDEs with identical steady-state distributions and real data from schooling fish polarization and confined cell migration illustrate the method's ability to uncover noise-driven vs deterministic dynamics and to recover interpretable stochastic models. The work advances mechanistic, data-driven modeling of biological systems by enabling interpretable SDE discovery from time-series with diagnostic safeguards and accessible software.

Abstract

Theoretical studies have shown that stochasticity can affect the dynamics of ecosystems in counter-intuitive ways. However, without knowing the equations governing the dynamics of populations or ecosystems, it is difficult to ascertain the role of stochasticity in real datasets. Therefore, the inverse problem of inferring the governing stochastic equations from datasets is important. Here, we present an equation discovery methodology that takes time series data of state variables as input and outputs a stochastic differential equation. We achieve this by combining traditional approaches from stochastic calculus with the equation-discovery techniques. We demonstrate the generality of the method via several applications. First, we deliberately choose various stochastic models with fundamentally different governing equations; yet they produce nearly identical steady-state distributions. We show that we can recover the correct underlying equations, and thus infer the structure of their stability, accurately from the analysis of time series data alone. We demonstrate our method on two real-world datasets -- fish schooling and single-cell migration -- which have vastly different spatiotemporal scales and dynamics. We illustrate various limitations and potential pitfalls of the method and how to overcome them via diagnostic measures. Finally, we provide our open-source codes via a package named PyDaDDy (Python library for Data Driven Dynamics).

Figures (6)

• Figure 1: Overall schematics of the PyDaDDy package. PyDaDDy takes as input uniformly sampled 1D or 2D time series, and computes the drift and diffusion components from the time series. Several functions are provided to visualize data as time series or histograms, fit drift and diffusion functions, diagnose whether underlying assumptions for drift-diffusion estimation are met, and to export data.
• Figure 2: Sparse regression on the extracted instantaneous jump moments reliably reconstructs governing equations from simulated data, for a wide range of dynamics. The mathematical expressions for the actual and estimated drift and diffusion functions are shown in the shaded boxes on the top (grey-shaded) and bottom (red-shaded), respectively. Row i shows a sample of the time series $x(t)$ of the models Eqs \ref{['eq:m1']}-\ref{['eq:m4']}, respectively, across four columns (A-D), while row ii shows the corresponding histograms of $x(t)$. The arrows in histograms in row ii mark the deterministic stable equilibria of the respective systems. Notice that the histograms in (A) and (B) are nearly identical, despite the underlying SDEs being very different. Similarly, the histograms in (C) and (D) are similar. (A) A deterministic unimodal system (from Eq \ref{['eq:m1']}. (B) A noise-induced unimodal state (from Eq. \ref{['eq:m2']}). (C) A deterministic bimodal system from Eq. \ref{['eq:m3']}), with noise facilitating transitions between two deterministic states. (D) A noise-induced bimodal state (from Eq. \ref{['eq:m4']}). Rows iii and iv compare the ground-truth (black) drift and diffusion functions with the estimated ones. The functions estimated using sparse regression are shown as dark red dashed lines. In all cases, the proposed method accurately recovers the drift and diffusion functions.
• Figure 3: Data driven SDE discovery for polarization dynamics in fish schools; data from jhawar2020fish. (A) Experimental setup and individual fish trajectories (in different colours) are extracted from video recordings of fish swimming in a tank (from jhawar2020fish.). (B) From the individual trajectories, a time series of the polarization vector, $\bm$, is computed. The trajectories shown are the $x$- (red) and $y$- (blue) components of $\bm$. (C, D) Histograms of the polarization vector $\bm$ and the net polarization, $|\bm|$, respectively. (E) The drift function discovered by the SDE discovery procedure. The $x$-component of the function, $f_x$ is shown. The surface plot shows the fitted drift function, and the points show the binwise averaged estimates. Inset shows a slice of this function along the $y=0$ plane, showing a single stable equilibrium at $x = 0$ (F) The discovered diffusion function. The inset shows a slice along the $y=0$ plane. The diffusion is maximum at $\bm = 0$, and decreases outwards. In (E) and (F), the colours redundantly represent $f_x$ and $G_{xx}$, respectively and are added only for visual clarity.
• Figure 4: Diagnostics of the discovered SDE for fish polarization dynamics. (A, B) Noise diagnostics. (A) The histogram of the residual noise $r$ is computed based on the discovered SDE. The $r_x$ and $r_y$ marginals are shown as red and blue. These should match a standard normal distribution, shown in black. (B) The autocorrelation functions of $r_x$ and $r_y$. The autocorrelation should decay within one sampling time step. (C, D) Model diagnostics. (C) Histogram of $|\bm|$ from the original time series (black) compared to that of a simulated time series generated from the discovered SDE (red). (D) Autocorrelation of the $m_x$ component of the real (black) and simulated (red) time series.
• Figure 5: Discovering SDEs governing the dynamics of confined migration of cancer cells; data from bruckner2019stochastic (A) The experimental setup, showing a micropattern with two landing pads, and the cell hopping between them (from bruckner2019stochastic). (B) Example time series of cell position. (C) Joint histogram of the position $x$ and velocity $v$ of the cell trajectories. (D, E) Marginal histograms of $x$ and $v$. (F, G) The discovered SDE for the dynamics of $v$. (F) The drift function, $f(x, v)$. The surface plot shows the fitted drift function, and the points show the binwise averaged estimates. The inset shows the direction field of the deterministic dynamics, with the background colour denoting $f$. (G) The diffusion function, $g(x, v)$. The inset shows the diffusion as a colour map. Since $g$ is relatively higher near the limit cycle, the actual observed cycles have slightly larger amplitude. In (F) and (G), the colours redundantly represent $f_x$ and $G_{xx}$, respectively and are added only for visual clarity.