Kinetics-Informed Neural Networks

How do you figure out the rules governing a system when all you have is noisy measurements? This is the inverse problem, and it arises everywhere: a chemist measuring reaction rates to infer activation energies, a biologist fitting population dynamics to recover growth constants, an engineer calibrating a turbulence model from wind-tunnel data. In each case, the underlying physics is described by differential equations whose parameters are unknown, and the goal is to recover those parameters from observations alone.

Neural networks can approximate any continuous function, which makes them powerful interpolators. But without constraints, they overfit to noise and produce predictions that violate basic physical laws: conservation of mass, non-negative concentrations, and thermodynamic consistency. Physics-informed neural networks (PINNs) [3] address this by encoding the governing differential equations directly into the training objective. The network must simultaneously fit the data and satisfy the physics. This dual constraint acts as a powerful regularizer: it prevents overfitting, enables extrapolation beyond the training window, and, crucially, allows the unknown parameters of the differential equations to be recovered alongside the solution itself.

We developed Kinetics-Informed Neural Networks (KINNs) [1] as a specialized PINN framework for chemical kinetics, embedding stoichiometric constraints, conservation laws, and mass-action rate expressions into the neural network architecture. KINNs operate in two complementary modes: the forward problem, where known rate constants yield smooth concentration profiles satisfying the ODE system; and the inverse problem, where the rate constants themselves become trainable parameters recovered from noisy data. We then reformulated the inverse problem using maximum-likelihood estimation (MLE) [2], eliminating the need for manual hyperparameter tuning and enabling automatic uncertainty quantification on every recovered parameter.

The interactive tool below lets you experience all of this in your browser. Configure the network architecture, toggle between forward and inverse modes, add noise, switch between the original and MLE-based formulations, and watch training unfold in real time. Every computation runs client-side in JavaScript, implementing the same mathematical framework described in the original publications [1], [2].

The Forward Problem

The forward problem consists of solving the kinetic ODE system given known rate constants ($k_1$, $k_{-1}$, $\ldots$) and initial conditions. The neural network $\hat{\mathbf{x}}(t; \boldsymbol{\theta})$ serves as a surrogate approximator for the concentration profiles over time. Under the mean-field approximation, the governing ODE relates species concentrations to reaction rates through the stoichiometric matrix:

$$\frac{d\mathbf{c}}{dt} = \mathbf{M} \cdot \mathbf{r}(\mathbf{c}, \boldsymbol{\theta})$$

For the simple reversible reaction $A + B \rightleftharpoons C$, the rate expression follows from power-law kinetics:

$$A + B \rightleftharpoons C, \quad r = k_1[A][B] - k_{-1}[C]$$

The training objective is formulated as the minimization of two complementary terms. The data loss ($j_d$) measures the mean squared error between the network output and observed concentrations at sampled time points. The physics loss ($j_m$) penalizes the residual of the ODE system at collocation points throughout the time domain, enforcing the constraint $\tfrac{d\hat{\mathbf{x}}}{dt} = \mathbf{f}(\hat{\mathbf{x}}, \mathbf{k})$, where $\mathbf{f}$ encodes the mass-action kinetics.

The physics loss acts as a regularizer: even with very few data points, the network cannot learn concentration profiles that violate conservation laws or kinetic constraints. This is fundamentally different from a purely data-driven fit, which may overfit to noise or produce physically impossible trajectories. Structural boundary condition operators further guarantee that the initial conditions are satisfied for any set of network weights, as described in the original KINNs formulation [1].

The Inverse Problem

Toggle the mode to Inverse in the tool above. In this configuration, the rate constants are unknown and become additional trainable parameters alongside the network weights. The optimizer must simultaneously learn the concentration profiles and recover the rate constants from data alone. Formally, we combine the residual between the surrogate approximator and observed data ($j_d$) with the physics regularization term ($j_m$):

$$\mathcal{L} = \frac{1}{d}\sum_{i=1}^{d} \|\dot{\mathbf{x}}_{NN}(t_i) - \mathbf{f}(\mathbf{x}_{NN}(t_i), \boldsymbol{\theta})\|^2 + \alpha \cdot \frac{1}{d}\sum_{i=1}^{d} \|\mathbf{x}_{NN}(t_i) - \tilde{\mathbf{x}}_i\|^2$$

This formulation addresses the inverse problem at the heart of experimental kinetics: given time-series measurements of species concentrations, determine the underlying rate law parameters. Traditional approaches require solving a stiff ODE system at every optimization step. KINNs circumvent this by having the neural network itself serve as the ODE solution, with automatic differentiation providing the time derivatives directly. The kinetic parameters are trained in an exponential mapping ($k = \exp(p)$) to ensure non-negativity and to place them on a similar scale as the network weights, enabling simultaneous optimization [1].

Observe the Parameter Convergence chart as training progresses. The estimated rate constants converge toward the true values (shown as dashed lines). The speed and accuracy of convergence depend on the noise level, the amount of available data, and the relative weighting between data fidelity and physics constraints, as controlled by the parameter $\alpha$.

Self-Adapting Uncertainty Through Error Propagation

In the original KINNs formulation, the total loss is a weighted sum: $\mathcal{L} = j_m + \alpha \cdot j_d$. The regularization hyperparameter $\alpha$ conveys the ratio between the variances of the data residuals (observed states) and the physics residuals (state derivatives) [1]. Selecting an appropriate value of $\alpha$ is notoriously difficult; if $\alpha$ is too small the physics constraints dominate and the data is ignored, while if $\alpha$ is too large the network overfits to noisy observations at the expense of physical consistency.

