Tailored parameter optimization methods for ordinary differential equation models with steady-state constraints

Background Ordinary differential equation (ODE) models are widely used to describe (bio-)chemical and biological processes. To enhance the predictive power of these models, their unknown parameters are estimated from experimental data. These experimental data are mostly collected in perturbation experiments, in which the processes are pushed out of steady state by applying a stimulus. The information that the initial condition is a steady state of the unperturbed process provides valuable information, as it restricts the dynamics of the process and thereby the parameters. However, implementing steady-state constraints in the optimization often results in convergence problems. Results In this manuscript, we propose two new methods for solving optimization problems with steady-state constraints. The first method exploits ideas from optimization algorithms on manifolds and introduces a retraction operator, essentially reducing the dimension of the optimization problem. The second method is based on the continuous analogue of the optimization problem. This continuous analogue is an ODE whose equilibrium points are the optima of the constrained optimization problem. This equivalence enables the use of adaptive numerical methods for solving optimization problems with steady-state constraints. Both methods are tailored to the problem structure and exploit the local geometry of the steady-state manifold and its stability properties. A parameterization of the steady-state manifold is not required. The efficiency and reliability of the proposed methods is evaluated using one toy example and two applications. The first application example uses published data while the second uses a novel dataset for Raf/MEK/ERK signaling. The proposed methods demonstrated better convergence properties than state-of-the-art methods employed in systems and computational biology. Furthermore, the average computation time per converged start is significantly lower. In addition to the theoretical results, the analysis of the dataset for Raf/MEK/ERK signaling provides novel biological insights regarding the existence of feedback regulation. Conclusion Many optimization problems considered in systems and computational biology are subject to steady-state constraints. While most optimization methods have convergence problems if these steady-state constraints are highly nonlinear, the methods presented recover the convergence properties of optimizers which can exploit an analytical expression for the parameter-dependent steady state. This renders them an excellent alternative to methods which are currently employed in systems and computational biology. Electronic supplementary material The online version of this article (doi:10.1186/s12918-016-0319-7) contains supplementary material, which is available to authorized users.

(iii) The functionsF ,G, and x s and their partial derivatives up to order 2 are bounded for (x s −x s (θ)) ∈ B ζ (Assumption 1.4).
The properties (i)-(v) are the prerequisites of (Khalil, 2002, Theorem 9.3), establishing the existence of ε * > 0 such that for all ε < ε * systems of type (1) are locally exponentially stable. As stability properties are conserved under the performed transformations, we obtain for ε = λ −1 the Theorem 1.1.
Stability and convergence are not affected by the approximation of the steady-state sensitivity viaŜ. (2) These reactions describe binding of NGF to TrkA (R 1 and R 2 ), TrkA:NGF-mediated ERK phosphorylation (R 3 ), basal ERK phosphorylation (R 4 ) and ERK dephosphorylation (R 5 ). Brackets indicate the concentration of a biochemical species. The input is the initial NGF concentration, u = [NGF] 0 , and the measured output is the relative pERK concentration, The experimental noise is assumed to be normally distributed with the unknown variance σ 2 , ∼ N (0, σ 2 ).
In accordance with previous publications, we assume conservation of mass ( (Hasenauer et al., 2014) and references therein). Under these assumptions the activities of the NGF receptor TrkA and ERK are captured by in which x 1 = k 3 [TrkA:NGF] and x 2 = s[pERK]. This model possesses a minimal number of parameters, θ = (k 1 , k 2 , k 3 [TrkA] 0 , k 4 , s[ERK] 0 , k 5 , σ 2 ) and is structurally identifiable. For details on the model, we refer to the original publication (Hasenauer et al., 2014).

