Conformal Prediction in Dynamic Biological Systems
Alberto Portela Julio R. Banga Marcos Matabuena


  
  
Conformal Prediction in Dynamic Biological Systems
    Alberto Portela
 
Computational Biology Lab
 
MBG-CSIC, Spanish National Research Council
 
Pontevedra, Spain 
 
 
Julio R. Banga
 
Computational Biology Lab
 
MBG-CSIC, Spanish National Research Council
 
Pontevedra, Spain 
 
 
Marcos Matabuena
 
Department of Biostatistics
 
Harvard University
 
Cambridge, MA, USA 
 
 

    September 9, 2024
=======================================================================================================================================================================================================================================================================================================================================




§ ABSTRACT      
 Uncertainty quantification (UQ) is the process of systematically determining and characterizing the degree of confidence in computational model predictions. In the context of systems biology, especially with dynamic models, UQ is crucial because it addresses the challenges posed by nonlinearity and parameter sensitivity, allowing us to properly understand and extrapolate the behavior of complex biological systems. Here, we focus on dynamic models represented by deterministic nonlinear ordinary differential equations. Many current UQ approaches in this field rely on Bayesian statistical methods. While powerful, these methods often require strong prior specifications and make parametric assumptions that may not always hold in biological systems. Additionally, these methods face challenges in domains where sample sizes are limited, and statistical inference becomes constrained, with computational speed being a bottleneck in large models of biological systems. As an alternative, we propose the use of conformal inference methods, introducing two novel algorithms that, in some instances, offer non-asymptotic guarantees, enhancing robustness and scalability across various applications. We demonstrate the efficacy of our proposed algorithms through several scenarios, highlighting their advantages over traditional Bayesian approaches. The proposed methods show promising results for diverse biological data structures and scenarios, offering a general framework to quantify uncertainty for dynamic models of biological systems.The software for the methodology and the reproduction of the results is available at https://zenodo.org/doi/10.5281/zenodo.13644870https://zenodo.org/doi/10.5281/zenodo.13644870.

Contact: mailto:mmatabuena@hsph.harvard.edummatabuena@hsph.harvard.edu; mailto:j.r.banga@csic.esj.r.banga@csic.es







§ INTRODUCTION









In the field of systems biology, we utilize mechanistic dynamic models as tools to analyze, predict and understand the intricate behaviors of complex biological processes <cit.>. These models are typically composed of sets of deterministic nonlinear ordinary differential equations, and are designed to provide a quantitative understanding of the dynamics which would be difficult to achieve through other means. In particular, this mechanistic approach offers several advantages over data-driven approaches <cit.>. First, it can generate more accurate predictions and can be applied to a broader range of situations. Second, it provides a deeper understanding of how the system works thanks to its mechanism-based nature, making it easier to interpret the reasons behind its behavior. Finally, it requires less data for training because it is based on established theories and principles that describe the underlying processes. Overall, mechanistic models can help in understanding the dynamics of biological systems, in predicting their behaviour under different conditions, in generating testable hypotheses, and in identifying knowledge gaps.

However, these benefits come with a trade-off. As the number of elements (species) and unknown variables in the system increases, the model becomes significantly more complex in terms of number of parameters and non-linear relationships. This complexity can make it difficult to interpret the model's results <cit.>, and damages its identifiability, i.e. the ability to uniquely determine the unknown parameters in a model from the available data. As the complexity increases with more species and unknown parameters, achieving full identifiability and observability becomes more difficult <cit.>. As a result, developing a reliable dynamic mechanistic model can be a demanding and error prone task, requiring considerable expertise and the use of comprehensive and systematic protocols <cit.>. Additionally, highly detailed models can lead to more uncertain predictions <cit.>. The uncertainty in parameters and their identifiability impacts model predictions, thus ideally we should be able to characterize such impact in an interpretable manner, aiming to make useful predictions despite poor identifiability <cit.>. Therefore, quantifying this uncertainty and how it affects different system states, a process known as Uncertainty Quantification (UQ) <cit.>, is an open and fundamental challenge <cit.>.

UQ plays a key role in enhancing the reliability and interpretability of mechanistic dynamic models <cit.>. It helps in understanding the underlying uncertainties in the model parameters and predictions, thereby improving the model’s predictive power and its utility in decision-making processes <cit.>. A lack of proper UQ can result in models that are too confident in their predictions, which can be misleading. 

Different approaches for uncertainty quantification and robustness analysis in the context of systems biology have been reviewed elsewhere <cit.>. Roughly speaking, we can distinguish between Bayesian and frequentist methods <cit.>. 
Lately, the most predominant approach in the literature is to use Bayesian methods, which treat model parameters as random variables. Bayesian methods can perform well even with small sample sizes, especially when informative priors are used. Frequentist methods often require larger sample sizes to achieve reliable estimates. In practice, Bayesian approaches require parametric assumptions to define likelihood equations and the specification of a prior to derive an approximate posterior distribution of parameters for approximate and analytical methods <cit.>. This process is crucial for model estimation and inference.

However, Bayesian approaches often demand significant computational resources. Moreover, in the particular case of systems of differential equations, it is also typical to encounter identifiability issues, which result in multimodal posterior distributions that are challenging to handle in practice. In fact, although non-identifiabilities poses challenges to both frequentist and Bayesian sampling approaches, the latter can be especially susceptible to convergence failures <cit.>. Prediction profile likelihood methods and variants <cit.> provide a competitive alternative by combining a frequentist perspective with a maximum projection of the likelihood by solving a sequence of optimization problems. However, they can be computationally demanding when a large number of predictions must be assessed.

Although various uncertainty quantification approaches are utilized in the domain of systems biology, comparative assessments of the strengths and weaknesses of state-of-the-art methods remain scarce. <cit.> recently presented a systematic comparison of four methods: Fisher information matrix (FIM), Bayesian sampling, prediction profile likelihood and ensemble modelling. The comparison was made considering case studies of increasing computational complexity. This assessment revealed an interplay between their applicability and statistical interpretability. The FIM method was not reliable, and the prediction profile likelihood did not scale up well, being very computationally demanding when a large number of predictions had to be assessed. An interesting trade-off between computational scalability and accuracy and statistical guarantees was found for the ensemble and Bayesian sampling approaches. The Bayesian method proved adequate for less complex scenarios; however, it faced scalability challenges and encountered convergence difficulties when applied to the more intricate problems. The ensemble approach presented better performance for large-scale models, but weaker theoretical justification.

Therefore there is a clear need of UQ methods with both good scalability and strong theoretical statistical properties. Recently, in the literature of statistics and machine learning research, the use of conformal prediction <cit.> to quantify the uncertainty of model outputs has become increasingly popular as an alternative to Bayesian methods and other asymptotic approximations <cit.>. One of the successful aspects of this methodology in practice is the  non-asymptotic guarantees that ensure the coverage of prediction regions is well-calibrated, at least from a global (marginal) perspective, as reviewed by <cit.>. However, to the best of our knowledge, their use in the systems biology and dynamical systems literature is not widespread, despite their expected promising properties for examining and making predictions in complex biological systems.

