0pt captionlabel systems biology and neuroscience

# UQSA - an R-package for uncertainty quantification and sensitivity analysis
for biochemical reaction network models

Andrei Kramer  Science for Life Laboratory, Department of Neuroscience,
Karolinska Institute, Solna, Sweden  Have contributed equally to this work as
first authors.  To whom correspondence should be addressed.  Federica
Milinanni  Department of Mathematics, KTH Royal Institute of Technology,
Stockholm, Sweden.  Have contributed equally to this work as first authors.
Pierre Nyquist  Department of Mathematics, KTH Royal Institute of Technology,
Stockholm, Sweden.  Alexandra Jauhiainen  Early Biometrics and Statistical
Innovation, Data Science & AI, BioPharmaceuticals R&D, AstraZeneca,
Gothenburg, Sweden  Have contributed equally to this work as last authors.
Olivia Eriksson  Science for Life Laboratory, School of Electrical Engineering
and Computer Science, KTH Royal Institute of Technology, Sweden  Have
contributed equally to this work as last authors.  To whom correspondence
should be addressed.

###### Abstract

Summary: We present an R-package developed for modeling of biochemical
reaction networks that especially includes tools for uncertainty
quantification and sensitivity analysis. Estimating parameters and quantifying
their uncertainty, especially in relationship with prediction uncertainty, is
required for data-driven systems biology modeling. The methods we chose
require efficient sampling from high-dimensional, correlated parameter
distributions. We have developed UQSA (Uncertainty Quantification and
Sensitivity Analysis) to be fast for the problem class and work well with
other tools for modelling. We aim for simplicity, and part of that is our use
of the SBtab format for the unified storage of model and data. Our tool-set is
modular enough, that parts can be replaced. We use intermediate formats that
are not hidden from the user to make this feasible. UQ is performed through
Markov chain Monte Carlo (MCMC) sampling in an Approximate Bayesian
Computation (ABC) setting. This can be followed by a variance-decomposition
based global sensitivity analysis. If needed, complex parameter distributions
can be described, evaluated, and sampled from, with the help of Vine-copulas
that are available in R. This approach is especially useful when new
experimental data become available, and a previously calibrated model needs to
be updated. An additional feature is that complicated experimental input data
can be easily fitted to corresponding model functions. This is important for
e.g. neuroscience applications, where transient time dependent input (spike-
trains) are common.

Availability and Implementation: The main part of the code is written in R, as
well as the user facing functions. R is a high-level language and allows the
use of sophisticated statistical methods. The ode solver is written in C in
order to be computationally efficient (from the GSL, odeiv2). We use the SBtab
table format for Systems Biology projects [20] for the model description as
well as the calibration data and an event system to be able to model
complicated transient input. The code has been tested on a one node 256 core
computer cluster, but smaller examples are included that can be run on a
laptop. The source code is hosted on github.com/icpm-kth/uqsa.

##  1 Introduction

UQSA has the goal to combine the preexisting statistics tools of R with
functions specific to the modeling of biological reaction networks. This R
package is well suited for _systems biology_ projects as it uses the SBtab
table format [20] for the model description and experimental data. It also has
functionalities to represent complicated transient input, like spike trains
within neuroscience. With _uncertainty quantification_ (UQ) we refer to any
method that characterizes the full region of the parameter space that result
in a good fit to experimental data. While (global) _sensitivity analysis_ (SA)
refers to the attribution of the uncertainty in outputs to the uncertainty in
the parameters over their entire (joint) probability distribution. In the
first case (UQ) quantitative experimental data is of utmost importance, in the
second case, experimental data is not necessarily needed.

Biochemical reaction network models within systems biology, in neuroscience
and elsewhere, have a tendency to be over-parameterized111for good reason, in
mechanistic models the parameters have meaning by themselves, so in this case
the data sets are too sparse rather than the model too large, see also [11].
We are confronted with large possible parameter regions, rather than a unique
optimal point. In such cases many points allow equally valid best fits to
experimental data [15, 24]. There is also measurement noise, as well as
intrinsic stochasticity of these systems to consider. Thus, it is important to
include the parameter uncertainty in model based predictions. It is also of
interest to assess where the uncertainty in the prediction stems from to guide
new experiments. It is therefore we include tools for both uncertainty
quantification and global sensitivity analysis (GSA) in this package [11]
(Figure 1).

