Penalized Likelihood Parameter Estimation for Differential Equation Models: A Computational Tutorial

Parameter estimation connects mathematical models to real-world data and decision making across many scientific and industrial applications. Standard approaches such as maximum likelihood estimation and Markov chain Monte Carlo estimate parameters by repeatedly solving the model, which often requires numerical solutions of differential equation models. In contrast, generalized profiling (also called parameter cascading) focuses directly on the governing differential equation(s), linking the model and data through a penalized likelihood that explicitly measures both the data fit and model fit. Despite several advantages, generalized profiling is relatively rarely used in practice. This tutorial-style article outlines a set of self-directed computational exercises that facilitate skills development in applying generalized profiling to a range of ordinary differential equation models. All calculations can be repeated using reproducible open-source Jupyter notebooks that are available on GitHub.

💡 Research Summary

The paper presents a comprehensive tutorial on generalized profiling (also known as parameter cascading) for estimating parameters in ordinary differential equation (ODE) models. Traditional approaches such as maximum likelihood estimation (MLE) and Markov chain Monte Carlo (MCMC) require solving the ODE at every iteration, which can be computationally expensive and sensitive to numerical truncation errors. Generalized profiling circumvents these issues by introducing smooth trial functions—specifically B‑splines—to approximate the solution of the ODE directly, without ever computing the exact numerical solution.

The authors begin by reviewing B‑spline theory, emphasizing their local support, recursive Cox–de Boor construction, and the ability to represent any spline of a given order as a linear combination of basis functions. Using the Julia package BSplineKit.jl, they demonstrate how to construct cubic B‑splines, fit them to noisy data, and compute their first derivatives efficiently via sparse banded matrices A (function values) and A′ (derivatives). This dual representation enables the formulation of two loss components: a data‑fit term measuring the discrepancy between observed data and the spline, and a model‑fit term measuring the discrepancy between the spline derivative and the right‑hand side of the ODE evaluated at the spline. The total penalized likelihood is L = ‖y – A c‖² + λ ‖A′ c – g(θ, A c)‖², where c are the spline coefficients, θ the ODE parameters, g the ODE right‑hand side, and λ a penalty weight that balances data fidelity against ODE compliance.

Three case studies illustrate the workflow. The first uses Newton’s law of cooling (a linear first‑order ODE) to estimate the initial temperature, ambient temperature, heat‑transfer coefficient, and noise variance from synthetic noisy measurements. The second tackles the logistic growth model, estimating the intrinsic growth rate and carrying capacity. The third addresses a two‑species Lotka‑Volterra competition system, demonstrating how multiple state variables can be handled with separate splines while sharing a common penalty structure. In each example, the authors initialize the spline by interpolating the data, then iteratively update the spline coefficients and the ODE parameters using gradient‑based optimizers (e.g., Gauss‑Newton, BFGS). The penalty parameter λ is selected via an L‑curve or cross‑validation, ensuring that the spline does not overfit the noise yet remains faithful to the underlying dynamics.

Results show that generalized profiling achieves comparable or superior parameter accuracy to MLE/MCMC while requiring far fewer ODE solves, leading to substantial reductions in computational time. Moreover, the method naturally regularizes against overfitting, especially when the data are noisy or sparsely sampled. The authors also draw a clear connection to recent physics‑informed neural networks (PINNs) and biologically‑informed neural networks (BINNs), noting that the data‑loss and physics‑loss terms in those frameworks are mathematically identical to the two components of the penalized likelihood used here. This highlights that generalized profiling provides a solid analytical foundation for PINNs and can inform choices of network architecture, loss weighting, and initialization.

All code is provided as open‑source Jupyter notebooks on GitHub, written in Julia, and relies on widely used packages (BSplineKit.jl, Optim.jl, Distributions.jl). The notebooks are fully reproducible, allowing readers to modify knot placement, spline order, boundary conditions, and penalty weights to explore their effects on estimation. By offering a step‑by‑step, hands‑on approach, the paper lowers the barrier to adopting generalized profiling in applied sciences, ranging from engineering to biology, and positions it as a viable alternative—or complement—to standard likelihood‑based and Bayesian inference methods for ODE models.