Considering systems biology applications, we must account for the typically limited number of observations. In order to increase the statistical modelling efficiency, we focus on conformal predictions based on the estimation of the conditional mean regression function <cit.>, and the semi-parametric distributional location-scale regression models <cit.>. To exploit the specific structure of location-scale regression models and maximize the limited available information, we propose two algorithms based on the jackknife methodology (see for example <cit.>). Recently, conformal prediction has also been extended to accommodate general statistical objects, such as graphs and functions that evolve over time, which can be very relevant in many biological problems <cit.>.


The main contributions of our work are:






  * We propose and explore two conformal prediction algorithms for dynamical systems, specifically designed to optimize statistical efficiency when the measurement error of the biological system is homoscedastic. Specifically:

    

  * The first algorithm focuses on achieving a specific quantile of calibration for each dimension of the dynamical system, making it more flexible when the homoscedasticity assumption is violated in any given coordinate.


  * The second algorithm is specially designed for large dynamical systems. The core idea is to perform a global standardization of the residuals, deriving a global quantile of calibration to return the final prediction regions.



  * We illustrate both algorithms using several case studies of increasing complexity and evaluate their performance in terms of statistical efficiency and computation time respect to traditional uncertainity quantification Bayesian methods.











§ METHODOLOGY



 §.§ Modeling framework and notation

We consider dynamic models described by deterministic nonlinear ordinary differential equations (ODEs),

    ẋ(t)    = f(x,θ,t ,x(t_0)) 
    
            x(t_0)    = x_0
    
            y(t)    = g(x(t),θ,t),




in which x(t) ∈ℝ^n_x is the vector of state variables at time t, y(t) ∈ℝ^n_y is the vector of observables at time t, and θ∈ℝ^n_θ is the vector of unknown parameters. 

The vector field f: ℝ^n_x×ℝ^n_θ×ℝ↦ℝ^n_x and the mappings g: ℝ^n_x×ℝ^n_θ×ℝ↦ℝ^n_y and x_0: ℝ^n_θ↦ℝ^n_x are possibly nonlinear.

The calibration of ODE models requires the estimation of parameter vector θ from theoretical observations of the output y(t) at n times, denoted as t_1, t_2, …, t_n. The total number of measurements is n × n_y. In practice, we observe a perturbed multidimensional random variable 



ỹ= [   ỹ_1,1   ỹ_1,2       ⋯   ỹ_1,n;   ỹ_2,1   ỹ_2,2       ⋯   ỹ_2,n;       ⋮       ⋮       ⋱       ⋮; ỹ_n_y,1 ỹ_n_y,2       ⋯ ỹ_n_y,n ], and maybe we apply over ỹ a specific data transformation to convert the underlying model in homocedasticity in the transformated space. To be more precise, each random observation from the random vector ỹ, as defined through the probabilistic model introduced below:





    h_k(ỹ_k,i,λ) = h_k(y_k(t_i),λ) + ϵ_k(t_i)= h_k(g_k(x(t_i),θ),λ) + ϵ_k,i, 
      i = 1, …, n; k = 1, …, n_y, λ∈ℝ^n_y,



where, for each coordinate of dynamical system k = 1, …, n_y, ϵ_k ∈ℝ^n denotes the measurement random noise of mean equal to zero, and h_k(·, λ) is the specific transformation function that is an increasing real-valued function depending on a shape parameter λ. Trivial examples of such functions include the logarithmic function, log(·), which is commonly used in the dynamical systems literature due to its equivalence with the log-normal probabilistic model <cit.>.
For identifiability purposes, model (<ref>) assumes that the variance of the random error ỹ_k,i is solely a function of the regression function g_k(x(t_i), θ). This assumption encompasses specific heteroscedastic cases, such as the log-normal model mentioned earlier, where the signal-noise ratio is homoscedastic in the transformed space but not in the original space.



Model (<ref>) is a generalization of the Box-Cox transformation models, adapted for dynamical systems, and is defined in the regression literature by <cit.>. The authors refer to this semi-parametric transformation family as the transform-both-sides (TBS) model.



For simplicity, we assume that the measurement random noise is independent across the different temporal points t_i i=1,…, n,  follows a normal distribution,  and is homoscedastic across dimensions, i.e., ϵ_k,i = ϵ_k(t_i) ∼𝒩(0, σ_k^2), where k = 1, …, n_y. σ_k(t_i) denotes the standard deviation for k state of the dynamical systems in the temporal instant t_i. Here, we also assume that h_k(s,λ) = s for all s ∈ℝ, meaning that h_k is the identity function, or that we are  modeling the ODE systems directly from the original observations collected from the dynamical system. For simplicity, throughout this manuscript, we present all the modelling steps of the algorithms directly in terms of the regression function g(·, ·) to aid the reader. The transformed version of the different algorithms introduced here only involves applying the specified data transformation to the original sample y and then running the non-transformed algorithm in the transformed (or image) space.




The maximum likelihood estimate (MLE) of the vector of unknown parameters θ, denoted as θ, for any dataset 𝒟_n can be found by minimizing the log-likelihood function,


    θ= _θ∈ℝ^n_θℒ(θ;𝒟_n) = 1/2∑^n_y_k=1∑^n_i=1[ log(2πσ_k^2) + ( ỹ_k,i - g_k(x(t_i),θ)/σ_k)^2 ],




which is the predominant optimization approach in the field of dynamical biological systems. Another popular approach, when no external information about the probabilistic mechanism of the random error ϵ is available, is to use the minimization of the mean squared error as the optimization criterion to estimate the parameters of the regression function g.





 §.§ Conformal prediction for dynamical systems


Here, our objective is to present the development of new algorithms of conformal predictors for the class of dynamical systems stated above. The key idea is to consider the solution of the dynamical systems as the regression function g (for example, the conditional mean function) and the biological signals observed as corrupted by a measurement error ϵ. Then, using the residuals, it is possible to derive prediction regions using different conformal inference strategies. The underlying challenge in this type of regression is that the signal of the time series is observed in only a few n observations. Using full conformal methods, which require fitting hundreds or even thousands of models with n+1 observations, is not advisable due to the prohibitive computational cost. This cost arises from the need to obtain model parameters through the mathematical optimization of large and high-dimensional systems of differential equations.

On the other hand, the split conformal methods (we fit each model component using random, disjoint splits) are not a good strategy from a statistical efficiency perspective when typically n < 20 in this type of biological problems. To mitigate the limitations of full and split conformal, we propose to use two new specific methods of conformal prediction for dynamical systems that increase statistical efficiency in different scenarios and are based on the application of jackknife techniques <cit.>.

Given a new random observation Y_n+1= g(x(t_n+1),θ)+ϵ_n+1, independent and identically distributed (i.i.d) with respect to 𝒟_n (the original sample of n observations collected from the dynamical system), and for a predefined confidence level α∈ [0,1], the goal of predictive inference is to provide a prediction region C^α(·) ⊂ℝ^n_y such that ℙ(Y_n+1∈ C^α(t_n+1)| T=t_n+1) = 1 - α, where the probability event is conditioned to the temporal instant t_n+1. We assume for practical purposes that such  prediction region exists and is unique.

