Motivation: Mechanistic kinetic models usually contain unknown parameters, which need to be estimated by optimizing the fit of the model to experimental data. This task can be computationally challenging due to the presence of local optima and ill-conditioning. While a variety of optimization methods have been suggested to surmount these issues, it is not obvious how to choose the best one for a given problem a priori, since many factors can influence their performance. A systematic comparison of methods that are suited to parameter estimation problems of sizes ranging from tens to hundreds of optimization variables is currently missing, and smaller studies indeed provided contradictory findings. Results: Here, we use a collection of benchmark problems to evaluate the performance of two families of optimization methods: (i) a multi-start of deterministic local searches; and (ii) a hybrid metaheuristic combining stochastic global search with deterministic local searches. A fair comparison is ensured through a collaborative evaluation, involving researchers applying each method on a daily basis, and a consideration of multiple performance metrics capturing the trade-off between computational efficiency and robustness. Our results show that, thanks to recent advances in the calculation of parametric sensitivities, a multi-start of gradient-based local methods is often a successful strategy, but a better performance can be obtained with a hybrid metaheuristic. The best performer is a combination of a global scatter search metaheuristic with an interior point local method, provided with gradients estimated with adjoint-based sensitivities. We provide an implementation of this novel method in an open-source software toolbox to render it available to the scientific community. Availability and Implementation: The code to reproduce the results is available at Zenodo https://doi.org/10.5281/

SummarySingle-cell time-lapse studies have advanced the quantitative understanding of cell-to-cell variability. However, as the information content of individual experiments is limited, methods to integrate data collected under different conditions are required.Here we present a multi-experiment nonlinear mixed effect modeling approach for mechanistic pathway models, which allows the integration of multiple single-cell perturbation experiments. We apply this approach to the translation of green fluorescent protein after transfection using a massively parallel read-out of micropatterned single-cell arrays. We demonstrate that the integration of data from perturbation experiments allows the robust reconstruction of cell-to-cell variability, i.e., parameter densities, while each individual experiment provides insufficient information. Indeed, we show that the integration of the datasets on the population level also improves the estimates for individual cells by breaking symmetries, although each of them is only measured in one experiment. Moreover, we confirmed that the suggested approach is robust with respect to batch effects across experimental replicates and can provide mechanistic insights into the nature of batch effects. We anticipate that the proposed multi-experiment nonlinear mixed effect modeling approach will serve as a basis for the analysis of cellular heterogeneity in single-cell dynamics.

Anti-cancer drugs commonly target signal transduction proteins activated by mutation. In patients with BRAF V600E melanoma, small molecule RAF and MEK kinase inhibitors cause dramatic but often transient tumor regression. Emerging evidence suggests that cancer cells adapting by non-genetic mechanisms constitute a reservoir for the development of drug-resistant tumors. Here, we show that few hours after exposure to RAF/MEK inhibitors, BRAF V600E melanomas undergo adaptive changes involving disruption of negative feedback and sporadic pulsatile reactivation of the MAPK pathway, so that MAPK activity is transiently high enough in some cells to drive proliferation. Quantitative proteomics and computational modeling show that pulsatile MAPK reactivation is possible due to the co-existence in cells of two MAPK cascades: one driven by BRAF V600E that is drug-sensitive and a second driven by receptors that is drug-resistant. Paradoxically, this may account both for the frequent emergence of drug resistance and for the tolerability of RAF/MEK therapy in patients.

