Sensitivity Approximation by the Peano-Baker Series

TL;DR

The paper tackles the challenge of efficiently computing the sensitivity matrix S for parameter-dependent ODEs in high-dimensional, potentially stiff settings where forward/adjoint methods are costly. It develops a PBS-based representation using the state-transition matrix Φ and constructs numerical approximations via a truncated Peano-Baker series combined with the trapezoidal rule, yielding a PBS formula; to address stiffness, it introduces PBSR (refined PBS) and Exp (exponential) variants. Theoretical results prove a global error of O(Δt_max^2) under standard regularity, with PBSR maintaining this order under mild time-derivative bounds. Numerical experiments on systems biology models (PKA, CaMKII), a random linear system, and Chua’s circuit demonstrate that PBSR achieves substantial accuracy and speed-ups over forward sensitivity, underscoring its practical impact for uncertainty quantification and Fisher-information-based analyses in large-scale ODE models. The work offers a scalable framework for incorporating sensitivity information into Bayesian inference and related tasks, potentially enabling faster convergence in MCMC methods that rely on sensitivity-derived metrics.

Abstract

In this paper we develop a new method for numerically approximating sensitivities in parameter-dependent ordinary differential equations (ODEs). Our approach, intended for situations where the standard forward and adjoint sensitivity analyses become too computationally costly for practical purposes, is based on the Peano-Baker series from control theory. Using this series, we construct a representation of the sensitivity matrix $\mathbf{S}$ and, from this representation, a numerical method for approximating $\mathbf{S}$. We prove that, under standard regularity assumptions, the error of our method scales as $\mathcal{O}(Δt^2_{\text{max}})$, where $Δt_{\text{max}}$ is the largest time step used when numerically solving the ODE. We illustrate the performance of the method in several numerical experiments, taken from both the systems biology setting and more classical dynamical systems. The experiments show the sought-after improvement in running time of our method compared to the forward sensitivity approach. In experiments involving a random linear system, the forward approach requires roughly $\sqrt{n}$ longer computational time, where $n$ is the dimension of the parameter space, than our proposed method.

Figures (8)

• Figure 1: Error propagation graph.
• Figure 2: Relative error of the sensitivity matrix, against time step, for the PBSR (solid line) and the Exp (dashed line) algorithms for the PKA model. The simulations are run over the time span $[0,600]$. The gray circles superimposed on the solid line indicate the time steps at which within the PBSR algorithm the exponential formula \ref{['solutionAtEquilibrium']} is used instead of the PBS formula \ref{['eq:Approx']} because the system is approximately at equilibrium; the small red triangles show when the exponential formula \ref{['solutionAtEquilibrium']} is used instead of the PBS formula \ref{['eq:Approx']} because the number of sub-intervals $n_{int}$ in the PBSR refinement is too large.
• Figure 3: Relative error of the sensitivity matrix, against time step, for the PBSR (solid line) and the Exp (dashed line) algorithms for the CaMKII model. The simulations are run over the time span $[0,600]$. The gray circles superimposed on the solid line indicate the time steps at which within the PBSR algorithm the exponential formula \ref{['solutionAtEquilibrium']} is used instead of the PBS formula \ref{['eq:Approx']} because the system is approximately at equilibrium; the small red triangles show when the exponential formula \ref{['solutionAtEquilibrium']} is used instead of the PBS formula \ref{['eq:Approx']} because the number of sub-intervals $n_{int}$ in the PBSR refinement is too large.
• Figure 4: The estimated relative error $\bar{\varkappa}$ against time for the PKA model. The simulations are run over the time span $[0,60000]$. The blue and red shaded areas both indicate that within the PBSR method the exponential formula \ref{['solutionAtEquilibrium']} is used: blue because the system is close to a (stable) steady state, red because refinement is numerically prohibitive. The error estimate $\bar{\varkappa}$ uses small parameter perturbations $\delta_{\mathbf{p}}$ to test how well the sensitivity matrix predicts the changes in the solution (averaged over $N=100$ small, random parameter perturbations).
• Figure 5: The estimated relative error $\bar{\varkappa}$ against time for the CaMKMIIs model. The simulations are run over the time span $[0,600]$. The blue and red shaded areas both indicate that within the PBSR method the exponential formula \ref{['solutionAtEquilibrium']} is used: blue because the system is close to a (stable) steady state, red because refinement is numerically prohibitive. The error estimate $\bar{\varkappa}$ uses small parameter perturbations $\delta_\mathbf{p}$ to test how well the sensitivity matrix predicts the changes in the solution (averaged over $N=100$ small, random parameter perturbations).

Theorems & Definitions (26)

• Definition 1
• Theorem 3.1
• proof
• Proposition 3.2
• Theorem 3.3
• proof
• Remark 1
• Remark 2
• Corollary 3.4
• Remark 3