Conformal prediction <cit.>
 is a general uncertainty quantification framework that, independent of the underlying regression function g employed, provides non-asymptotic marginal (global) guarantees of the type ℙ(Y_n+1∈C^α(t_n+1)) ≥ 1 - α. In this case, the probability ℙ is over the random sample 𝒟_n ∪ (t_n+1,Y_n+1). In the literature, there are three variants of conformal predictions related to the partitioning of the original random sample 𝒟_n <cit.>: full, split, and jackknife conformal.

To determine whether a single prediction for the fixed pair (t_n+1, Y_n+1) is within the prediction region for any confidence level α∈ [0, 1], full conformal prediction requires using all n+1 observations and estimating the regression model g. For each fixed pair (t_n+1, Y_n+1), this involves optimizing the parameter vector θ and evaluating whether the point lies inside or outside the prediction region by numerically approximating the predictions across a large-scale grid of observations, both in the temporal domain and in the states space. Such a procedure is computationally intensive, often requiring the estimation of hundreds of models, making it impractical for applications involving dynamical systems, especially in high-dimensional settings.

As a alternative, split conformal is valid for obtaining prediction regions for any new data point (t_n+1,Y_n+1), and often divides the sample into ⌊ n/2 ⌋ observations to estimate the function g and the remaining n-⌊ n/2 ⌋ to obtain the quantile for calibration to derive the final predictive regions.
Finally, the jackknife approach is an intermediate method for making predictions for any data point (t_n+1,Y_n+1) and derive the prediction regions without sacrificing statistical efficiency. The jackknife approach requires the fitting of n predictive models, excluding in the i-th iteration the i^th observation. Here, in order to improve the robustness of the statistical properties of conformal jackknife, we derive two algorithms based on <cit.> and the , a new jackknife conformal strategy proposed in <cit.>.




 §.§ Algorithms

Algorithms <ref> and <ref> describe the core steps of our conformal UQ strategies for dynamical systems. In both algorithms, the first step involves excluding each i-th observation and fitting the regression functions to obtain the jackknife residuals. 

In the first algorithm, , we apply a version of conformal  in each coordinate of the dynamical system, introducing flexibility in the modeling in cases where uncertainty shape varies differently across each dimension. However, the theoretical convergence rates are often slower and require a larger number of observed data points in comparison with our second algorithm. 

To alleviate this issue and create a more efficient algorithm in some special homoscedastic cases, the second algorithm, , uses the hypothesis that the model is homoscedastic along each coordinate. In its second step, we standardize the residuals using a prior estimation of the standard deviation. Then, we consider the global quantile of calibration. Finally, by re-scaling with the specific standard deviation, we obtain the final prediction interval.






Further implementation details (as a Matlab package) are given in the Supplementary Information.





 §.§ Theory


For this theorem, and all results that follow, all probabilities are stated with respect to the distribution of the training data points (t_1, Y_1), …,(t_n, Y_n) and the test data point (t_n+1, Y_n+1) drawn i.i.d. from an arbitrary distribution P, and we assume implicitly that the regression method g is invariant to the ordering of the data–invariant to permutations. We will treat the sample size n ≥ 2 and the target coverage level α∈[0,1] as fixed throughout.

Theorem:    The conformal jacknife algorithms  and   satisfies
ℙ(Y_n+1∈C^α(t_n+1))
}≥ 1-2 α.

Consequence from <cit.>. The target in the inequality is 1-α that is reached often except in some non-trivial mathematical examples.






§ RESULTS

To assess the performance of the Uncertainty Quantification methods outlined in the previous section, we apply it to four case studies based on dynamic models of increasing complexity (Table 1 shows a summary): (i) a simple logistic growth model, (ii) a Lotka-Volterra classical model considering 2 species, (iii) the well-known α-pinene kinetic model, widely used in parameter estimation studies, and (iv) the challenging NFKB signalling pathway. For each one of these models we formulated several parameter estimation problems generating synthetic data sets considering different scenarios for data cardinality and noise levels. All problems are fully observed except the NFKB case study. Due to space constraints, details about the models and the different estimation problems are given in the Supplementary Information.







Below, we present and discuss these case studies where we have considered synthetic datasets that have been generated considering a noise model described by Equation (<ref>). For simplicity, we assumed the errors to be normally distributed, centered around the noise-free data sample, and we have adopted a homoscedastic model, i.e. the variance remains constant across each dimension of the dataset. More specifically,

    ỹ_k,i∼𝒩(y_k(t_i), σ_k^2),    i = 1, …, n;  k = 1, …, n_y,

where, σ_k=ϵ·μ_k, with μ_k=∑_i=1^ny_k(t_i)/n capturing the mean value of state k, and ϵ representing the percentage of added noise. The ODE parameters θ are estimated by MLE equations for gaussian data as is specified in the Equation (<ref>).  Our results are compared with a Bayesian method (STAN), as described in the Supplementary Information.



 §.§ Case I: Logistic growth model

To evaluate the performance of our methods on this case study, we considered various scenarios with different noise levels (0%, 1%, 5% and 10%) and dataset sizes (10, 20, 50 and 100 data points). For each combination of noise level and dataset size, we generated 50 different synthetic datasets, totaling 800 unique datasets. By generating multiple datasets for each scenario, we were able to obtain a robust estimate of the methods' behavior and assess their consistency across different realizations of the data.

The comparative analysis of the logistic growth model, as shown in Figure <ref>, highlights the robustness of the proposed methods  and  compared to conventional methodologies such as the Bayesian approach implemented with STAN. For a 10-point synthetic dataset with a 10 percent noise level, the predictive regions obtained by both conformal methods showed good coverage without requiring prior calibration of the models, unlike the Bayesian approach. Moreover, both  and  yield predictive regions comparable to those generated by the jackknife+ method; however, in this particular case, the  method shows superior performance.

In terms of computational efficiency, the conformal methods proved to be marginally faster than STAN, even for a problem of this small size, with differences on the order of a few seconds. This makes them more suitable for real-time applications.

To examine the marginal coverage ℙ(Y_n+1∈C^α(X_n+1)) for α=0.05, 0.1, 0.5 of our first algorithm , see Figure <ref>
 for different signal noises and sample sizes. The figure indicates the good empirical performance of our algorithm, achieving the desired nominal level in expectation as is expected for the non-assympotic guarantees of conformal prediction methods.











 §.§ Case II: Lotka-Volterra model

For this case study we generated datasets with the same noise levels (0%, 1%, 5% and 10%) as in the previous example and three different sizes (30, 60 and 120 points). Additionally, for each combination of noise level and dataset size, we generated 50 different synthetic datasets, resulting in a total of 600 unique datasets. 

Figure <ref>
 shows the results in a 30-point Lotka-Volterra