The uncertainty of the parameter estimates are quantified directly during the
parameter estimation process through a Bayesian approach. The aim is to get a
representative sample of the parameter region (the posterior distribution)
which allow a good fit to data (Figure 1). This done through Markov chain
Monte Carlo (MCMC) sampling [13] an an Approximate Bayesian Computation (ABC)
[21] framework. The standard Bayesian framework is likelihood-based (see e.g.
[18]), but it is common to use Approximate Bayesian Computation (ABC) [21]
when the likelihood is intractable or complex, as for stochastic models. For
practical reasons we use ABC for all our models, both stochastic and
deterministic ODE models (were a likelihood exist), with the aim to also
implement likelihood based methods in the future. The use of ABC with
deterministic ordinary differential equation (ODE) models (see e.g. [27]) have
been discussed recently [1] with the recommendation to use likelihood-based
Bayesian methods for these, when possible, because of computational
efficiency. It has also been emphasized that the measurement noise has to be
accounted for if ABC methods are used with deterministic ODE models, either
directly when the model is simulated or through _probabilistic acceptance
criteria_ ([1, 26, 28]).

Bayesian methods provide a good framework for uncertainty analysis and MCMC is
a good method to use the Bayes theorem for this problem class in practice;
but, MCMC and ABC222both are sampling methods have costly disadvantages.
Sampling (using either method) means to solve the model many thousands of
times (forward), accepting and rejecting moves; this is computationally costly
and forces us to use fast and efficient ODE solvers. Sampling explores the
volume of the parameter space within the limits of acceptable fits and thus
scales exponentially with increasing numbers of parameters333this is worse
than with optimization - a big optimization problem differs from a big
sampling problem by many orders of magnitude. The samples have intrinsic
(integrated) auto-correlation; every MCMC sample has an _effictive sample
size_ , which is lower if the auto-correlation is higher, see [2]. Bad MCMC
configurations (transition kernel choice, etc.) have large auto-correlations.
Furthermore, MCMC methods typically need a good initial location to begin
sampling and have a convergence phase, where the target distribution is
approached. Many implementations have been developed specifically to solve
these difficulties generally: parallel tempering [10], Hamiltonian Monte Carlo
[4], Fisher information as metric tensors for movement [14], and automatic
tuning of acceptance rates [6]. However, for some model classes444including
ordinary differential equation models these solutions can further increase the
computational costs dramatically; the urgent need for numerically efficient
implementations remains, regardless of the chosen algorithm. Other approaches,
such as the profile likelihood method [19], address the scaling problem
directly by exploring the likelihood’s shape555in profile rather than sampling
from a target distribution. Our package aims to offer a use-friendly tool with
good sampling methods specifically for models in biology.

In previous work, we introduced the use of Vine-copulas to model complex
multivariate probability distributions between sampling runs, for use in model
calibration [11]. Vine copulas are flexible tools to model high-dimensional
probability distributions [3]. By separating the dependency structure from the
marginal distributions and using bivariate copulas with Vines, complex
distributions can be described. We recommend using Vine copulas especially in
the setting when an already calibrated model needs to be adapted to new
experimental data sets. As new data become available, the posterior
distribution of the parameters from the previous model calibration can be
approximated via Vine copulas and used as the prior for the fitting of the new
data sets, avoiding the need to re-calibrate the model with all available data
sets.

In [11] we also introduced the possibility perform a global sensitivity
analysis [25] on the posterior distribution from an uncertainty quantification
run (see Figure 1). This could, for example, be done in order to find out
which parameter that if it was known, would reduce the uncertainty of the
predictions the most; to guide further experiments. This required adaptions to
standard global sensitivity analysis algorithms were the input factors (e.g.
model parameters) are assumed to be orthogonal, in order to handle the non-
orthogonal inputs [25, 11] resulting from a uncertainty quantification run.

