Scalable Parameter Estimation for Genome-Scale Biochemical Reaction Networks

Abstract

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.

Fig 1. Illustration of gradient calculation using finite differences, forward sensitivity analysis and adjoint sensitivity equations for a model of mRNA transfection.

(a) Sketch and mathematical formulation of the mathematical model of mRNA transfection presented by [36]. The intracellular release of mRNA at time point t_r is modeled using the Dirac delta distribution δ. (b) Illustration of finite differences, forward sensitivity analysis and adjoint sensitivity analysis for the model of mRNA transfection: (top) Step 1: simulation of model; (middle) Step 2: intermediate step for gradient calculation; and (bottom) Step 3: calculation of gradient from intermediate results. For all methods, Step 1 and 2 involve numerical simulation (the direction indicated by the arrow) and are computationally demanding, while Step 3 is computationally negligible.

Fig 2. Comparison of gradient computation times for finite differences and forward and adjoint sensitivity analysis.

(a) Scaling of computation time with respect to the number of parameters for the model of ErbB signaling (BM3). Computation time for finite differences and forward sensitivity equations increases roughly linearly. Computation time for adjoint sensitivity analysis is almost independent of the number of parameters but possesses a higher initial cost. Adjoint sensitivity analysis is 48 times faster than forward sensitivity analysis when considering all parameters. (b,c) Speedup when using adjoint sensitivity analysis over forward sensitivity analysis for gradient computation evaluated for all investigated models compared against n_θ and n_x ⋅ n_θ. Regression curves (dashed lines) have been fitted to the results of all models excluding B3, which seems to be an outlier. All computations were performed on a MacBook Pro with an 2.9 GHz Intel Core i7 processor.

Fig 3. Comparison of the gradients computed using adjoint sensitivity equations with gradients computed using finite differences and forward sensitivity equations with default accuracies (absolute error < 10⁻¹⁶, relative error < 10⁻⁸).

Each point represents the absolute value of one gradient element. Points on the diagonal indicate a good agreement. (a) Forward finite differences with ϵ = 10⁻³ vs. adjoint sensitivities. (b) Forward sensitivities vs. adjoint sensitivities. (c) Adjoint sensitivities with high accuracies (absolute error < 10⁻³², relative error < 10⁻¹⁶) and default accuracies (absolute error < 10⁻¹⁶, relative error < 10⁻⁸).

Fig 4. Comparison of optimization speed using forward and adjoint sensitivities for the model of ErbB signaling.

For local optimization using forward sensitivity analysis (trust-region method) and local optimization using adjoint sensitivity analysis (interior-point method) we quantified the computation time across 100 local optimization runs with different initial conditions. For 41 out of 100 initial points the gradient could not be evaluated due to numerical problems. These optimization runs are omitted in all further analysis. (a,c) Comparison of objective function value with respect to iteration number and computation time. The hulls and medians computed for both methods are depicted as shaded areas and solid lines. (b,d) Pairwise comparison of objective function value after 10 iterations and 5 hours for both methods. Each dot corresponds to one initial point for the optimization. The coloring indicates which method performed better. (e) Pairwise comparison of the time required to reach the final objective function value achieved in the forward approach. For the adjoint approach the equivalent time is the minimal time to reach the same objective function value. Each dot corresponds to one initial point for the optimization. (f) Histogram of speedup by using adjoint sensitivity analysis over forward sensitivity analysis for individual initial points, computed from (e). All computations were performed on a linux cluster. Runs with same initial conditions were carried out on the same computation node.