dataset, indicating that the predictive regions generated by
the conformal methods and STAN are similar in terms of
coverage. However, as in the previous case,  and 
offer the advantage of not requiring extensive hyperparameter
tuning, while also being more computationally efficient. In this
particular example, while the bayesian method obtains results
within a timeframe on the order of minutes, both conformal
methods achieve this in a significantly shorter span, on the order of seconds.






 §.§ Case III: Isomerization of α-Pinene


The dataset generation procedure for this case study mirrored that used for the Logistic model, employing the same noise levels and dataset sizes. Although we generated synthetic datasets to assess the method's behavior, we illustrated this behavior with a real dataset from <cit.>.

Figure <ref>
 shows the resulting regions of the isomerization of α-Pinene by applying the different algorithms to the 9-point real dataset. The results are once again consistent between both conformal algorithms and closely align with the regions obtained using STAN. In terms of computational cost, the conformal algorithms are notably more efficient, requiring less than a minute to compute the regions, whereas the Bayesian approach takes several minutes.





 §.§ Case IV: NFKB signaling pathway


The dataset generation process was similar to the one used in the previous case studies. The results of applying different algorithms to a 13-point synthetic dataset are illustrated in the Supplementary Information. As shown, both methods based on conformal inference yield results that are in close agreement with each other. However, in this case, the computational cost difference becomes particularly significant. While the  and  algorithms compute the regions within a few minutes, the Bayesian approach requieres several hours to achieve the same task.




§ DISCUSSION


In this study, we present two algorithms,  and , using conformal methods to perform uncertainty quantification. These methods allow the computation of prediction regions for nonlinear dynamic models of biological systems. We assume that the signal-to-noise ratio is homoscedastic in the measurements or any transformation of the original data. We successfully compared the performance of these new methods with Bayesian approaches using a set of problems of increasing complexity. The main conclusions from the numerical results are summarized below.

Our algorithms were significantly faster than the Bayesian method (STAN) for the case studies examined. Our methods, which do not require tuning of hyperparameters, performed well in agreement with the Bayesian approach for smaller case studies and larger datasets (more than 50 temporal points). However, for high-dimensional biological systems, as illustrated with the NFKB case study, our conformal methods exhibited better accuracy, while STAN encountered convergence issues.

Furthermore, our methods, although not based on specified regression models, achieved good marginal coverage due to their non-asymptotic properties. In contrast, Bayesian methods, exemplified by STAN, showed a more critical impact of poor calibration on marginal coverage, especially for small sample sizes, due to the lack of non-asymptotic properties.

Our study also revealed that obtaining good coverage properties with a Bayesian method requires careful tuning of the prior, which can be challenging even for well-known small problems and may be very difficult for new, larger problems arising in real applications. We encountered convergence issues with the MCMC strategy in STAN, likely due to the multimodal nature of posterior distributions and identifiability issues, consistent with previous reports <cit.>.

A primary limitation of our new methods is that the prediction regions might occasionally take negative values when observed states are very close to zero, which may be unrealistic from a mechanistic perspective. This issue is observed in both the STAN Bayesian implementation and the conformal methods presented here. One potential cause is the assumption of homoscedasticity (i.e., consistent signal-to-noise ratio across the entire domain). Moreover, the underlying model may not be correctly specified across all parts of the domain.

A straightforward solution to this issue is to apply a prior data transformation, such as the log-normal transformation, which allows for modeling heteroscedastic scenarios, or a more general family of possible transformations introduced in the model (<ref>). To enhance the usability of our proposed methods within the scientific community, we have made the code and data from this study available in a public repository. Moving forward, we plan to offer updates, including optimized code for high-performance computers, novel validations, and automatic transformation approaches for various error structures in the modeling.

Overall, our study presents a new framework for uncertainty quantification in dynamical models using conformal prediction methods. This framework provides an alternative to classical Bayesian methods. Notably, our new methods are computationally scalable, which is crucial for large biological models. From a mathematical statistics perspective, they offer non-asymptotic guarantees and avoid the technical difficulties of calibrating prior functions necessary in Bayesian statistics. For future work, we suggest exploring conformal quantile algorithms for massive dynamical biological systems <cit.>. These algorithms typically provide better conditional coverage than other conformal algorithms <cit.> and do not require the assumption of symmetric random errors. However, applying quantile conformal algorithms in practice may require collecting more temporal observations of dynamical systems, which might not always be feasible in real-world scenarios.



§ ACKNOWLEDGMENTS

JRB acknowledges support from grant PID2020-117271RB-C22 (BIODYNAMICS) funded by MCIN/AEI/10.13039/01100011033, and from CSIC grant PIE 202470E108 (LARGO). The authors also wish to thank Javier Enrique Aguilar Romero for his assistance with the use of Stan. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript





















§ APPENDIX - SUPPLEMENTARY MATERIAL





§ MATLAB IMPLEMENTATION OF THE ALGORITHMS


We implemented our  and  algorithms in Matlab. Parameter estimations were formulated as the minimization of a least squares cost function subject the dynamics (described by the model ODEs) and parameter bounds. These problems are non-convex and were solved using a global hybrid method, enhanced scatter search () due to its good performance an robustness <cit.>.  is available in Matlab as part of the  optimization toolbox <cit.>. Our code also has dependencies with the the Optimization Toolbox and the Parallel Computing Toolbox.  

The software for the methodology and the reproduction of the results is available at https://zenodo.org/doi/10.5281/zenodo.13644870https://zenodo.org/doi/10.5281/zenodo.13644870. All computations were carried out on a PC DELL Precision 7920 workstation with dual Intel Xeon Silver 4210R processors.



§ COMPARISON WITH A BAYESIAN METHOD


Bayesian methods are a classical approach for performing automated uncertainty quantifications by estimating the posterior distribution P(θ|𝒟_n), where θ represents the parameter of interest and 𝒟_n = {X_i}_i=1^n denotes the observed data. The key components in Bayesian analysis are the prior distribution P(θ), which encapsulates our initial beliefs about θ, and the likelihood function P(𝒟_n |θ), which represents the probability of observing the data 𝒟_n given the parameter θ.

In many practical scenarios, computing the posterior distribution analytically is challenging. Markov Chain provide Monte Carlo (MCMC) methods provides a general and powerful techniques used to estimate the posterior distribution by generating samples from it. Notable MCMC algorithms include Metropolis-Hastings and Gibbs sampling.

Nowadays, there are general software tools available for implementing Bayesian inference and MCMC methods. One such tool is STAN. To use STAN, one writes a model in the STAN modeling language, which involves defining the data, parameters, and model (i.e., prior and likelihood). STAN can be seamlessly integrated with  through the  package <cit.>, allowing users to perform Bayesian analyses within the R environment. The  package provides functions to compile STAN models, fit them to data, and extract samples for posterior analysis. Our implementations of the different case studies are also available in the Zenodo link above.




§ CASE STUDIES





 





 




 §.§ Case I: Logistic growth model


As our initial case study we considered the well-known logistic model <cit.>, governed by a single differential equation with two unknown parameters. This model is frequently used in population growth and epidemic spread modeling.


    ẋ = rx(1-x/K).