MotivationParameter estimation methods for ordinary differential equation (ODE) models of biological processes can exploit gradients and Hessians of objective functions to achieve convergence and computational efficiency. However, the computational complexity of established methods to evaluate the Hessian scales linearly with the number of state variables and quadratically with the number of parameters. This limits their application to low-dimensional problems.ResultsWe introduce second order adjoint sensitivity analysis for the computation of Hessians and a hybrid optimization-integration based approach for profile likelihood computation. Second order adjoint sensitivity analysis scales linearly with the number of parameters and state variables. The Hessians are effectively exploited by the proposed profile likelihood computation approach. We evaluate our approaches on published biological models with real measurement data. Our study reveals an improved computational efficiency and robustness of optimization compared to established approaches, when using Hessians computed with adjoint sensitivity analysis. The hybrid computation method was more than two-fold faster than the best competitor. Thus, the proposed methods and implemented algorithms allow for the improvement of parameter estimation for medium and large scale ODE models.AvailabilityThe algorithms for second order adjoint sensitivity analysis are implemented in the Advance MATLAB Interface CVODES and IDAS (AMICI, https://github.com/ICB-DCM/AMICI/). The algorithm for hybrid profile likelihood computation is implemented in the parameter estimation toolbox (PESTO, https://github.com/ICB-DCM/PESTO/). Both toolboxes are freely available under the BSD license.Contactjan.hasenauer@helmholtz-muenchen.deSupplementary informationSupplementary data are available at Bioinformatics online.

Motivation: Because they effectively represent mass action kinetics, ordinary differential equation models are widely used to describe biochemical processes. Optimization-based calibration of these models on experimental data can be challenging, even for low-dimensional problems. However, reliable model calibration is a prerequisite for many subsequent analysis steps, including uncertainty analysis, model selection and biological interpretation. Although multiple hypothesis have been advanced to explain why optimization based calibration of biochemical models is challenging, there are few comprehensive studies that test these hypothesis and tools for performing such studies are also lacking. Results: We implemented an established trust-region method as a modular python framework (fides) to enable structured comparison of different approaches to ODE model calibration involving Hessian approximation schemes and trust-region subproblem solvers. We evaluate fides on a set of benchmark problems that include experimental data. We find a high variability in optimizer performance among different implementations of the same algorithm, with fides performing more reliably that other implementations investigated. Our investigation of possible sources of poor optimizer performance identify shortcomings in the widely used Gauss-Newton approximation. We address these shortcomings by proposing a novel hybrid Hessian approximation scheme that enhances optimizer performance.

Motivation: Mechanistic models of biochemical reaction networks facilitate the quantitative understanding of biological processes and the integration of heterogeneous datasets.However, some biological processes require the consideration of comprehensive reaction networks and therefore large-scale models. Parameter estimation for such models poses great challenges, in particular when the data are on a relative scale.Results: Here, we propose a novel hierarchical approach combining (i) the efficient analytic evaluation of optimal scaling, offset, and error model parameters with (ii) the scalable evaluation of objective function gradients using adjoint sensitivity analysis. We evaluate the properties of the methods by parameterizing a pan-cancer ordinary differential equation model (>1000 state variables, >4000 parameters) using relative protein, phospho-protein and viability measurements. The hierarchical formulation improves optimizer performance considerably. Furthermore, we show that this approach allows estimating error model parameters with negligible computational overhead when no experimental estimates are available, providing an unbiased way to weight heterogeneous data. Overall, our hierarchical formulation is applicable to a wide range of models, and allows for the efficient parameterization of largescale models based on heterogeneous relative measurements.Contact: jan.hasenauer@helmholtz-muenchen.de Supplementary information: Supplementary information are available at bioRxiv online.Supplementary code and data are available online at

Anti-cancer drugs commonly target signal transduction proteins activated by mutation. In patients with BRAFV600E melanoma, small molecule RAF and MEK kinase inhibitors cause dramatic but often transient tumor regression. Emerging evidence suggests that cancer cells adapting by non-genetic mechanisms constitute a reservoir for the development of drug-resistant tumors. Here, we show that few hours after exposure to RAF/MEK inhibitors, BRAFV600E melanomas undergo adaptive changes involving disruption of negative feedback and sporadic pulsatile reactivation of the MAPK pathway, so that MAPK activity is transiently high enough in some cells to drive proliferation. ERK pulses are locally triggered by autocrine/paracrine factors in the local microenvironment so that ERK activity levels are low on average, but transiently high enough in patches of cells to explain the infrequent division characteristic of persister cells. Quantitative proteomics and computational modeling show that pulsatile MAPK reactivation is possible due to the co-existence in cells of two MAPK cascades: one driven by BRAFV600E that is drug-sensitive and a second driven by receptors that is drug-resistant. Paradoxically, this may account both for the frequent emergence of drug resistance and for the tolerability of RAF/MEK therapy in patients. Citation Format: Luca Gerosa, Christopher Chidley, Fabian Froehlich, Gabriela Sanchez, Sang Kyun Lim, Jeremy Muhlich, Jia-Yun Chen, Gregory J. Baker, Denis Schapiro, Tujin Shi, Lian Yi, Carrie D. Nicora, Wei-Jun Qian, H. Steven Wiley, Peter Sorger. ERK pulses drive non-genetic resistance in drug-adapted BRAFV600E melanoma cells [abstract]. In: Proceedings of the AACR-NCI-EORTC International Conference on Molecular Targets and Cancer Therapeutics; 2019 Oct 26-30; Boston, MA. Philadelphia (PA): AACR; Mol Cancer Ther 2019;18(12 Suppl):Abstract nr LB-B09. doi:10.1158/1535-7163.TARG-19-LB-B09