Impilict Runge-Kutta based sparse identification of governing equations in biologically motivated systems

TL;DR

This work addresses the challenge of identifying governing differential equations from sparse, noisy data in biological and physical systems. It introduces IRK-SINDy, which couples high-order implicit Runge-Kutta methods with sparse regression, and adds a deep IRK-SINDy variant that uses a neural network to predict IRK stage values for efficient, derivative-free learning. The approach yields parsimonious, interpretable models and shows superior robustness to data scarcity and noise across a suite of benchmarks including linear/cubic oscillators, Lorenz, predator-prey, logistic growth, and FitzHugh–Nagumo, outperforming conventional SINDy and RK4-SINDy. The work advances data-driven model discovery for biology by leveraging A-stable IRKs and neural stage predictors, with plans to extend to other high-order implicit methods and to provide code.

Abstract

Identifying governing equations in physical and biological systems from datasets remains a long-standing challenge across various scientific disciplines, providing mechanistic insights into complex system evolution. Common methods like sparse identification of nonlinear dynamics (SINDy) often rely on precise derivative estimations, making them vulnerable to data scarcity and noise. This study presents a novel data-driven framework by integrating high order implicit Runge-Kutta methods (IRKs) with the sparse identification, termed IRK-SINDy. The framework exhibits remarkable robustness to data scarcity and noise by leveraging the lower stepsize constraint of IRKs. Two methods for incorporating IRKs into sparse regression are introduced: one employs iterative schemes for numerically solving nonlinear algebraic system of equations, while the other utilizes deep neural networks to predict stage values of IRKs. The performance of IRK-SINDy is demonstrated through numerical experiments on benchmark problems with varied dynamical behaviors, including linear and nonlinear oscillators, the Lorenz system, and biologically relevant models like predator-prey dynamics, logistic growth, and the FitzHugh-Nagumo model. Results indicate that IRK-SINDy outperforms conventional SINDy and the RK4-SINDy framework, particularly under conditions of extreme data scarcity and noise, yielding interpretable and generalizable models.

Figures (9)

• Figure 1: Overview of the IRK-SINDy framework: a. For each benchmark problem, we perform measurements that incorporate noise and, thereafter form a dataset. Our objective is to construct a model that is parsimonious, interpretable, and possesses generalizability, capable of accurately forecasting reference dynamics. b. Given an appropriate initial guess (e.g., $X(t_{k})$), the stage values of the IRKs are approximated by solving the system of nonlinear equations \ref{['eqrkb']} through iterative schemes. In this context, we employ two iterative approaches: (i) fixed point iteration and (ii) Newton's method. c. With the stage values estabilished the subsequent step values are computed according to eq.\ref{['Firk']}. This computational process is depicted as the systematic IRK network. d. Within this structured representation of IRK-SINDy, the dataset is classified into two categories: forward and backward, followed by the formation of a symbolic features library comprising candidate nonlinear functions. To solve a nonlinear sparse regression problem using the forward and backward predictions illustrated in (b) and (c), an IRK step is applied, and the loss function is minimized by choosing a suitable optimizer. Following a certain number of epochs, a sparsity-promoting algorithm is employed. Finally, every non-zero element in the coefficient matrix $\xi^{*}$ signifies an active term within the feature library, thereby representing the resultant discovered model.
• Figure 2: Overview of the deep IRK-SINDy framework: a. The dataset is prepared for the purpose of training the neural network. b. The inputs to the neural network are assigned into two distinct variables: time and state variables. The neurons located in the output layer of the network are partitioned into $s$ segments, each containing $d$ neurons. The $i$'th segment predicts the $d$ stage values corresponding to the $\chi_{i}$. c. Through the process of forward propagation within the DNN, the stage values are predicted, and these predictions are subsequently employed in the IRK steps, i.e. eq.\ref{['Firk']}, facilitating both forward and backward predictions. d. By comparing the predictions against the data, the loss is computed, followed by the optimization step. Upon reaching a specified number of epochs, at which point the loss is sufficiently minimized, the sparsity-promotion algorithm is exclusively applied to the coefficient matrix $\xi$. Finally, the non-zero coefficients of the $\xi$ denote the active terms in the nonlinear feature library.
• Figure 3: Linear damped oscillator: Comparing identified models under various levels of data scarcity with reference model. a. Data, b. sample size $m=801$, c. sample size $m=201$, d. sample size $m=51$, e. sample size $m=31$.
• Figure 4: Linear damped oscillator: Comparing the response of identified models under various noise levels in measurements with reference model. a. Noisy data, b. noise level $\sigma = 0.01$, c. nise level $\sigma = 0.04$, d. noise level $\sigma=0.08$, e. noise level $\sigma=0.16$.
• Figure 5: Cubic damped oscillator: Comparing identified models under various levels of data scarcity with reference model. a. Phase portraits, b. coefficient matrices for $m\in\{ 801, 401, 101, 51 \}$. IRK-SINDy provides a more parsimonious and generalizable model compared to RK4-SINDy and Conv-SINDy.