Here, r represents the growth rate, and K denotes the carrying capacity. The initial condition considered in the generation of the datasets was x(0)=10. Additionally, the values of the parameters used were r=0.1 and K=100. The initial condition is assumed to be known across all case studies considered. Since this logistic model has an analytical solution, it facilitates the comparison of our methods' performance with other established conformal methods for algebraic models, such as the jackknife+ <cit.>.

To evaluate the performance of our methods on this case study, we considered various scenarios with different noise levels (0%, 1%, 5% and 10%) and dataset sizes (10, 20, 50 and 100 data points). For each combination of noise level and dataset size, we generated 50 different synthetic datasets, totaling 800 unique datasets. By generating multiple datasets for each scenario, we were able to obtain a robust estimate of the methods' behavior and assess their consistency across different realizations of the data.

The comparative analysis of the logistic growth model, as shown in Figure <ref>, highlights the robustness of the proposed methods  and  compared to conventional methodologies such as the Bayesian approach implemented with STAN. For a 10-point synthetic dataset with a 10 percent noise level, the predictive regions obtained by both conformal methods showed good coverage without requiring prior calibration of the models, unlike the Bayesian approach. Moreover, both  and  yield predictive regions comparable to those generated by the jackknife+ method; however, in this particular case, the  method shows superior performance.

In terms of computational efficiency, the conformal methods proved to be marginally faster than STAN, even for a problem of this small size, with differences on the order of a few seconds. This makes them more suitable for real-time applications.

To examine the marginal coverage ℙ(Y_n+1∈C^α(X_n+1)) for α=0.05, 0.1, 0.5 of our first algorithm , see Figure <ref>
 for different signal noises and sample sizes. The figure indicates the good empirical performance of our algorithm, achieving the desired nominal level in expectation.






[!htbp]
Numerical results corresponding to the comparative analysis of the Logistic model predictive regions. This table presents the 95% lower and upper predictive bounds (LPB and UPB, respectively) obtained for each time point (t) in a 10-point dataset subjected to 10% noise. The observed data (y) and the true state value (x_nom) are also shown. The results are reported for four methodologies: the proposed  and  methods, the original jackknife+ approach, and a Bayesian method implemented with STAN.


            2cCUQDyn1     2cCUQDyn2     2cJackknife+     2cSTAN

4-56-78-910-11
t     y     x_nom     LPB     UPB     LPB     UPB     LPB     UPB     LPB     UPB

0     1.000e+01     1.000e+01     1.000e+01     1.000e+01     1.000e+01     1.000e+01     1.000e+01     1.000e+01     1.000e+01     1.000e+01     

10     2.226e+01     2.320e+01     7.698e+00     4.440e+01     7.213e+00     4.357e+01     9.493e+00     4.507e+01     5.545e+00     4.756e+01     

20     6.066e+01     4.509e+01     3.092e+01     7.100e+01     3.318e+01     6.954e+01     3.140e+01     7.043e+01     2.904e+01     7.781e+01     

30     7.327e+01     6.906e+01     5.593e+01     9.594e+01     5.916e+01     9.552e+01     5.578e+01     9.442e+01     5.324e+01     1.021e+02     

40     9.123e+01     8.585e+01     7.284e+01     1.095e+02     7.453e+01     1.109e+02     7.296e+01     1.083e+02     7.015e+01     1.138e+02     

50     8.895e+01     9.428e+01     7.816e+01     1.149e+02     8.116e+01     1.175e+02     7.888e+01     1.143e+02     7.681e+01     1.198e+02     

60     9.703e+01     9.782e+01     7.997e+01     1.167e+02     8.357e+01     1.199e+02     8.105e+01     1.169e+02     7.905e+01     1.215e+02     

70     1.004e+02     9.919e+01     8.055e+01     1.177e+02     8.441e+01     1.208e+02     8.181e+01     1.189e+02     7.917e+01     1.232e+02     

80     1.174e+02     9.970e+01     8.073e+01     1.183e+02     8.468e+01     1.210e+02     8.207e+01     1.198e+02     7.982e+01     1.229e+02     

90     1.135e+02     9.989e+01     8.079e+01     1.185e+02     8.477e+01     1.211e+02     8.215e+01     1.201e+02     7.953e+01     1.234e+02     

100     9.390e+01     9.996e+01     8.081e+01     1.185e+02     8.480e+01     1.212e+02     8.218e+01     1.203e+02     8.060e+01     1.230e+02     








 §.§ Case II: Lotka-Volterra model


As a second case study, we considered a two species Lotka-Volterra model <cit.>, often referred to as the predator-prey model. This model provides a fundamental framework for studying the dynamics between two interacting species. In its simplest form, it describes the interactions between a predator species and a prey species through a set of coupled differential equations with four unknown parameters:



    ẋ_1    = x_1(α-β x_2), 
    ẋ_2    = -x_2(γ-δ x_1).


Here, x_1 and x_2 represent the populations of the prey and predator, respectively. The parameters α, β, γ and δ are positive constants representing the interactions between the two species. Specifically, this parameters dictate the growth rates and interaction strengths, capturing the essence of biological interactions such as predation and competition. The initial conditions considered in the generation of the datasets were x(0)=(10,5). Additionally, the values of the parameters used were α=γ=0.5 and β=δ=0.02.

For this case study we generated datasets with the same noise levels (0%, 1%, 5% and 10%) as in the previous example and three different sizes (30, 60 and 120 points). Additionally, for each combination of noise level and dataset size, we generated 50 different synthetic datasets, resulting in a total of 600 unique datasets. 

Figure <ref>
 shows the results in a 30-point Lotka-Volterra dataset, indicating that the predictive regions generated by the conformal methods and STAN are similar in terms of coverage. However, as in the previous case,  and  offer the advantage of not requiring extensive hyperparameter tuning, while also being more computationally efficient. In this particular example, while the bayesian method obtains results within a timeframe on the order of minutes, both conformal methods achieve this in a significantly shorter span, on the order of seconds.










 §.§ Case III: Isomerization of α-Pinene


As a third case study, we examined the α-pinene isomerization model. The isomerization process of α-pinene is significant in industry, especially in the production of synthetic fragrances and flavors. These complex biochemical reactions can be effectively modeled using a system of five differential equations with five unknown parameters. The resulting kinetic model has been a classical example in the analysis of multiresponse data <cit.>. The kinetic equations encapsulate the transformation dynamics of α-pinene into its various isomers through a series of reaction steps: 


    ẋ_1   =-(p_1+p_2)x_1, 
    ẋ_2   =p_1x_1,
    ẋ_3   =p_2x_1-(p_3+p_4)x_3+p_5x_5,
    ẋ_4   =p_3x_3,
    ẋ_5   =p_4x_3-p_5x_5.