To address this hyperparameterization issue, we reformulated the inverse problem in terms of maximum-likelihood estimators (MLE), giving rise to robust-KINNs (rKINNs) [2]. Instead of a fixed scalar $\alpha$, rKINNs learn full covariance matrices $\boldsymbol{\Sigma}_x$ (data) and $\boldsymbol{\Sigma}_f$ (physics) that are continuously updated based on the incumbent mismatch between observations and network interpolations. The loss function becomes the negative log-likelihood:

$$\mathcal{L}_{MLE} = \frac{1}{d}\sum_{i=1}^{d} \left[\boldsymbol{\varepsilon}_{\dot{x},i}^T \boldsymbol{\Omega}_{\dot{x},i} \boldsymbol{\varepsilon}_{\dot{x},i} + \boldsymbol{\varepsilon}_{x,i}^T \boldsymbol{\Omega}_x \boldsymbol{\varepsilon}_{x,i}\right]$$

A key insight of the rKINNs formulation is that the physics-residual covariance is not independent of the measurement uncertainty. It propagates from the interpolation residuals through the Jacobians of the kinetic model at each time point:

$$\boldsymbol{\Sigma}_{\dot{x},i} = \mathbf{J}_x^{(i)} \boldsymbol{\Sigma}_x \mathbf{J}_x^{(i)T} + \mathbf{J}_p^{(i)} \boldsymbol{\Sigma}_p \mathbf{J}_p^{(i)T}$$

The relative weighting between data and physics thus emerges naturally from the data itself, and the hyperparameter $\alpha$ is no longer needed. Furthermore, the MLE construction provides a statistical backdrop for the interplay between data interpolation and physical model agreement, which enables parameter uncertainty estimation based on second-order analysis (Hessian-based uncertainty quantification) [2].

To observe this in practice, switch the algorithm to Auto α in the tool and note the α auto badge. The alpha slider is disabled because the covariance matrices handle the weighting automatically. The Covariance Evolution chart shows how $\text{tr}(\boldsymbol{\Sigma}_x)$ and $\text{tr}(\boldsymbol{\Sigma}_f)$ evolve during training.

Extrapolation Through Physical Constraints

One of the most compelling advantages of embedding physical constraints into the neural network is the ability to extrapolate beyond the training domain. A purely data-driven model will diverge rapidly outside the distribution of observed data, whereas a KINNs model is constrained by the ODE system at every collocation point, including regions where no training data is available.

To observe this, reduce the training region to $t \in [0, 2.5]$ using the Training Region sliders in the configuration panel. Then examine the predictions in the $t \in [2.5, 5.0]$ range. With physics constraints active, the network continues to produce physically consistent extrapolations that respect the governing differential equations. Without them (set $\alpha = 0$), the predictions quickly become nonsensical.

The mechanism behind this extrapolation is worth emphasizing. The data loss $j_d$ is evaluated only at observation times within the training region, but the physics loss $j_m$ is evaluated at collocation points spanning the entire time domain, including regions where no measurements exist. The neural network basis must satisfy $\tfrac{d\hat{\mathbf{x}}}{dt} = \mathbf{M} \cdot \mathbf{r}(\hat{\mathbf{x}}, \mathbf{k})$ everywhere, not just where data is available. This is what allows the network to produce physically consistent predictions beyond the training window: the ODE constraint propagates information from the data-rich region into the data-free region through the governing equations.

This capability is particularly relevant for experimental kinetics, where measurements are expensive and time-limited. KINNs can leverage sparse early-time data to predict late-time behavior, which is precisely the scenario where traditional regression methods struggle most. More broadly, it illustrates the benefit of physics-informed approaches: the inductive bias provided by the kinetic model enables generalization that would be unattainable with flexible but unconstrained function approximators.

Generalization to Coupled Nonlinear Systems

While we demonstrate the framework on chemical reaction networks, the underlying mathematical structure is far more general. The essential ingredients are: a system of coupled nonlinear ODEs, a fixed linear map (the stoichiometric matrix $\mathbf{M}$) relating state derivatives to a nonlinear rate vector, partial observability of the state, and unknown parameters to be recovered from noisy measurements.

This structure appears across many domains. In systems biology, gene regulatory networks and metabolic pathways obey analogous stoichiometric constraints with unknown kinetic parameters. In pharmacokinetics, compartmental models describe drug absorption and elimination through coupled ODEs with rate constants inferred from blood concentration data. In atmospheric chemistry, hundreds of coupled radical reactions govern ozone depletion and pollutant formation, with rate constants that must be extracted from spectroscopic measurements. In reactor engineering, microkinetic models with dozens of elementary steps and surface intermediates require simultaneous estimation of activation energies from transient experiments.

The key generalization is this: any high-dimensional process governed by coupled nonlinear ODEs, where the state evolves on a manifold defined by a linear constraint structure (conservation laws, stoichiometric balance, site balance), and where parameters must be recovered from partial, noisy observations, falls within the scope of this framework. The neural network provides a universal basis for the solution manifold. The stoichiometric SVD decomposes the problem into range and null space components, separating the dynamics from the conserved quantities. The MLE formulation propagates measurement uncertainty through the model Jacobians, automatically weighting the data-physics trade-off at every point in time and in every direction of the stoichiometric subspace.

The interactive tool above implements this complete pipeline for simple reaction systems. The same mathematical machinery, implemented in JAX with automatic differentiation and XLA compilation [1], [2], [4], scales to the high-dimensional microkinetic models encountered in heterogeneous catalysis, where the parameter space spans hundreds of rate constants and the state includes both gas-phase and surface species.