SBtab (table format for Systems Biology) is designed to support automated data
integration and model building [20]. The format can be adapted to new types of
data and can be used for data exchange in Systems Biology. SBtab relies on the
structure of spreadsheets and includes predefined table types for experimental
data and SBML-compliant model constituents. SBtab is well suited to deal with
model calibration projects as the model, the experimental data, and the prior
assumptions made on the possible parameter distributions can be defined within
SBtab. These prior assumptions (prior distribution) can for example correspond
to ranges for all parameters taken from previous knowledge, or to known
physical constraints.

The UQSA tool allows the user to chose whether to simulate the model as a
stochastic or deterministic system. Several other implementations of model
fitting via MCMC and/or ABC in R exist (e.g. [16], [9]), however, the UQSA
tool is more comprehensive, offering a full framework for UQSA on biochemical
network models.

Figure 1: Full model calibration requires the uncertainty of the parameter
estimates to be quantified through a Bayesian approach and that this
uncertainty (the posterior distribution) is propagated to the predictions made
from the model. Finally a global sensitivity analysis can be performed on the
posterior distribution to guide further experiments. Figure modified from [11]

##  2 Implementation

In our implementation, functionality that can be more broadly useful, has been
turned into a standalone package when feasible. We use these packages as
dependencies, if needed, especially in examples. But, they use standard R
variable types (no special classes), and can be replaced with alternative
methods to achieve the same effect. The packages we use in our examples can be
found in the same GitHub account; we consider them a part of the UQSA tool-
set.

###  2.1 Features of the UQSA tool-set

  * •

An easy to use, human, and machine readable format for reaction based models
and calibration data in the form of SBtab.

  * •

Functions for converting models and experimental data from SBtab spreadsheets
into R variables (icpm-kth/SBtabVFGEN). While parsing a model, these functions
also create R and C files666we use the solvers in the GNU Scientific Library
that contain functions with the ODE vector field and its Jacobian, as well as
other functions related to the model, [12, 17].

  * •

Functions777interfaces for simulating the model either as an ODE system in R
or C [12, 17] or as a stochastic model using the GillespieSSA2 package [7].

  * •

An MCMC algorithm to sample from the posterior distribution using an ABC
scheme. For ODE models measurement noise is modeled.

  * •

An ABC distance function in the form of weighted euclidean distance, with
weights given by (the reciprocals of) the experimental measurement errors (if
available) [23]. The ABC distance function can be also user-defined.

  * •

Interfacing a Vine copula based approach to sample from a non-trivial, non-
orthogonal prior, most relevant for sequential fitting of data sets to a
model.

  * •

Functionality for global sensitivity analysis on orthogonal and non-orthogonal
input factors.

  * •

An approach to model simulation that is oriented around _experiments_ ; one
experiment may be composed of several simulations.

  * •

An event system for scheduled events within an experiment, to model sudden
activation (or other changes) within an experimental protocol.

###  2.2 Modelling examples

We provide a range of different models and corresponding experimental data to
demonstrate the use of the UQSA software. These examples are named: AKAR4,
AKAP79, CaMKII after some important species of the models. The AKAR4 model is
the smallest example and available as a Jupyter notebook in the format of a
tutorial.

  * •

AKAR4 [8], a small model with two simulation backends

    * –

ODE model with stochastic noise, also provided as a Jupyter notebook.

    * –

stochastic model, available as an R-script, using a Gillespie solver.

  * •

AKAP79 [5, 8], larger ODE model example provided as an R-script. The
experimental data correspond to time-series.

  * •

CaMKII [22, 11], larger ODE model example provided as an R-script. The
experimental data consist of dose-response curves.

##  3 Conclusion and outlook

The UQSA tool addresses a need for better validation of biochemical network
models occurring for example within neuroscience or systems biology projects.
With this package the uncertainty of parameter estimates can be provided at
the same time as the parameter estimation is performed. This uncertainty can
also be propagated to the predictions made by the model and a global
sensitivity analysis can be performed to guide further experiments (Figure 1).
Other tools for uncertainty quantification exist (see e.g. pyABC [26]), the
advantages with this toolbox is that it is tailored for biochemical network
models (through the use of SBtab) and the ability to preform uncertainty
quantification _as well as_ global sensitivity analysis. The users get a
choice between different numerical solvers for the simulations, in R and C.
This means that you can run your model either as an ODE or a stochastic
system. It is also fast and should be relatively easy to use.