In the equations above, each p_i∈ [0,1], i=1,…,5 represents a different rate of reaction, defining the conversion speed from one isomer to another. The initial conditions considered in the generation of the datasets were x(0)=(100,0,0,0,0). Additionally, the values of the parameters used were p=(5.93e-05, 2.96e-05,2.05e-05,2.75e-04,4.00e-05). The dataset generation procedure for this case study mirrored that used for the Logistic model, employing the same noise levels and dataset sizes. Although we generated synthetic datasets to assess the method's behavior, we illustrated this behavior with a real dataset from <cit.>.

Figure <ref>
 shows the resulting regions of the isomerization of α-Pinene by applying the different algorithms to the 9-point real dataset. The results are once again consistent between both conformal algorithms and closely align with the regions obtained using STAN. In terms of computational cost, the conformal algorithms are notably more efficient, requiring less than a minute to compute the regions, whereas the Bayesian approach takes several minutes.



    
    < g r a p h i c s >





    
    < g r a p h i c s >





    
    < g r a p h i c s >





    
    < g r a p h i c s >



















[!htbp]
Parameter estimation comparison for the first three case studies. Since the parameter estimation process is identical for both algorithms proposed in this paper, for simplicity, we will refer to the parameters obtained with them as CUQDyn. For the case study on the isomerization of α-pinene, a real dataset was considered, and thus the nominal parameters for this problem are unknown.
=0pt


    2@c@Logistic     4@c@Lotka-Volterra     5@c@α-Pinene 

2-34-78-12
     r     K     α     β     γ     δ     p_1     p_2     p_3     p_4     p_5 

True      1.000e-01     1.000e+02     5.000e-01     2.000e-02     5.000e-01     2.000e-02     NA     NA     NA     NA     NA 

STAN      1.152e-01     1.021e+02     4.992e-01     2.004e-02     5.000e-01     2.001e-02     5.998e-05     2.797e-05     1.860e-05     2.826e-04     5.366e-05

CUQDyn      1.112e-01     1.027e+02     4.984e-01     2.001e-02     5.001e-01     2.000e-02     6.307e-05     2.841e-05     1.609e-05     2.729e-04     4.365e-05






 §.§ Case IV: NFKB signaling pathway


The Nuclear Factor Kappa-light-chain-enhancer of activated B cells (NFKB) signaling pathway plays a key role in the regulation of immune response, inflammation and cell survival. This pathway is activated in response to various stimuli, including cytokines, stress and microbial infections, leading yo the transcription of target genes involved in immune and inflammatory responses.
Here we consider the dynamics of this pathway as described by a system of differential equations <cit.>:


    İK̇K̇ṅ   = kprod - kdeg · IKKn
       - Tr · k1 · IKKn, 
    İK̇K̇ȧ   = Tr · k1 · IKKn - k3 · IKKa 
       - Tr · k2 · IKKa · A20 - kdeg · IKKa
       - a2 · IKKa · IkBa + t1 · IKKaIkBa
       - a3 · IKKa · IkBaNFkB
       + t2 · IKKaIkBaNFkB, 
    İK̇K̇i̇   = k3 · IKKa + Tr · k2 · IKKa · A20
       - kdeg · IKKi, 
    İK̇K̇ȧİk̇Ḃȧ   = a2 · IKKa · IkBa - t1 · IKKaIkBa, 
    İK̇K̇ȧİk̇ḂȧṄḞk̇Ḃ   = a3 · IKKa · IkBaNFkB
       - t2 · IKKaIkBaNFkB, 
    ṄḞk̇Ḃ   = c6a · IkBaNFkB - a1 · NFkB · IkBa
       + t2 · IKKaIkBaNFkB - i1 · NFkB, 
    ṄḞk̇Ḃṅ   = i1 · kv · NFkB - a1 · IkBan · NFkBn, 
    Ȧ2̇0̇   = c4 · A20t - c5 · A20, 
    Ȧ2̇0̇ṫ   = c2 + c1 · NFkBn - c3 · A20t, 
    İk̇Ḃȧ   = -a2 · IKKa · IkBa
       - a1 · IkBa · NFkB
       + c4a · IkBat - c5a · IkBa - i1a · IkBa
       + e1a · IkBan, 
    İk̇Ḃȧṅ   = -a1 · IkBan · NFkBn + i1a · kv · IkBa
       - e1a · kv · IkBan, 
    İk̇Ḃȧṫ   = c2a + c1a · NFkBn - c3a · IkBat, 
    İk̇ḂȧṄḞk̇Ḃ   = a1 · IkBa · NFkB - c6a · IkBaNFkB
        - a3 · IKKa · IkBaNFkB
       + e2a · IkBanNFkBn, 
    İk̇ḂȧṅṄḞk̇Ḃṅ   = a1 · IkBan · NFkBn
       - e2a · kv · IkBanNFkBn, 
    ċġėṅṫ   = c2c + c1c · NFkBn - c3c · cgent.


In agreement with the scenario considered by <cit.>, we assume that the available measurements are determined by the observation function g: ℝ^15→ℝ^6, which is defined as follows:


    g(·) =    (NFkBn(·), IkBa(·) + IkBaNFkB(·), A20t(·),
       IKKn(·) + IKKa(·) + IKKi(·), IKKa(·), IkBat(·)).


Out of the system of 15 equations, only 6 observables, defined by the function g, are available. The parameter values used in the generation of the datasets are as follows:


    a1   = 5e-01,    a2   = 2e-01,    t1   = 1e-01, 
    a3   = 1e+00,    t2   = 1e-01,   c1a   = 5e-07, 
    c2a   = 0e+00,    c3a   = 4e-04,    c4a   = 5e-01, 
    c5a   = 1e-04,    c6a   = 2e-05,    c1   = 5e-07,
    c2   = 0e+00,    c3   = 4e-04,    c4   = 5e-01, 
    c5   = 3e-04,    k1   = 2.5e-03,    k2   = 1e-01, 
    k3   = 1.5e-03,    kprod   = 2.5e-05,    kdeg   = 1.25e-04, 
    kv   = 5e+00,    i1   = 2.5e-03,    e2a   = 1e-02, 
    i1a   = 1e-03,    e1a   = 5e-04,    c1c   = 5e-07,
    c2c   = 0e+00,    c3c   = 4e-04.



It should be noted that in this problem, which involves 29 unknown parameters and 15 state variables, we only have access to 6 observable outputs. This discrepancy between the number of parameters and the available observables presents a challenge in the context of parameter identifiability, and is very commmon in systems biology applications. Identifiability refers to the ability to uniquely determine the model parameters based on the available data. When a system lacks identifiability, inferring unique parameter values from observable data becomes challenging, if not impossible. However, as mentioned in the introduction, by characterizing the impact of this lack of identifiability using appropriate uncertainty quantification (UQ) methods, it might still be possible to make useful predictions.




As shown, both methods based on conformal inference yield results that are in close agreement with each other. However, in this case, the computational benefits of our conformal strategies are especially pronounced. The  and  algorithms can compute the regions in just a few minutes, whereas the Bayesian approach requires several hours to accomplish the same task.




