Scalable Parameter Estimation for Genome-Scale Biochemical Reaction Networks
In the life sciences , the abundance of experimental data is rapidly increasing due to the advent of novel measurement devices ., Genome and transcriptome sequencing , proteomics and metabolomics provide large datasets 1 at a steadily decreasing cost ., While these genome-scale datasets allow for a variety of novel insights 2 , 3 , a mechanistic understanding on the genome scale is limited by the scalability of currently available computational methods ., For small- and medium-scale biochemical reaction networks mechanistic modeling contributed greatly to the comprehension of biological systems 4 ., Ordinary differential equation ( ODE ) models are nowadays widely used and a variety of software tools are available for model development , simulation and statistical inference 5–7 ., Despite great advances during the last decade , mechanistic modeling of biological systems using ODEs is still limited to processes with a few dozens biochemical species and a few hundred parameters ., For larger models rigorous parameter inference is intractable ., Hence , new algorithms are required for massive and complex genomic datasets and the corresponding genome-scale models ., Mechanistic modeling of a genome-scale biochemical reaction network requires the formulation of a mathematical model and the inference of its parameters , e . g . reaction rates , from experimental data ., The construction of genome-scale models is mostly based on prior knowledge collected in databases such as KEGG 8 , REACTOME 9 and STRING 10 ., Based on these databases a series of semi-automatic methods have been developed for the assembly of the reaction graph 11–13 and the derivation of rate laws 14 , 15 ., As model construction is challenging and as the information available in databases is limited , in general , a collection of candidate models can be constructed to compensate flaws in individual models 16 ., For all these model candidates the parameters have to be estimated from experimental data , a challenging and usually ill-posed problem 17 ., To determine maximum likelihood ( ML ) and maximum a posteriori ( MAP ) estimates for model parameters , high-dimensional nonlinear and non-convex optimization problems have to be solved ., The non-convexity of the optimization problem poses challenges , such as local minima , which have to be addressed by the selection of optimization methods ., Commonly used global optimization methods are multi-start local optimization 18 , evolutionary and genetic algorithms 19 , particle swarm optimizers 20 , simulated annealing 21 and hybrid optimizers 22 , 23 ( see 18 , 24–26 for a comprehensive survey ) ., For ODE models with a few hundred parameters and state variables multi-start local optimization methods 18 and related hybrid methods 27 have proven to be successful ., These optimization methods use the gradient of the objective function to establish fast local convergence ., While the convergence of gradient based optimizers can be significantly improved by providing exact gradients ( see e . g . 18 , 28 , 29 ) , the gradient calculation is often the computationally most demanding step ., The gradient of the objective function is usually approximated by finite differences ., As this method is neither numerically robust nor computationally efficient , several parameter estimation toolboxes employ forward sensitivity analysis ., This decreases the numerical error and computation time 18 ., However , the dimension of the forward sensitivity equations increases linearly with both the number of state variables and parameters , rendering its application for genome-scale models problematic ., In other research fields such as mathematics and engineering , adjoint sensitivity analysis is used for parameter estimation in ordinary and partial differential equation models ., Adjoint sensitivity analysis is known to be superior to the forward sensitivity analysis when the number of parameters is large 30 ., Adjoint sensitivity analysis has been used for inference of biochemical reaction networks 31–33 ., However , the methods were never picked up by the systems and computational biology community , supposedly due to the theoretical complexity of adjoint methods , a missing evaluation on a set of benchmark models , and an absence of an easy-to-use toolbox ., In this manuscript , we provide an intuitive description of adjoint sensitivity analysis for parameter estimation in genome-scale biochemical reaction networks ., We describe the end value problem for the adjoint state in the case of discrete-time measurement and provide an user-friendly implementation to compute it numerically ., The method is evaluated on seven medium- to large-scale models ., By using adjoint sensitivity analysis , the computation time for calculating the objective function gradient becomes effectively independent of the number of parameters with respect to which the gradient is evaluated ., Furthermore , for large-scale models adjoint sensitivity analysis can be multiple orders of magnitude faster than other gradient calculation methods used in systems biology ., The reduction of the time for gradient evaluation is reflected in the computation time of the optimization ., This renders parameter estimation for large-scale models feasible on standard computers , as we illustrate for a comprehensive kinetic model of ErbB signaling ., We consider ODE models for biochemical reaction networks ,, x ˙ = f ( x , θ ) , x ( t 0 ) = x 0 ( θ ) , ( 1 ), in which x ( t , θ ) ∈ R n x is the concentration vector at time t and θ ∈ R n θ denotes the parameter vector ., Parameters are usually kinetic constants , such as binding affinities as well as synthesis , degradation and dimerization rates ., The vector field f : R n x × R n θ ↦ R n x describes the temporal evolution of the concentration of the biochemical species ., The mapping x 0 : R n θ ↦ R n x provides the parameter dependent initial condition at time t0 ., As available experimental techniques usually do not provide measurements of the concentration of all biochemical species , we consider the output map h : R n x × R n θ ↦ R n y ., This map models the measurement process , i . e . the dependence of the output ( or observables ) y ( t , θ ) ∈ R n y at time point t on the state variables and the parameters ,, y ( t , θ ) = h ( x ( t , θ ) , θ ) ., ( 2 ), The i-th observable yi can be the concentration of a particular biochemical species ( e . g . yi = xl ) as well as a function of several concentrations and parameters ( e . g . yi = θm ( xl1 + xl2 ) ) ., We consider discrete-time , noise corrupted measurements, y ¯ i j = y i ( t j , θ ) + ϵ i j , ϵ i j ∼ N ( 0 , σ i j 2 ) , ( 3 ), yielding the experimental data D = { ( ( y ¯ i j ) i = 1 n y , t j ) } j = 1 N . The number of time points at which measurements have been collected is denoted by N . Remark: For simplicity of notation we assume throughout the manuscript that the noise variances , σ i j 2 , are known and that there are no missing values ., However , the methods we will present in the following as well as the respective implementations also work when this is not the case ., For details we refer to the S1 Supporting Information ., We estimate the unknown parameter θ from the experimental data D using ML estimation ., Parameters are estimated by minimizing the negative log-likelihood , an objective function indicating the difference between experiment and simulation ., In the case of independent , normally distributed measurement noise with known variances the objective function is given by, J ( θ ) = 1 2 ∑ i = 1 n y ∑ j = 1 N y ¯ i j - y i ( t j , θ ) σ i j 2 , ( 4 ), where yi ( tj , θ ) is the value of the output computed from Eqs ( 1 ) and ( 2 ) for parameter value θ ., The minimization ,, θ * = arg min θ ∈ Θ J ( θ ) , ( 5 ), of this weighted least squares J yields the ML estimate of the parameters ., The optimization problem Eq ( 5 ) is in general nonlinear and non-convex ., Thus , the objective function can possess multiple local minima and global optimization strategies need to be used ., For ODE models multi-start local optimization has been shown to perform well 18 ., In multi-start local optimization , independent local optimization runs are initialized at randomly sampled initial points in parameter space ., The individual local optimizations are run until the stopping criteria are met and the results are collected ., The collected results are visualized by sorting them according to the final objective function value ., This visualization reveals local optima and the size of their basin of attraction ., For details we refer to the survey by Raue et al . 18 ., In this study , initial points are generated using latin hypercube sampling and local optimization is performed using the interior point and the trust-region-reflective algorithm implemented in the MATLAB function fmincon . m ., Gradients are computed using finite differences , forward sensitivity analysis or adjoint sensitivity analysis ., A näive approximation to the gradient of the objective function with respect to θk is obtained by finite differences ,, ∂ J ∂ θ k ≈ J ( θ + a e k ) - J ( θ - b e k ) a + b , ( 6 ), with a , b ≥ 0 and the kth unit vector ek ., In practice forward differences ( a = ϵ , b = 0 ) , backward differences ( a = 0 , b = ϵ ) and central differences ( a = ϵ , b = ϵ ) are widely used ., For the computation of forward finite differences , this yields a procedure with three steps: In theory , forward and backward differences provide approximations of order ϵ while central differences provide more accurate approximations of order ϵ2 , provided that J is sufficiently smooth ., In practice the optimal choice of a and b depends on the accuracy of the numerical integration 18 ., If the integration accuracy is high , an accurate approximation of the gradient can be achieved using a , b ≪ 1 . For lower integration accuracies , larger values of a and b usually yield better approximations ., A good choice of a and b is typically not clear a priori ( cf . 34 and the references therein ) ., The computational complexity of evaluating gradients using finite differences is affine linear in the number of parameters ., Forward and backward differences require in total nθ + 1 function evaluations ., Central differences require in total 2nθ function evaluations ., As already a single simulation of a large-scale model is time-consuming , the gradient calculation using finite differences can be limiting ., State-of-the-art systems biology toolboxes , such as the MATLAB toolbox Data2Dynamics 7 , use forward sensitivity analysis for gradient evaluation ., The gradient of the objective function is, ∂ J ∂ θ k = ∑ i = 1 n y ∑ j = 1 N y ¯ i j - y i ( t j , θ ) σ i j 2 s i , k y ( t j ) , ( 7 ), with s i , k y ( t ) :  t 0 , t N  ↦ R denoting the sensitivity of output yi at time point t with respect to parameter θk ., Governing equations for the sensitivities are obtained by differentiating Eqs ( 1 ) and ( 2 ) with respect to θk and reordering the derivatives ., This yields, s ˙ k x = ∂ f ∂ x s k x + ∂ f ∂ θ k , s k x ( t 0 ) = ∂ x 0 ∂ θ k s i , k y = ∂ h i ∂ x s k x + ∂ h i ∂ θ k ( 8 ), with s k x ( t ) :  t 0 , t N  ↦ R n x denoting the sensitivity of the state x with respect to θk ., Note that here and in the following , the dependencies of f , h , x0 and their ( partial ) derivatives on t , x and θ are not stated explicitly but have the to be assumed ., For a more detailed presentation we refer to the S1 Supporting Information Section 1 . Forward sensitivity analysis consists of three steps: Step 1 and 2 are often combined , which enables simultaneous error control and the reuse of the Jacobian 30 ., The simultaneous error control allows for the calculation of accurate and reliable gradients ., The reuse of the Jacobian improves the computational efficiency ., The number of state and output sensitivities increases linearly with the number of parameters ., While this is unproblematic for small- and medium-sized models , solving forward sensitivity equations for systems with several thousand state variable bears technical challenges ., Code compilation can take multiple hours and require more memory than what is available on standard machines ., Furthermore , while forward sensitivity analysis is usually faster than finite differences , in practice the complexity still increases roughly linearly with the number of parameters ., In the numerics community , adjoint sensitivity analysis is frequently used to compute the gradients of a functional with respect to the parameters if the function depends on the solution of a differential equation 35 ., In contrast to forward sensitivity analysis , adjoint sensitivity analysis does not rely on the state sensitivities s k x ( t ) but on the adjoint state p ( t ) ., The calculation of the objective function gradient using adjoint sensitivity analysis consists of three steps: Step 1 and 2 , which are usually the computationally intensive steps , are independent of the parameter dimension ., The complexity of Step 3 increases linearly with the number of parameters , yet the computation time required for this step is typically negligible ., The calculation of state and output trajectories ( Step 1 ) is standard and does not require special methods ., The non-trivial element in adjoint sensitivity analysis is the calculation of the adjoint state p ( t ) ∈ R n x ( Step 2 ) ., For discrete-time measurements—the usual case in systems and computational biology—the adjoint state is piece-wise continuous in time and defined by a sequence of backward differential equations ., For t > tN , the adjoint state is zero , p ( t ) = 0 . Starting from this end value the trajectory of the adjoint state is calculated backwards in time , from the last measurement t = tN to the initial time t = t0 ., At the time points at which measurements have been collected , tN , … , t1 , the adjoint state is reinitialised as, p ( t j ) = lim t → t j + p ( t ) + ∑ i = 1 n y ∂ h i ∂ x T y ¯ i j - y i ( t j ) σ i j 2 , ( 9 ), which usually results in a discontinuity of p ( t ) at tj ., Starting from the end value p ( tj ) as defined in Eq ( 9 ) the adjoint state evolves backwards in time until the next measurement point tj−1 or the initial time t0 is reached ., This evolution is governed by the time-dependent linear ODE, p ˙ = - ∂ f ∂ x T p ., ( 10 ), The repeated evaluation of Eqs ( 9 ) and ( 10 ) until t = t0 yields the trajectory of the adjoint state ., Given this trajectory , the gradient of the objective function with respect to the individual parameters is, ∂ J ∂ θ k = - ∫ t 0 t N p T ∂ f ∂ θ k d t - ∑ i , j ∂ h i ∂ θ k y ¯ i j - y i ( t j ) σ i j 2 - p ( t 0 ) T ∂ x 0 ∂ θ k ., ( 11 ), Accordingly , the availability of the adjoint state simplifies the calculation of the objective function to nθ one-dimensional integration problems over short time intervals whose union is the total time interval t0 , tN ., Algorithm 1: Gradient evaluation using adjoint sensitivity analysis % State and output Step 1 Compute state and output trajectories using Eqs ( 1 ) and ( 2 ) ., % Adjoint state Step 2 . 1 Set end value for adjoint state , ∀t > tN: p ( t ) = 0 . for j = N to 1 do Step 2 . 2 Compute end value for adjoint state according to the jth measurement using Eq ( 9 ) ., Step 2 . 3 Compute trajectory of adjoint state on time interval t = ( tj−1 , tj by solving Eq ( 10 ) ., end % Objective function gradient for k = 1 to nθ do Step 3 Evaluation of the sensitivity ∂J/∂θk using Eq ( 11 ) ., end Pseudo-code for the calculation of the adjoint state and the objective function gradient is provided in Algorithm 1 . We note that in order to use standard ODE solvers the end value problem Eq ( 10 ) can be transformed in an initial value problem by applying the time transformation τ = tN − t ., The derivation of the adjoint sensitivities for discrete-time measurements is provided in the S1 Supporting Information Section 1 . The key difference of the adjoint compared to the forward sensitivity analysis is that the derivatives of the state and the output trajectory with respect to the parameters are not explicitly calculated ., Instead , the sensitivity of the objective function is directly computed ., This results in practice in a computation time of the gradient which is almost independent of the number of parameters ., A visual summary of the different sensitivity analysis methods is provided in Fig 1 . Besides the procedures also the computational complexity is indicated ., The implementation of adjoint sensitivity analysis is non-trivial and error-prone ., To render this method available to the systems and computational biology community , we implemented the Advanced Matlab Interface for CVODES and IDAS ( AMICI ) ., This toolbox allows for a simple symbolic definition of ODE models ( 1 ) and ( 2 ) as well as the automatic generation of native C code for efficient numerical simulation ., The compiled binaries can be executed from MATLAB for the numerical evaluation of the model and the objective function gradient ., Internally , the SUNDIALS solvers suite is employed 30 , which offers a broad spectrum of state-of-the-art numerical integration of differential equations ., In addition to the standard functionality of SUNDIALS , our implementation allows for parameter and state dependent discontinuities ., The toolbox and a detailed documentation can be downloaded from http://ICB-DCM . github . io/AMICI/ ., For the comparison of different gradient calculation methods , we consider a set of standard models from the Biomodels Database 37 and the BioPreDyn benchmark suite 27 ., From the biomodels database we considered models for the regulation of insulin signaling by oxidative stress ( BM1 ) 38 , the sea urchin endomesoderm network ( BM2 ) 39 , and the ErbB sigaling pathway ( BM3 ) 40 ., From BioPreDyn benchmark suite we considered models for central carbon metabolism in E . coli ( B2 ) 41 , enzymatic and transcriptional regulation of carbon metabolism in E . coli ( B3 ) 42 , metabolism of CHO cells ( B4 ) 43 , and signaling downstream of EGF and TNF ( B5 ) 44 ., Genome-wide kinetic metabolic models of S . cerevisiae and E . coli ( B1 ) 45 contained in the BioPreDyn benchmark suite and the Biomodels Database 15 , 45 were disregarded due to previously reported numerical problems 27 , 45 ., The considered models possess 18-500 state variable and 86-1801 parameters ., A comprehensive summary regarding the investigated models is provided in Table 1 ., To obtain realistic simulation times for adjoint sensitivities realistic experimental data is necessary ( see S1 Supporting Information Section 3 ) ., For the BioPreDyn models we used the data provided in the suite , for the ErbB signaling pathway we used the experimental data provided in the original publication and for the remaining models we generated synthetic data using the nominal parameter provided in the SBML definition ., In the following , we will compare the performance of forward and adjoint sensitivities for these models ., As the model of ErbB signaling has the largest number of state variables and is of high practical interest in the context of cancer research , we will analyze the scalability of finite differences and forward and adjoint sensitivity analysis for this model in greater detail ., Moreover , we will compare the computational efficiency of forward and adjoint sensitivity analysis for parameter estimation for the model of ErbB signaling ., The evaluation of the objective function gradient is the computationally demanding step in deterministic local optimization ., For this reason , we compared the computation time for finite differences , forward sensitivity analysis and adjoint sensitivity analysis and studied the scalability of these approaches at the nominal parameter θ0 which was provided in the SBML definitions of the investigated models ., For the comprehensive model of ErbB signaling we found that the computation times for finite differences and forward sensitivity analysis behave similarly ( Fig 2a ) ., As predicted by the theory , for both methods the computation time increased linearly with the number of parameters ., Still , forward sensitivities are computationally more efficient than finite differences , as reported in previous studies 18 ., Adjoint sensitivity analysis requires the solution to the adjoint problem , independent of the number of parameters ., For the considered model , solving the adjoint problem a single time takes roughly 2-3-times longer than solving the forward problem ., Accordingly , adjoint sensitivity analysis with respect to a small number of parameter is disadvantageous ., However , adjoint sensitivity analysis scales better than forward sensitivity analysis and finite differences ., Indeed , the computation time for adjoint sensitivity analysis is almost independent of the number of parameters ., While computing the sensitivity with respect to a single parameter takes on average 10 . 09 seconds , computing the sensitivity with respect to all 219 parameters takes merely 14 . 32 seconds ., We observe an average increase of 1 . 9 ⋅ 10−2 seconds per additional parameter for adjoint sensitivity analysis which is significantly lower than the expected 3 . 24 seconds for forward sensitivity analysis and 4 . 72 seconds for finite differences ., If the sensitivities with respect to more than 4 parameters are required , adjoint sensitivity analysis outperforms both forward sensitivity analysis and finite differences ., For 219 parameters , adjoint sensitivity analysis is 48-times faster than forward sensitivities and 72-times faster than finite differences ., To ensure that the observed speedup is not unique to the model of ErbB signaling ( BM3 ) we also evaluated the speedup of adjoint sensitivity analysis over forward sensitivity analysis on models B2-5 and BM1-2 ., The results are presented in Fig 2b and 2c ., We find that for all models , but model B3 , gradient calculation using adjoint sensitivity is computationally more efficient than gradient calculation using forward sensitivities ( speedup > 1 ) ., For model B3 the backwards integration required a much higher number of integration steps ( 4 ⋅ 106 ) than the forward integration ( 6 ⋅ 103 ) , which results to a poor performance of the adjoint method ., One reason for this poor performance could be that , in contrast to other models , the right hand side of the differential equation of model B3 consists almost exclusively of non-linear , non-mass-action terms ., Excluding model B3 we find an polynomial increase in the speedup with respect to the number of parameters nθ ( Fig 2b ) , as predicted by theory ., Moreover , we find that the product nθ ⋅ nx , which corresponds to the size of the system of forward sensitivity equations , is an even better predictor ( R2 = 0 . 99 ) than nθ alone ( R2 = 0 . 83 ) ., This suggest that adjoint sensitivity analysis is not only beneficial for systems with a large number of parameters , but can also be beneficial for systems with a large number of state variables ., As we are not aware of any similar observations in the mathematics or engineering community , this could be due to the structure of biological reaction networks ., Our results suggest that adjoint sensitivity analysis is an excellent candidate for parameter estimation in large-scale models as it provides good scaling with respect to both , the number of parameters and the number of state variables ., Efficient local optimization requires accurate and robust gradient evaluation 18 ., To assess the accuracy of the gradient computed using adjoint sensitivity analysis , we compared this gradient to the gradients computed via finite differences and forward sensitivity analysis ., Fig 3 visualizes the results for the model of ErbB signaling ( BM3 ) at the nominal parameter θ0 which was provided in the SBML definition ., The results are similar for other starting points ., The comparison of the gradients obtained using finite differences and adjoint sensitivity analysis revealed small discrepancies ( Fig 3a ) ., The median relative difference ( as defined in S1 Supporting Information Section, 2 ) between finite differences and adjoint sensitivity analysis is 1 . 5 ⋅ 10−3 ., For parameters θk to which the objective function J was relatively insensitive , ∂J/∂θk < 10−2 , there are much higher discrepancies , up to a relative error of 2 . 9 ⋅ 103 ., Forward and adjoint sensitivity analysis yielded almost identical gradient elements over several orders of magnitude ( Fig 3b ) ., This was expected as both forward and adjoint sensitivity analysis exploit error-controlled numerical integration for the sensitivities ., To assess numerical robustness of adjoint sensitivity analysis , we also compared the results obtained for high and low integration accuracies ( Fig 3c ) ., For both comparisons we found the similar median relative and maximum relative error , namely 2 . 6 ⋅ 10−6 and 9 . 3 ⋅ 10−4 ., This underlines the robustness of the sensitivitity based methods and ensures that differences observed in Fig 3a indeed originate from the inaccuracy of finite differences ., Our results demonstrate that adjoint sensitivity analysis provides objective function gradients which are as accurate and robust as those obtained using forward sensitivity analysis ., As adjoint sensitivity analysis provides accurate gradients for a significantly reduced computational cost , this can boost the performance of a variety of optimization methods ., Yet , in contrast to forward sensitivity analysis , adjoint sensitivities do not yield sensitivities of observables and it is thus not possible to approximate the Hessian of the objective function via the Fisher Information Matrix 46 ., This prohibits the use of possibly more efficient Newton-type algorithms which exploit second order information ., Therefore , adjoint sensitivities are limited to quasi-Newton type optimization algorithms , e . g . the Broyden-Fletcher-Goldfarb-Shanno ( BFGS ) algorithm 47 , 48 , for which the Hessian is iteratively approximated from the gradient during optimization ., In principle , the exact calculation of the Hessian and Hessian-Vector products is possible via second order forward and adjoint sensitivity analysis 49 , 50 , which possess similar scaling properties as the first order methods ., However , both forward and adjoint approaches come at an additional cost and are thus not considered in this study ., To assess whether the use of adjoint sensitivities for optimization is still viable , we compared the performance of the interior point algorithm using adjoint sensitivity analysis with the BFGS approximation of the Hessian to the performance of the trust-region reflective algorithm using forward sensitivity analysis with Fisher Information Matrix as approximation of the Hessian ., For both algorithms we used the MATLAB implementation in fmincon . m ., The employed setup of the trust-region algorithm is equivalent to the use of lsqnonlin . m which is the default optimization algorithm in the MATLAB toolbox Data2Dynamics 7 , which was employed to win several DREAM challenges ., For the considered model the computation time of forward sensitivities is comparable in Data2Dynamics and AMICI ., Therefore , we expect that Data2Dynamics would perform similar to the trust-region reflective algorithm coupled to forward sensitivity analysis ., We evaluated the performance for the model of ErbB signaling based on 100 multi-starts which were initialized at the same initial points for both optimization methods ., For 41 out of 100 initial points the gradient could not be evaluated due numerical problems ., These optimization runs are omitted in all further analysis ., To limit the expected computation to a bearable amount we allowed a maximum of 10 iterations for the forward sensitivity approach and 500 iterations for the adjoint sensitivity approach ., As the previously observed speedup in gradient computation was roughly 48 fold , we expected this setup should yield similar computation times for both approaches ., We found that for the considered number of iterations , both approaches perform similar in terms of objective function value compared across iterations ( Fig 4a and 4b ) ., However , the computational cost of one iteration was much cheaper for the optimizer using adjoint sensitivity analysis ., Accordingly , given a fixed computation time the interior-point method using adjoint sensitivities outperforms the trust-region method employing forward sensitivities and the FIM ( Fig 4c and 4d ) ., In the allowed computation time , the interior point algorithm using adjoint sensitivities could reduce the objective function by up to two orders of magnitude ( Fig 4c ) ., This was possible although many model parameters seem to be non-identifiable ( see S1 Supporting Information Section 4 ) , which can cause problems ., To quantify the speedup of the optimization using adjoint sensitivity analysis over the optimization using forward sensitivity analysis , we performed a pairwise comparison of the minimal time required by the adjoint sensitivity approach to reach the final objective function value of the forward sensitivity approach for the individual points ( Fig 4e ) ., The median speedup achieved across all multi-starts was 54 ( Fig 4f ) , which was similar to the 48 fold speedup achieved in the gradient computation ., The availability of the Fisher Information Matrix for forward sensitivities did not compensate for the significantly reduced computation time achieved using adjoint sensitivity analysis ., This could be due to the fact that adjoint sensitivity based approach , being able to carry out many iterations in a short time-frame , can build a reasonable approximation of the Hessian approximation relatively fast ., In summary , this application demonstrates the applicability of adjoint sensitivity analysis for parameter estimation in large-scale biochemical reaction networks ., Possessing similar accuracy as forward sensitivities , the scalability is improved which results in an increased optimizer efficiency ., For the model of ErbB signaling , optimization using adjoint sensitivity analysis outperformed optimization using forward sensitivity analysis ., Mechanistic mathematical modeling at the genome scale is an important step towards a holistic understanding of biological processes ., To enable modeling at this scale , scalable computational methods are required which are applicable to networks with thousands of compounds ., In this manuscript , we present a gradient computation method which meets this requirement and which renders parameter estimation for large-scale models significantly more efficient ., Adjoint sensitivity analysis , which is extensively used in other research fields , is a powerful tool for estimating parameters of large-scale ODE models of biochemical reaction networks ., Our study of several benchmark models with up to 500 state variables and up to 1801 parameters demonstrated that adjoint sensitivity analysis provides accurate gradients in a computation time which is much lower than for established methods and effectively independent of the number of parameters ., To achieve this , the adjoint state is computed using a piece-wise continuous backward differential equation ., This backward differential equation has the same dimension as the original model , yet the computation time required to solve it usually is slightly larger ., As a result , finite differences and forward sensitivity analysis might be more efficient if the sensitivities with respect to a few parameters are required ., The same holds for alternatives like complex-step derivative approximation techniques 51 and forward-mode automatic differentiation 28 , 52 ., For systems with many parameters , adjoint sensitivity analysis is advantageous ., A scalable alternative might be reverse-mode automatic differentiation 28 , 53 , which remains to be evaluated for the considered class of problems .,
Mechanistic mathematical modeling of biochemical reaction networks using ordinary differential equation ( ODE ) models has improved our understanding of small- and medium-scale biological processes ., While the same should in principle hold for large- and genome-scale processes , the computational methods for the analysis of ODE models which describe hundreds or thousands of biochemical species and reactions are missing so far ., While individual simulations are feasible , the inference of the model parameters from experimental data is computationally too intensive ., In this manuscript , we evaluate adjoint sensitivity analysis for parameter estimation in large scale biochemical reaction networks ., We present the approach for time-discrete measurement and compare it to state-of-the-art methods used in systems and computational biology ., Our comparison reveals a significantly improved computational efficiency and a superior scalability of adjoint sensitivity analysis ., The computational complexity is effectively independent of the number of parameters , enabling the analysis of large- and genome-scale models ., Our study of a comprehensive kinetic model of ErbB signaling shows that parameter estimation using adjoint sensitivity analysis requires a fraction of the computation time of established methods ., The proposed method will facilitate mechanistic modeling of genome-scale cellular processes , as required in the age of omics .
In this manuscript , we introduce a scalable method for parameter estimation for genome-scale biochemical reaction networks ., Mechanistic models for genome-scale biochemical reaction networks describe the behavior of thousands of chemical species using thousands of parameters ., Standard methods for parameter estimation are usually computationally intractable at these scales ., Adjoint sensitivity based approaches have been suggested to have superior scalability but any rigorous evaluation is lacking ., We implement a toolbox for adjoint sensitivity analysis for biochemical reaction network which also supports the import of SBML models ., We show by means of a set of benchmark models that adjoint sensitivity based approaches unequivocally outperform standard approaches for large-scale models and that the achieved speedup increases with respect to both the number of parameters and the number of chemical species in the model ., This demonstrates the applicability of adjoint sensitivity based approaches to parameter estimation for genome-scale mechanistic model ., The MATLAB toolbox implementing the developed methods is available from http://ICB-DCM . github . io/AMICI/ .
applied mathematics, simulation and modeling, algorithms, optimization, genomic databases, mathematics, genome analysis, research and analysis methods, genome complexity, biological databases, differential equations, biochemistry, biochemical simulations, database and informatics methods, genetics, biology and life sciences, physical sciences, genomics, computational biology