Learning physically consistent differential equation models from data using group sparsity

We propose a statistical learning framework based on group-sparse regression that can be used to (i) enforce conservation laws, (ii) ensure model equivalence, and (iii) guarantee symmetries when learning or inferring differential-equation models from data. Directly learning interpretable mathematical models from data has emerged as a valuable modeling approach. However, in areas such as biology, high noise levels, sensor-induced correlations, and strong intersystem variability can render data-driven models nonsensical or physically inconsistent without additional constraints on the model structure. Hence, it is important to leverage prior knowledge from physical principles to learn biologically plausible and physically consistent models rather than models that simply fit the data best. We present the group iterative hard thresholding algorithm and use stability selection to infer physically consistent models with minimal parameter tuning. We show several applications from systems biology that demonstrate the benefits of enforcing priors in data-driven modeling.

Figure 1

Core module of the JAK-STAT signaling pathway. The hormone EPO binding to the EpoR receptor results in activation of the receptor [activated form with concentration $c\left(t\right)]$ by transphosphorylation of JAK2 and subsequent tyrosine phosphorylation (P) of JAK2 and the EpoR cytoplasmic domain. Phosphotyrosine residues 343 and 401 in EpoR mediate recruitment of monomeric STAT-5 (concentration $x_1$). Upon receptor recruitment, monomeric STAT-5 is tyrosine phosphorylated ($x_2$), dimerizes ($x_3$), and translocates to the nucleus ($x_4$), where it binds to the promoters of target genes, is dephosphorylated, and is exported again to the cytoplasm [52]. The inset plot shows an experimentally measured time course of EpoR activation (data from [53]).

Figure 2

Dictionary design and simulated data. (a) Dictionary construction and coefficient grouping. The identical dictionaries ${\mathit{\Theta }}_i={\mathit{\Theta }}_b$ for each $x_i, i=1,2,3,4$, are stacked in a block-diagonal matrix for joint learning. Vertical lines indicate the coefficient groups $g_1,...,g_4$, corresponding to the four biochemical processes named in the legend. (b) Time-series data for different concentrations in the JAK-STAT pathway obtained by numerically integrating the deterministic ODE model in Eqs. (8) using the ode45 matlab solver and adding 10% Gaussian noise ($\sigma =0.1$).

Figure 3

Inferring JAK-STAT signaling models from noisy data. (a) Stability plot using grouping based on mass conservation. In the gray shaded range of $\lambda$ values, stability selection with ${\pi }_{\text{th}}=0.8$ identifies the correct model. The red solid lines show the behavior of the true components of the ODE model; the black dashed lines are all other $P=76$ dictionary terms. (b) Stability plot without grouping. There is no value of $\lambda$ for which the true model is found. In both (a) and (b) Gaussian noise with $\sigma =0.1$ is added to the simulated data before inference and $n=200$ time points for each component $x_i, i=1,2,3,4$, are used. (c) Achievability plot for model selection with mass conservation prior. (d) Achievability plot for model selection without mass conservation prior. In (c) and (d) the success probabilities of inferring the correct model over 20 independent trials with different noise realizations and different random data subsampling are shown as a function of the number of data points used in each trial. Colored bands are Bernoulli standard deviations for different amounts of noise added to the simulated data prior to inference (see the legend).

Figure 4

Simulated data used to learn spatiotemporal models of one-dimensional advection-diffusion dynamics. Visualization of the data is shown for (a) $u(x,t)$ and (b) $v(x,t)$ with 15% additive Gaussian noise ($\sigma =0.15$). Spatial and temporal discretization use 256 and 200 regularly spaced grid points, respectively. The solution is obtained via spectral differentiation and fourth-order Runge-Kutta time integration. The diffusion constants of the species are $D_u=0.25$ and $D_v=0.50$ in nondimensional units. The equations are solved with periodic boundary conditions in the domain $x\in [-\pi ,\pi )$ of length $L=2\pi$ over the time horizon $t\in [0,3]$ with initial conditions $u(x,t=0)=cos\left(\frac{2\pi x}{L}\right)$ and $v(x,t=0)=-cos\left(\frac{2\pi x}{L}\right)$ for species $u$ and $v$, respectively.

Figure 5