Apart from being necessary in order to falsify a specific model structure,
uncertainty quantification is also important for the reproducibility of the
model calibration process. As the aim is to provide a representative sample of
all statistically viable parameter sets, different calibration results can be
compared. In comparison, one (locally) optimal value can differ between
different optimization runs (depending on method, or starting locations); the
parameters may also be unidentifiable (even structurally) making
reproducibility more difficult in the case of optimization.

It is however important to note that when using UQSA only the parameter
uncertainty is quantified. The uncertainty of the model structure is not
investigated.

This toolbox is under further development and our aims are to provide a
likelihood based MCMC algorithm, parallel tempering, as well as improve user-
friendliness.

## Acknowledgements

The computations were enabled by resources provided by the National Academic
Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National
Infrastructure for Computing (SNIC) at PDC Centre for High Performance
Computing partially funded by the Swedish Research Council through grant
agreements no. 2022-06725 and no. 2018-05973.

## Funding

This work has been supported by The Swedish e-Science Research centre (SeRC)
and the Horizon 2020 Framework Programme (945539 (The Human Brain Project))

## References

  * [1] Amani A Alahmadi, Jennifer A Flegg, Davis G Cochrane, Christopher C Drovandi, and Jonathan M Keith.  A comparison of approximate versus exact techniques for bayesian parameter inference in nonlinear ordinary differential equation models.  Royal Society open science, 7(3):191315, 2020.
  * [2] Benjamin Ballnus, Sabine Hug, Kathrin Hatz, Linus Görlitz, Jan Hasenauer, and Fabian J Theis.  Comprehensive benchmarking of markov chain monte carlo methods for dynamical systems.  BMC Systems Biology, 11(1):1–18, 2017.
  * [3] Tim Bedford and Roger M Cooke.  Vines–a new graphical model for dependent random variables.  The Annals of Statistics, 30(4):1031–1068, 2002.
  * [4] Michael Betancourt, Simon Byrne, Sam Livingstone, and Mark Girolami.  The geometric foundations of Hamiltonian Monte Carlo.  Bernoulli, 23(4A):2257 – 2298, 2017.
  * [5] Joseph Daniel Buxbaum and Yadin Dudai.  A quantitative model for the kinetics of camp-dependent protein kinase (type ii) activity: Long-term activation of the kinase and its possible relevance to learning and memory.  Journal of Biological Chemistry, 264(16):9344–9351, 1989.
  * [6] Mylène Bédard.  Optimal acceptance rates for metropolis algorithms: Moving beyond 0.234.  Stochastic Processes and their Applications, 118(12):2198–2222, 2008.
  * [7] Robrecht Cannoodt, Wouter Saelens, Louise Deconinck, and Yvan Saeys.  Spearheading future omics analyses using dyngen, a multi-modal simulator of single cells.  Nature Communications, 2021.
  * [8] Timothy W Church, Parul Tewatia, Saad Hannan, João Antunes, Olivia Eriksson, Trevor G Smart, Jeanette Hellgren Kotaleski, and Matthew G Gold.  Akap79 enables calcineurin to directly suppress protein kinase a activity.  Elife, 10:e68164, 2021.
  * [9] Katalin Csilléry, Olivier François, and Michael GB Blum.  abc: an r package for approximate bayesian computation (abc).  Methods in ecology and evolution, 3(3):475–479, 2012.
  * [10] David J. Earl and Michael W. Deem.  Parallel tempering: Theory, applications, and new perspectives.  Phys. Chem. Chem. Phys., 7:3910–3916, 2005.
  * [11] Olivia Eriksson, Alexandra Jauhiainen, Sara Maad Sasane, Andrei Kramer, Anu G Nair, Carolina Sartorius, and Jeanette Hellgren Kotaleski.  Uncertainty quantification, propagation and characterization by bayesian analysis combined with global sensitivity analysis applied to dynamical intracellular pathway models.  Bioinformatics, 35(2):284–292, 2019.
  * [12] Mark Galassi, Jim Davies, James Theiler, Brian Gough, Gerard Jungman, Patrick Alken, Michael Booth, and Fabrice Rossi.  Gnu scientific library.  No. Release, 2, 1996.
  * [13] Andrew Gelman, John B Carlin, Hal S Stern, and Donald B Rubin.  Bayesian data analysis.  Chapman and Hall/CRC, 1995.
  * [14] Mark Girolami and Ben Calderhead.  Riemann Manifold Langevin and Hamiltonian Monte Carlo Methods.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 73(2):123–214, 03 2011.
  * [15] Ryan N Gutenkunst, Joshua J Waterfall, Fergal P Casey, Kevin S Brown, Christopher R Myers, and James P Sethna.  Universally sloppy parameter sensitivities in systems biology models.  PLoS computational biology, 3(10):e189, 2007.
  * [16] Franck Jabot, Thierry Faure, and Nicolas Dumoulin.  Easy abc: performing efficient approximate bayesian computation sampling schemes using r.  Methods in Ecology and Evolution, 4(7):684–687, 2013.
  * [17] Karline Soetaert, Thomas Petzoldt, and R. Woodrow Setzer.  Solving differential equations in R: Package deSolve.  Journal of Statistical Software, 33(9):1–25, 2010.
  * [18] Andrei Kramer, Jan Hasenauer, Frank Allgöwer, and Nicole Radde.  Computation of the posterior entropy in a bayesian framework for parameter estimation in biological networks.  In 2010 IEEE International Conference on Control Applications, pages 493–498. IEEE, 2010.
  * [19] Clemens Kreutz, Andreas Raue, Daniel Kaschek, and Jens Timmer.  Profile likelihood in systems biology.  The FEBS Journal, 280(11):2564–2571, 2013.
  * [20] Timo Lubitz, Jens Hahn, Frank T. Bergmann, Elad Noor, Edda Klipp, and Wolfram Liebermeister.  SBtab: a flexible table format for data exchange in systems biology.  Bioinformatics, 32(16):2559–2561, 2016.
  * [21] Paul Marjoram, John Molitor, Vincent Plagnol, and Simon Tavaré.  Markov chain monte carlo without likelihoods.  Proceedings of the National Academy of Sciences, 100(26):15324–15328, 2003.
  * [22] Anu G Nair, Omar Gutierrez-Arenas, Olivia Eriksson, Alexandra Jauhiainen, Kim T Blackwell, and Jeanette H Kotaleski.  Modeling intracellular signaling underlying striatal function in health and disease.  Progress in molecular biology and translational science, 123:277–304, 2014.
  * [23] Dennis Prangle.  Adapting the ABC Distance Function.  Bayesian Analysis, 12(1):289 – 309, 2017.
  * [24] Andreas Raue, Clemens Kreutz, Thomas Maiwald, Julie Bachmann, Marcel Schilling, Ursula Klingmüller, and Jens Timmer.  Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood.  Bioinformatics, 25(15):1923–1929, 2009.
  * [25] Andrea Saltelli, Marco Ratto, Terry Andres, Francesca Campolongo, Jessica Cariboni, Debora Gatelli, Michaela Saisana, and Stefano Tarantola.  Global sensitivity analysis: the primer.  John Wiley & Sons, 2008.
  * [26] Yannik Schälte, Emmanuel Klinger, Emad Alamoudi, and Jan Hasenauer.  pyabc: Efficient and robust easy-to-use approximate bayesian computation.  Journal of Open Source Software, 7(74):4304, 2022.
  * [27] Tina Toni, David Welch, Natalja Strelkowa, Andreas Ipsen, and Michael PH Stumpf.  Approximate bayesian computation scheme for parameter inference and model selection in dynamical systems.  Journal of the Royal Society Interface, 6(31):187–202, 2009.
  * [28] Richard David Wilkinson.  Approximate Bayesian computation (ABC) gives exact results under the assumption of model error.  Statistical Applications in Genetics and Molecular Biology, 12(2):129–141, 2013.