[Aldridge et al., 2006]albridge-burke-lauffenburger-sorger:2006
Aldridge, B., Burke, J., Lauffenburger, D., and Sorger, P. (2006).
Physicochemical modelling of cell signalling pathways.
Nature Cell Biology, 8(11):1195–1203.

[Angelopoulos et al., 2023]angelopoulos2023conformal
Angelopoulos, A. N., Bates, S., Lei, J., Fannjiang, C., Romano, Y., and Candès, E. (2023).
Conformal prediction: A gentle introduction.
Foundations and Trends® in Machine Learning, 16(4):494–591.

[Babtie and Stumpf, 2017]Babtie2017-rs
Babtie, A. C. and Stumpf, M. P. H. (2017).
How to deal with parameters for whole-cell modelling.
Journal of the Royal Society Interface, 14(133).

[Baker et al., 2018]Baker2018
Baker, R. E., Peña, J.-M., Jayamohan, J., and Jérusalem, A. (2018).
Mechanistic models versus machine learning, a fight worth fighting for the biological community?
Biology Letters, 14(5):20170660.

[Barber et al., 2021]10.1214/20-AOS1965
Barber, R. F., Candès, E. J., Ramdas, A., and Tibshirani, R. J. (2021).
Predictive inference with the jackknife+.
The Annals of Statistics, 49(1):486–507.

[Bayarri and Berger, 2004]Bayarri2004
Bayarri, M. J. and Berger, J. O. (2004).
The interplay of Bayesian and frequentist analysis.
Statistical Science, 19(1):58–80.

[Begoli et al., 2019]begoli2019need
Begoli, E., Bhattacharya, T., and Kusnezov, D. (2019).
The need for uncertainty quantification in machine-assisted medical decision making.
Nature Machine Intelligence, 1(1):20–23.

[Box et al., 1973]box1973some
Box, G., Hunter, W., MacGregor, J., and Erjavec, J. (1973).
Some problems associated with the analysis of multiresponse data.
Technometrics, 15(1):33–51.

[Carroll and Ruppert, 1984]carroll1984power
Carroll, R. J. and Ruppert, D. (1984).
Power transformations when fitting theoretical models to data.
Journal of the American Statistical Association, 79(386):321–328.

[Cedersund, 2012]cedersund:2012
Cedersund, G. (2012).
Conclusions via unique predictions obtained despite unidentifiability—new definitions and a general method.
FEBS Journal, 279:3513–3527.

[Cedersund, 2016]cedersund2016prediction
Cedersund, G. (2016).
Prediction uncertainty estimation despite unidentifiability: an overview of recent developments.
In Geris, L. and Gomez-Cabrero, D. (Eds.), Uncertainty in Biology, pages 449–466. Springer.

[Cedersund and Roll, 2009]cedersund-roll:2009
Cedersund, G. and Roll, J. (2009).
Model-based evaluation and comparison of potential explanations for given biological data.
FEBS Journal, 276(4):903–922.

[Coveney et al., 2016]Coveney2016
Coveney, P. V., Dougherty, E. R., and Highfield, R. R. (2016).
Big data need big theory too.
Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 374(2080):20160153.

[DiStefano III, 2015]distefano2015dynamic
DiStefano III, J. (2015).
Dynamic Systems Biology Modeling and Simulation.
Academic Press.

[Egea et al., 2014]egea2014meigo
Egea, J. A., Henriques, D., Cokelaer, T., Villaverde, A. F., MacNamara, A., Danciu, D.-P., Banga, J. R., and Saez-Rodriguez, J. (2014).
MEIGO: An open-source software suite based on metaheuristics for global optimization in systems biology and bioinformatics.
BMC Bioinformatics, 15(1):136.

[Eriksson et al., 2019]Eriksson2019
Eriksson, O., Jauhiainen, A., Maad Sasane, S., Kramer, A., Nair, A. G., Sartorius, C., and Hellgren Kotaleski, J. (2019).
Uncertainty quantification, propagation and characterization by Bayesian analysis combined with global sensitivity analysis applied to dynamical intracellular pathway models.
Bioinformatics, 35(2):284–292.

[Geris and Gomez-Cabrero, 2016]geris2016uncertainty
Geris, L. and Gomez-Cabrero, D., editors (2016).
Uncertainty in Biology: A Computational Modeling Approach.
Springer.

[Guo et al., 2020]guo2020package
Guo, J., Gabry, J., Goodrich, B., and Weber, S. (2020).
Package ‘rstan’.
R package version 2.21.2. Available at <https://cran.r-project.org/web/packages/rstan/>.

[Hass et al., 2016]hass2016fast
Hass, H., Kreutz, C., Timmer, J., and Kaschek, D. (2016).
Fast integration-based prediction bands for ordinary differential equation models.
Bioinformatics, 32(8):1204–1210.

[Hines et al., 2014]hines2014determination
Hines, K. E., Middendorf, T. R., and Aldrich, R. W. (2014).
Determination of parameter identifiability in nonlinear biophysical models: A Bayesian approach.
Journal of General Physiology, 143(3):401–416.

[Hinkley, 1979]Hinkley1979
Hinkley, D. (1979).
Predictive likelihood.
The Annals of Statistics, 7(4):718–728.

[Ingalls, 2013]ingalls2013mathematical
Ingalls, B. P. (2013).
Mathematical Modeling in Systems Biology: An Introduction.
MIT Press.

[Kaltenbach et al., 2009]kaltenbach2009systems
Kaltenbach, H.-M., Dimopoulos, S., and Stelling, J. (2009).
Systems analysis of cellular networks under uncertainty.
FEBS Letters, 583(24):3923–3930.

[Kirk et al., 2015]kirk2015systems
Kirk, P. D., Babtie, A. C., and Stumpf, M. P. (2015).
Systems biology (un)certainties.
Science, 350(6259):386–388.

[Kreutz et al., 2012]kreutz2012likelihood
Kreutz, C., Raue, A., and Timmer, J. (2012).
Likelihood-based observability analysis and confidence intervals for predictions of dynamic models.
BMC Systems Biology, 6(1):120.

[Lei et al., 2018]lei2018distribution
Lei, J., G'Sell, M., Rinaldo, A., Tibshirani, R. J., and Wasserman, L. (2018).
Distribution-free predictive inference for regression.
Journal of the American Statistical Association, 113(523):1094–1111.

[Liepe et al., 2014]liepe2014framework
Liepe, J., Kirk, P., Filippi, S., Toni, T., Barnes, C. P., and Stumpf, M. P. (2014).
A framework for parameter estimation and model selection from experimental data in systems biology using approximate Bayesian computation.
Nature Protocols, 9(2):439–456.

[Linden et al., 2022]Linden2022
Linden, N. J., Kramer, B., and Rangamani, P. (2022).
Bayesian parameter estimation for dynamical models in systems biology.
PLOS Computational Biology, 18(10):e1010651.

[Lipniacki et al., 2004]lipniacki-paszek-brasier-luxon-kimmel:2004
Lipniacki, T., Paszek, P., Brasier, A., Luxon, B., and Kimmel, M. (2004).
Mathematical model of NFκB regulatory module.
Journal of Theoretical Biology, 228:195–215.