Inferring advection-diffusion dynamics with unknown spatially varying velocity field. (a) Stability plot with groups to encode both spatially varying coefficients and model equivalence between the species. The gray shaded region is the range of $\lambda$ for which model selection with ${\pi }_{\text{th}}=0.8$ identifies the correct model. (b) Stability plot with groups only to encode spatially varying coefficients, but no grouping for model equivalence. (c) Stability plot with no groupings at all. In (a)–(c) the red solid lines correspond to the true components of the PDE for the field $u$, with symbols referring to the differential operators as given in the legends. The dictionary block size is $n=75$ and $p=15$ with $p_b=10$ blocks and 15% Gaussian noise ($\sigma =0.15$) added to the simulated data. (d) Achievability plot for model selection using both priors for different levels of noise in the data. Each point is averaged over 20 independent trials. The colored bands correspond to the Bernoulli standard deviation.

Figure 6

Simulated data used to learn reaction-diffusion dynamics. Visualization of the two-dimensional concentration fields (a) $u(x,y)$ and (b) $v(x,y)$ is shown at time $t=7.5$ from numerical solution of the model in Eqs. (12) with $D_u=D_v=0.1$ and 10% additive Gaussian noise ($\sigma =0.1$). The solution is obtained via spectral differentiation in the domain $(x,y)\in {[-10,10]}^2$ and fourth-order Runge-Kutta time integration with time step size 0.05 on a Cartesian grid of $128\times 128$ points with initial conditions $u(x,y,0)=\text{tanh}\{\sqrt{x^2+y^2}cos[3\angle (x+iy)-\sqrt{x^2+y^2}]\}$ and $v(x,y,0)=\text{tanh}\{\sqrt{x^2+y^2}sin[3\angle (x+iy)-\sqrt{x^2+y^2}]\}$.

Figure 7

Inferring reaction-diffusion models from noisy spatiotemporal data. Stability plots (a) with and (b) without symmetry priors for noise level $\sigma =0.1$. Achievability plots (c) with and (d) without symmetry priors for different noise levels in the data (legends). We show the stability plots for ${\lambda }_{\text{min}}=0.01{\lambda }_{\text{max}}$ to illustrate the difference between (a) and (b), but still run our algorithm with ${\lambda }_{\text{min}}=0.1{\lambda }_{\text{max}}$. Each point is averaged over 20 independent trials. The colored bands correspond to the Bernoulli standard deviation.

Figure 8

Relative errors in the coefficients inferred for the JAK-STAT pathway reactions. The plots show the relative errors $\frac{|\xi -{\xi }^{*}|}{{\xi }^{*}}$ (vs ground truth ${\xi }^{*}$) in the reaction rate estimates of the JAK-STAT pathway as inferred by the GIHT algorithm for different noise levels $\sigma$. Symbols show estimated means with bars indicating estimation standard deviations over 20 independent trials. The ground-truth values are $k_1^{\pm }=0.021, k_2^{\pm }=2.46, k_3^{\pm }=0.2066$, and $k_4^{\pm }=0.10658$. The closed and open symbols correspond to the independently estimated rate constants of different signs, which should be identical.

Figure 9

Spatially varying velocity field and its gradient for the advection-diffusion example. The plots show the estimates for the latent spatially varying velocity $c\left(x\right)$ (left column) and its gradient ${\partial }_{x}c\left(x\right)$ (right column) from the GIHT algorithm. The rows correspond to the inference from data with different noise levels $\sigma$ (shown also in the legend). Symbols show estimated means with bars indicating estimation standard deviations over 20 independent trials. Black solid lines are the ground truth.

Figure 10

Relative errors in the coefficient estimation for the $\lambda \text{−}\omega$ reaction-diffusion system. The plots show the relative errors $\frac{|\xi -{\xi }^{*}|}{{\xi }^{*}}$ (vs ground truth ${\xi }^{*}$) in the GIHT estimates of reaction coefficients and diffusion constants for the species $u$ (left) and $v$ (right) as a function of the noise level $\sigma$ in the data. The ground-truth coefficients for the species $u$ and $v$ are as given in Eq. (12).

Figure 11

Reconstructed spatially varying velocity field and its gradient for the advection-diffusion example with additional second-order smoothness prior. The plots show the estimates for the latent spatially varying velocity $c\left(x\right)$ (left column) and its gradient ${\partial }_{x}c\left(x\right)$ (right column) when imposing an additional smoothness prior over each group recovered by stability selection using GIHT. The rows correspond to data with different noise levels $\sigma$ (shown in the legend). For $\sigma =0$ (top row), we use a smoothness regularization ${\lambda }_{\text{f}}=1$, and for $\sigma =0.05$ and 0.1 (middle and bottom rows), we use ${\lambda }_{\text{f}}=20$. Black solid lines are the ground truth.