Derivation of the mathematical model for Raf/MEK/ERK signaling
We consider a model for Raf/MEK/ERK signaling which accounts for the core proteins and considers the six reactions: The upstream signaling is summarized in the time-dependent rate constant k 1,max (t) for which the flexible parameterization (as proposed in the Data2Dynamics toolbox (Raue et al., 2015)) is used. The effects of sorafenib and UO126 are captured by a reduction in the kinase activity of pRaf and pMEK (R 4 and R 6 ).
Experimental studies proved an inhibition of Raf phosphorylation by pERK (Fritsche-Guenther et al., 2011). This feedback is however context-dependent (Santos et al., 2007). To study the importance of this feedback during the G1/S phase transition, we considered two model hypotheses: with experimental input u = ([sora], [UO126]). The four collected Western blots, b = 1, . . . , 4, provide time-resolved relative data for different experimental inputs, with unknown, blot-dependent scaling constants s 1,b and s 2,b . The experimental noise for y 1,b and y 2,b is assumed to be normally distributed and proportional to the scaling constant, 1,b ∼ N (0, s 2 1,b σ 2 1 ) and 2 ∼ N (0, s 2 2,b σ 2 2 ). The unknown variances are denoted by σ 2 1 and σ 2 2 . for H1 and H2, model (7)-(8) possesses more than 20 parameters. Several of these parameters are structurally non-identifiable, including the absolute abundances of Raf, MEK and ERK. To circumvent these non-identifiablilities, we reformulate the model in terms of the fractions of phosphorylated proteins: ). The reformulated model does not depend explicitly on the total abundances [Raf] 0 , [MEK] 0 and [ERK] 0 but only on products and ratios of these parameters with other parameters, e.g., k 3 [Raf] 0 . Defining these products and ratios as new parameters eliminates non-identifiabilities and reduces the number of parameters. As all parameter are non-negative, a log-parameterization is used for parameter estimation (Raue et al., 2013). The states of the reformulated model are between 0 and 1.

Retraction factor
As an example, we studied the influence of the retraction factor λ on the convergence and computation time for Application example 1: NGF-induced ERK phosphorylation. To this end we considered retraction factors varying over several orders of magnitude and performed for each of these λ a multi-start optimization with 100 starts using simulation-based optimization with gradient descent and newton-type descent directions. We used the same setup as in Application example 1. The results are illustrated in Figure 1. In the case of gradient descent, the computation time decreases with increasing λ, which also results in decreasing average computation time per converged start. Choosing λ > 2000 cannot decrease the average computation time further. It can also be observed that there seems to be an upper limit on the maximal computation time. This is caused by a restriction on the maximal possible function evaluations allowed. In the case of a newtontype descent, increasing λ also decreases the computation time, however, not as drastically as when using gradient descent. Again, also the average computation time per converged start decreases, but plateaus at much smaller values of λ than in the gradient descent case.

Behavior in the presence of bistability, bifurcations and oscilations
Biological systems often possess multiple stable steady states and non-trivial ω-limit sets, e.g. stable limit cycles. To evaluate the simulation-based optimization method and the hybrid optimization method, we consider a bistable system (with hysteresis) and a system with a Hopf bifurcation.

Bistable system
To study the performance of the proposed methods in the presence of bistability we considered the systeṁ as described in (Müller & Kuttler, 2015, Chapter 5, p. 532). We assume direct observation of the state and the objective function J(x s , θ) = (x s − x s ) 2 withx s ≈ 1.5. The measurementx s is the larger of the two steady states for θ = 0.5. The resulting minimization problem is The optimization was performed using the hybrid and the simulation-based optimization method.
This is the normal form for systems with Hopf bifurcation as described in (Müller & Kuttler, 2015, Chapter 2, p. 227). The measurement was taken to be x s = (0, 0) T which corresponds to the steady state for θ ≤ 0. For θ > 0 this system exhibits stable limit cycles and no stable steady state. We considered a least squares objective function yielding the optimization problem min xs,θ x 2 s,i s.t. 0 = x s,2 + x s,1 (θ − x 2 s,1 − x 2 s,2 ) 0 = x s,1 + x s,2 (θ − x 2 s,1 − x 2 s,2 ).
The optimization was performed using only the simulation-based optimization method as the simulation in the hybrid optimization step is not converging for θ > 0 as no stable steady state exists. Figure 3 illustrates example optimizer trajectories. The results illustrate a good convergence of the simulation-based optimization method to the set of parameters θ with steady state x s = (0, 0) T , independent of the starting point.
In summary, the analysis of these simple systems indicates that in particular the simulation-based optimization method also facilitates the analysis of models with multiple stable steady states and limit cycles.