[Lugosi and Matabuena, 2024]lugosi2024uncertainty
Lugosi, G. and Matabuena, M. (2024).
Uncertainty quantification in metric spaces.
arXiv preprint arXiv:2405.05110.

[Massonis et al., 2022]Massonis2022
Massonis, G., Villaverde, A. F., and Banga, J. R. (2022).
Improving dynamic predictions with ensembles of observable models.
Bioinformatics, 39(1).

[Matabuena et al., 2024]matabuena2024conformal
Matabuena, M., Ghosal, R., Mozharovskyi, P., Padilla, O. H. M., and Onnela, J.-P. (2024).
Conformal uncertainty quantification using kernel depth measures in separable Hilbert spaces.
arXiv preprint arXiv:2405.13970.

[Mitra and Hlavacek, 2019]Mitra2019
Mitra, E. D. and Hlavacek, W. S. (2019).
Parameter estimation and uncertainty quantification for systems biology models.
Current Opinion in Systems Biology, 18:9–18.

[Mikovičá and Hatzimanikatis, 2010]Mikovi2010
Mikovičá, L. and Hatzimanikatis, V. (2010).
Modeling of uncertainties in biochemical reactions.
Biotechnology and Bioengineering, 108(2):413–423.

[Murphy et al., 2024]Murphy2024
Murphy, R. J., Maclaren, O. J., and Simpson, M. J. (2024).
Implementing measurement error models with mechanistic mathematical models in a likelihood-based framework for estimation, identifiability analysis, and prediction in the life sciences.
Journal of the Royal Society Interface, 21(210).

[National Research Council et al., 2012]National_Research_Council2012-wq
National Research Council, Division on Engineering and Physical Sciences, Board on Mathematical Sciences and Their Applications, and Committee on Mathematical Foundations of Verification, Validation, and Uncertainty Quantification (2012).
Assessing the Reliability of Complex Models: Mathematical and Statistical Foundations of Verification, Validation, and Uncertainty Quantification.
National Academies Press.

[Plank and Simpson, 2024]plank2024structured
Plank, M. J. and Simpson, M. J. (2024).
Structured methods for parameter inference and uncertainty quantification for mechanistic models in the life sciences.
arXiv preprint arXiv:2403.01678.

[Prybutok et al., 2022]Prybutok2022-kz
Prybutok, A. N., Cain, J. Y., Leonard, J. N., and Bagheri, N. (2022).
Fighting fire with fire: Deploying complexity in computational modeling to effectively characterize complex biological systems.
Current Opinion in Biotechnology, 75:102704.

[Puy et al., 2022]Puy2022-yc
Puy, A., Beneventano, P., Levin, S. A., Lo Piano, S., Portaluri, T., and Saltelli, A. (2022).
Models with higher effective dimensions tend to produce more uncertain estimates.
Science Advances, 8(42):eabn9450.

[Quenouille, 1956]quenouille1956notes
Quenouille, M. H. (1956).
Notes on bias in estimation.
Biometrika, 43(3/4):353–360.

[Raue et al., 2013]raue2013joining
Raue, A., Kreutz, C., Theis, F. J., and Timmer, J. (2013).
Joining forces of Bayesian and frequentist methodology: A study for inference in the presence of non-identifiability.
Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(1984):20110544.

[Romano et al., 2019]romano2019conformalized
Romano, Y., Patterson, E., and Candès, E. (2019).
Conformalized quantile regression.
Advances in Neural Information Processing Systems, 32.

[Sesia and Candès, 2020]sesia2020comparison
Sesia, M. and Candès, E. J. (2020).
A comparison of some conformal quantile regression methods.
Stat, 9(1):e261.

[Shafer and Vovk, 2008]shafer2008tutorial
Shafer, G. and Vovk, V. (2008).
A tutorial on conformal prediction.
Journal of Machine Learning Research, 9(3):371–421.

[Sharp et al., 2022]sharp2022parameter
Sharp, J. A., Browning, A. P., Burrage, K., and Simpson, M. J. (2022).
Parameter estimation and uncertainty quantification using information geometry.
Journal of the Royal Society Interface, 19(189):20210940.

[Siegfried et al., 2023]siegfried2023distribution
Siegfried, S., Kook, L., and Hothorn, T. (2023).
Distribution-free location-scale regression.
The American Statistician, 77(4):345–356.

[Simpson and Maclaren, 2023]Simpson2023
Simpson, M. J. and Maclaren, O. J. (2023).
Profile-wise analysis: A profile likelihood-based workflow for identifiability analysis, estimation, and prediction with mechanistic mathematical models.
PLOS Computational Biology, 19(9):e1011515.

[Smith, 2013]Smith2013
Smith, R. C. (2013).
Uncertainty Quantification: Theory, Implementation, and Applications.
Society for Industrial and Applied Mathematics.

[Streif et al., 2016]streif2016robustness
Streif, S., Kim, K.-K. K., Rumschinski, P., Kishida, M., Shen, D. E., Findeisen, R., and Braatz, R. D. (2016).
Robustness analysis, prediction, and estimation for uncertain biochemical networks: An overview.
Journal of Process Control, 42:14–34.

[Tsoularis and Wallace, 2002]Tsoularis2002
Tsoularis, A. and Wallace, J. (2002).
Analysis of logistic growth models.
Mathematical Biosciences, 179(1):21–55.

[Vanlier et al., 2013]vanlier-tiemann-hilbers-vanriel:2013
Vanlier, J., Tiemann, C., Hilbers, P., and van Riel, N. (2013).
Parameter uncertainty in biochemical models described by ordinary differential equations.
Mathematical Biosciences, 246:305–314.

[Villaverde and Banga, 2014]villaverde2014reverse
Villaverde, A. F. and Banga, J. R. (2014).
Reverse engineering and identification in systems biology: Strategies, perspectives, and challenges.
Journal of the Royal Society Interface, 11(91):20130505.

[Villaverde et al., 2019]villaverde2018benchmarking
Villaverde, A. F., Froehlich, F., Weindl, D., Hasenauer, J., and Banga, J. R. (2019).
Benchmarking optimization methods for parameter estimation in large kinetic models.
Bioinformatics, 35(5):830–838.

[Villaverde et al., 2022]Villaverde2022protocol
Villaverde, A. F., Pathirana, D., Fröhlich, F., Hasenauer, J., and Banga, J. R. (2022).
A protocol for dynamic model calibration.
Briefings in Bioinformatics, 23(1):1–12.

[Villaverde et al., 2023]Villaverde2023_comparisonUQ
Villaverde, A. F., Raimúndez, E., Hasenauer, J., and Banga, J. R. (2023).
Assessment of prediction uncertainty quantification methods in systems biology.
IEEE/ACM Transactions on Computational Biology and Bioinformatics, 20(3):1725–1736.

[Wangersky, 1978]wangersky1978lotka
Wangersky, P. J. (1978).
Lotka-volterra population models.
Annual Review of Ecology and Systematics, 9:189–218.






