 Methodology Article
 Open Access
Tailored parameter optimization methods for ordinary differential equation models with steadystate constraints
 Anna Fiedler^{1, 2},
 Sebastian Raeth^{3},
 Fabian J. Theis^{1, 2},
 Angelika Hausser^{3} and
 Jan Hasenauer^{1, 2}Email author
https://doi.org/10.1186/s1291801603197
© The Author(s) 2016
Received: 22 March 2016
Accepted: 12 July 2016
Published: 22 August 2016
Abstract
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 steadystate constraints in the optimization often results in convergence problems.
Results
In this manuscript, we propose two new methods for solving optimization problems with steadystate 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 steadystate constraints. Both methods are tailored to the problem structure and exploit the local geometry of the steadystate manifold and its stability properties. A parameterization of the steadystate 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 stateoftheart 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 steadystate constraints. While most optimization methods have convergence problems if these steadystate constraints are highly nonlinear, the methods presented recover the convergence properties of optimizers which can exploit an analytical expression for the parameterdependent steady state. This renders them an excellent alternative to methods which are currently employed in systems and computational biology.
Keywords
Background
Gene regulation, signal transduction, metabolism and many other biological processes are nowadays analyzed using mathematical models [1]. Mathematical models allow for the integration of available knowledge, providing mechanistic insights and an understanding of design principles [2]. The spectrum of employed modeling approaches ranges from qualitative Boolean models [3, 4] to quantitative stochastic spatial models [5]. In particular ODE models, such as reaction rate equations, are used on a range of spatial and temporal scales [6]. These models describe the temporal evolution of the concentration of biochemical species in cellular compartments as a function of initial conditions, parameters and stimuli.
ODE models are flexible and allow for the mechanistic description of a range of processes (see, e.g., [7–9] and references therein). Similar to other quantitative modeling approaches, ODE models rely on accurate values for initial conditions and parameters, e.g., binding affinities, synthesis and degradation rates. Initial conditions and parameters are often unknown and have to be inferred from experimental data [10].
The estimation of the parameters of ODE models from data collected during perturbation experiments requires the solution of optimization problems. These optimization problems are in general nonlinear, nonconvex and computationally demanding. This establishes the need for efficient and robust optimization methods [13, 14]. In the literature, multistart local optimization methods [15] and global optimization methods [16, 17] have been employed. If the stationarity of the initial condition does not have to be enforced or if analytical expressions of the steady state are available, these methods mostly work well [15, 17]. However, neglecting steadystate constraints often results in a loss of information and potentially implausible predictions. Furthermore, an analytical expression of the parameter and inputdependent steady state is rarely available.
Steadystate constraints are nonlinear equality constraints, which restrict the solution space to a manifold of feasible points. Enforcing these equality constraints causes efficiency and convergence problems for standard optimization methods [18]. Deterministic local optimization methods have to move along the manifold, resulting in small stepsizes or stagnation. Stochastic local and global optimization methods are only allowed to propose update steps on the manifold, which is not possible in stateoftheart toolboxes (see, e.g., [19, 20]). The elimination of the equality constraints using analytical expressions for the parameterdependent steady states resolves these problems [18].
The first method to derive analytical expressions for steady states has been proposed by [21] for networks of enzymecatalyzed reactions. This method was later implemented [22] and subsequently extended using results from graph theory [23, 24]. Furthermore, tailored methods for models with bilinear rate laws have been introduced [25]. For models with MichaelisMenten and HillType kinetics, pysubstitution has been developed [26]. This method solves the steadystate constraint for a combination of states and parameters. An important recent extension ensured positivity of the solution [18]. Pysubstitution and its extension are however only applicable if the model possesses sufficiently many degrees of freedom [18], which is difficult to assess a priori.
We propose two novel methods for solving optimization problems with a single or multiple steadystate constraints. These methods do not rely on an analytical expression for the steady state. Instead, the geometry of the steadystate manifold and the stability of the steady state are exploited.
The first method we introduce borrows ideas from optimization algorithms on matrix manifolds [27, 28]. These optimization algorithms employ retraction operators which map a point onto the manifold [28]. These retraction operators – usually analytical functions – facilitate an effective movement of optimizers on the manifold. Retraction operators are however problemspecific and their construction is usually nontrivial [28], which limits the application of established algorithms for optimization on manifolds. We exploit therefore a simulationbased retraction operator which exploits the stability of the steady state. This retraction operator can be used within stateoftheart optimizers to eliminate the equality constraint and reduce the problem dimension. As the method uses both, a simulationbased retraction operator and a stateoftheart optimizer, we will in the following refer to it as a hybrid optimization method.
The second method uses continuous analogues of local optimizers [29]. Continuous analogues are dynamical systems whose trajectories converge to a locally optimal point of an optimization problem. These dynamical systems often possess a larger basin of attraction [30] and can be solved using sophisticated numerical methods. This promises more robust convergence than simple stepsize controls used in existing optimization methods. Continuous analogues have been constructed for a series of linear and nonlinear problems [30–32]. We introduce continuous gradient descent and NewtonRaphson methods for solving optimization problems with steadystate constraints. The manifold is stabilized using a continuous retraction derived from the model. This method is purely simulationbased and will therefore be referred to as a simulationbased optimization method.
The proposed optimization methods are illustrated using a simulation example. This is followed by a rigorous evaluation of the methods and comparison with alternative methods. To this end we consider two applications: NGFinduced activation of ERK in primary sensory neurons; and Raf/MEK/ERK signaling in HeLa cells. Using novel data for the second application, new insights into Raf/MEK/ERK signaling upon release from Sphase arrest are discovered.
Methods
In this section we introduce the model class and the optimization problem. The differential geometry of the steadystate manifold is described and two optimization methods exploiting this geometry are introduced. The properties of these methods along with their implementation are discussed.
The optimization methods are used in the Results section to study Raf/MEK/ERK signaling in HeLa cells after release from Sphase arrest. The biological materials and the setups used to study this process experimentally are described below.
Mathematical modeling of perturbation experiments
with states \(x(t) \in \mathbb {R}^{n_{x}}\), observables \(y(t) \in \mathbb {R}^{n_{y}}\), parameters \(\theta \in \mathbb {R}^{n_{\theta }}\) and inputs \(u(t) \in \mathbb {R}^{n_{u}}\). The states x(t) are the concentrations of biochemical species at time t. The observables y(t) are the values of measurable quantities. The parameters θ are biochemical reaction rates, total abundances of biochemical species (in the presence of conservation relations) and experimental parameters (e.g. scaling and offset). The inputs u encode the experimental conditions applied to the system, e.g., extracellular concentration of ligands. To ensure existence and uniqueness of the solution of (1), the vector field \(f:\mathbb {R}^{n_{x}} \times \mathbb {R}^{n_{\theta }} \times \mathbb {R}^{n_{u}} \to \mathbb {R}^{n_{x}}\) is assumed to be Lipschitz continuous. The mapping \(h:\mathbb {R}^{n_{x}} \times \mathbb {R}^{n_{\theta }} \times \mathbb {R}^{n_{u}} \to \mathbb {R}^{n_{y}}\) describes the observation process and the mapping \(x_{0}:\mathbb {R}^{n_{\theta }} \to \mathbb {R}^{n_{x}}\) provides the initial conditions.
The stability of \({x_{s}^{e}}\) can be assessed using Lyapunov theory [33]. We denote the collection of all parameterstate pairs (θ,x _{0}) which fulfill the steadystate constraint (2) for a given input u as the manifold of steady states. For simplicity, we assume that (2) possesses a unique, exponentially stable steady state for every combination of parameters θ and inputs u. In this case, there exists a function \(x_{s}: \mathbb {R}^{n_{\theta }} \times \mathbb {R}^{n_{u}} \rightarrow \mathbb {R}^{n_{x}}\) which maps the parameters to the corresponding steady state, i.e., x(0)=x _{ s }(θ,u) (see Fig. 1b). An analytical expression of the function x _{ s }(θ,u) is usually not available.
with \(y^{e}\left (t,\theta,{x_{s}^{e}}\right)\) denoting the solution of (1) for input u=u ^{ e }, parameters θ and initial condition \({x_{s}^{e}}\) at time t. The measurement noise \(\epsilon _{ij}^{e}\) is assumed to be normally distributed, \(\epsilon _{ij}^{e} \sim \mathcal {N}\left (0,\left (\sigma _{ij}^{e}\right)^{2}\right)\), the methods presented in the following are however not limited to this case.
Parameter estimation
is used. This likelihood function depends on θ and \({x_{s}^{1}},\ldots,{x_{s}^{E}}\), variables which are coupled via the steadystate constraint (2).
in which the objective function denotes the negative loglikelihood function, \(J\!\!\left (\!\theta \!,{x_{s}^{1}},\ldots,{x_{s}^{E}}\right) \,=\, \! \log p\!\left (\!\mathcal {D}\theta \!,{x_{s}^{1}},\ldots,{x_{s}^{E}}\right)\). The solution of (6) provides parameterstate pairs \(\left (\hat \theta,\hat {x}_{s}^{e}\right)\) on the steadystates manifold which maximize the likelihood function (5). For these pairs it holds that \(\hat {x}_{s}^{e} = x_{s}\left (\hat \theta,{u_{c}^{e}}\right)\).
In general, optimization problem (6) is nonlinear and possesses local minima. To compute the optimum of (6), we employ multistart local optimization. This approach proved to be efficient in a variety of related problems (see, e.g., [15, 34]). Furthermore, sophisticated local optimizers allow for the consideration of nonlinear equality constraints (see [35] and references therein). The consideration of nonlinear equality constraints is not possible for most evolutionary and genetic algorithms [36], particle swarm optimizers [37], simulated annealing [38] and hybrid optimizers [19, 39]. Alternative methods are metaheuristics which combine ideas from local and global optimization methods [40], facilitating the analysis of optimization problems with nonlinear constraints and multiple local minima.

Unconstrained optimization method: An analytical expression of the steady state as a function of the parameter, x _{ s }(θ), is used to eliminate the constraint and x _{0} from optimization problem (6). This yields the reduced optimization problem$$ \begin{aligned} \min_{\theta} &\; J\left(\theta,x_{s}\left(\theta,{u_{c}^{1}}\right),\ldots,x_{s}\left(\theta,{u_{c}^{E}}\right)\right) \end{aligned} $$(7)
which does not possess any equality constraints. While this method is rarely applicable – analytical expressions for the x _{ s }(θ) are difficult to compute – it provides a gold standard.

Constrained optimization method: An interior point optimization method is used to solve the optimization problem (6). This is the stateoftheart method and mostly used in practice.

Hybrid optimization method: The optimization problem (6) is reduced to the manifold of steady states by computing x _{ s }(θ) numerically.

Simulationbased optimization method: A dynamical system is formulated whose trajectories converge to local optima of the optimization problem (6). The dynamical system is solved using stateoftheart numerical methods.
and the gradients of the constraint. The sensitivities of the observable with respect to the parameters and initial conditions, \(\frac {\partial y}{\partial \theta }\) and \(\frac {\partial y}{\partial {x_{s}^{e}}}\), are computed using the forward sensitivity equations [41].
The optimization problem considered (6) is rather general and allows for the consideration of multiple steadystate constraints, as well as steadystate dose response curves. In the next section the geometry of the steadystate manifold is discussed. Subsequently, the hybrid optimization method and the simulationbased optimization method are introduced.
Manifold of steady states

solving a feasibility problem (find \(x_{s} \in \mathbb {R}^{n_{x}}\) subject to 0=f(x _{ s },θ,u)),

simulating the dynamical system until the steady state is reached, or

combining the simulation of the dynamical system with finetuning using the NewtonRaphson method [42].
evaluated at (θ,x _{ s }(θ,u),u). The invertibility of the Jacobian (∂ f/∂ x) follows from local exponential stability of the steady state. This result can also be obtained using the implicit function theorem.
Given an update direction Δ θ and a length r, (14) provides the steady state for parameter θ+r Δ θ up to the accuracy of the chosen ODE solver. Hence, (14) enables moves on the steadystate manifold, similar to results in [28] for other manifolds.
Hybrid optimization method
The gold standard for solving optimization problem (6) is to determine an analytical expression for the parameterdependent steady state. While such an expression is not always available, a straightforward approach is to compute the steady state numerically for given parameters θ. This is computationally more demanding than using an analytical expression, but it also yields the reduced optimization problem (7).

Step 2: The steady states \(x_{s}^{1,l+1},\ldots, x_{s}^{E,l+1}\) for the parameters θ ^{ l+1} are computed using one of the methods discussed in the Manifold of steady states section (with starting points \(x_{s}^{1,l},\ldots, x_{s}^{E,l}\)). This yields the point \(\left (\theta ^{l+1},x_{s}^{1,l+1},\ldots,x_{s}^{E,l+1}\right)\) on the steadystate manifold. (Represented as a dotted arrow in Fig. 2 a.)

Step 3: The objective function \(J^{l+1} = J\left (\theta ^{l+1},{\newline } x_{s}^{1,l+1},\ldots, x_{s}^{E,l+1}\right)\) is computed for parameters θ ^{ l+1} and numerically calculated steady state \(x_{s}^{1,l+1},\ldots,x_{s}^{E,l+1}\). This objective function is provided to the local optimizer. (Not represented in Fig. 2 a.)
with the sensitivity of the steady states with respect to the parameters, \((\partial {x_{s}^{e}}/\partial \theta) = S(\theta,{x_{s}^{e}},{u_{c}^{e}})\), as defined in (12).
The proposed hybrid optimization method possesses all properties and options of the employed local optimizer. In addition, the retraction accuracy ε _{tol} of the convergence criteria \(f(x_{s}^{e,l},\theta ^{l},{u_{c}^{e}})_{2} < \epsilon _{\text {tol}}\) has to be selected.
Simulationbased optimization method
with the steadystate sensitivity \(S\left (\theta,{x_{s}^{e}},{u_{c}^{e}}\right)\). Initialization of (16) in a point \(\left (\theta _{0},x_{s,0}^{1},\ldots,x_{s,0}^{E}\right)\) fulfilling (2) yields a trajectory evolving on the steadystate manifold, along which the objective function decreases.
in which (∂ f/∂ x)^{+} denotes the Moore–Penrose pseudoinverse of (∂ f/∂ x) at \((\theta,{x_{s}^{e}},{u_{c}^{e}})\). On the steadystate manifold, the Jacobian is invertible and we recover the standard steadystate sensitivity. For a large retraction factor λ≫0, the state \((\theta,{x_{s}^{1}},\ldots,{x_{s}^{E}})\) is retracted quickly to the steadystate manifold.

Gradient descent: \(g(\theta,x_{s}) =  \frac {dJ}{d\theta }\) and

Newtontype descent: \(g(\theta,x_{s}) =  \left (F + \mu I \right)^{1} \frac {dJ}{d\theta }\).
The Newtontype descent exploits the Fisher Information Matrix F [43]. The Fisher Information matrix is an approximation to the Hessian of the objective function. This approximation can be computed from firstorder sensitivities and is positive semidefinite.
For the gradient descent we established local exponential stability of local optima for appropriate choice of λ for a broad class of problems (see Additional file 1: Section 1). This implies that the trajectories of the system converge to the local optima of the optimization problem (6). We expect that a similar result can be derived for the Newtontype descent. Though this is not yet available, we included the Newton type descent as we expect – similar to classical optimizers – faster convergence.
For local optimization of (6), the dynamical system (17) has to be simulated for r→∞. For this, implicit methods with adaptive stepsize selection and error control should be employed as (17) might be stiff. Appropriate numerical methods are implemented among others in MATLAB and the SUNDIALS package [44]. These simulations are stopped as soon as the convergence criterion max{∥d θ/d r∥,∥d x _{ s }/d r∥}<ε _{tol} is met.
Implementation
The different methods are implemented in MATLAB and provided in an Additional file 2. The local optimization is performed using the MATLAB routine fmincon.m. fmincon is supplied with the objective function value and the values of the constraint, as well as the respective analytical derivatives. The continuous analogue used for simulationbased optimization is simulated using the MATLAB ODE solver ode15s. The computationally intensive simulation of the perturbation experiments and the numerical calculation of the steady state is performed using the SUNDIALS package [44] which was accessed using the MATLAB Toolbox AMICI (https://github.com/AMICIdeveloper/AMICI). Default settings are used for fmincon.m and the simulation routines unless stated otherwise. The convergence tolerance for the hybrid optimization method is set to ε _{tol}=10^{−9} in the simulation example and the first application example and to ε _{tol}=10^{−13} in the second application example. For the simulationbased optimization methods it is set to ε _{tol}=10^{−6}.
Experimental data
To evaluate the performance of the proposed methods, we use two application examples with experimental data. In the first application example, we consider a dataset for NGFinduced ERK signaling in primary sensory neurons which was published by Andres et al. [45]. In the second application example, we use novel data for Raf/MEK/ERK signaling in HeLa cells after release from Sphase arrest. This novel experimental data for Raf/MEK/ERK signaling in HeLa cells was acquired using the following methods.
Cell culture. HeLa cells were obtained from the American Type Culture Collection (Manassas, VA). Cells were maintained in RPMI 1640 supplemented with 10 % fetal bovine serum.
Cell synchronization at the G1/S border. HeLa cells were synchronized at the G1/S border using an aphidicolin treatment. In brief, cells grown on Petri dishes were incubated in medium supplemented with aphidicolin (0,3 g/ml; Calbiochem) for 18 h. Afterward, cells were released from the Sphase arrest by washing with serumfree medium and were refed with growth medium.
Protein extraction of cells. Wholecell extracts were obtained by solubilizing cells in hot protein sample buffer (95 °C). After 10 min of incubation at 95 °C, extracts were placed on ice and centrifuged (16,000 μg, 15 min, and 4 °C). Samples were subjected to SDSPAGE.
Western blotting. Equal amounts of proteins were separated by SDSPAGE and blotted onto nitrocellulose membranes (Pall, Dreieich, Germany). After blocking with 0.5 % blocking reagent (Roche Diagnostics), filters were probed with specific antibodies. Proteins were visualized with IRdyecoupled secondary antibodies using the LiCOR Odyssey system. Protein bands were quantified using ImageJ.
Antibodies Commercially available antibodies used in this study were: antiERK rabbit polyclonal, antiphosphop44/42 MAPK (ERK1/2) (Thr202/Tyr204) rabbit polyclonal, and antiphosphoMEK1/2 (Ser217/221) rabbit polyclonal (all from Cell Signaling). The tubulinspecific mouse monoclonal antibody was from Millipore. The IRdyecoupled secondary antibodies were from LiCOR Biosciences.
Results
In the following, we illustrate the behavior of the proposed optimization method. Furthermore, the performance of the proposed methods will be compared to standard constrained and unconstrained optimization methods. For this purpose, we consider a simulation example for which the ground truth is known. Furthermore, we test the methods on two application examples using real data.
Simulation example: Conversion process
In this section, we illustrate the proposed optimization methods by studying parameter estimation for a conversion process from steadystate data. Conversion processes are among the most common motifs in biological systems, therefore particularly interesting, and provide a simple test case.
in which \(\bar {y}(t_{j})\) denotes the measured concentration of A at time point t _{ j } and y(t _{ j },θ,x _{ s }) denotes the solution of (20) for initial conditions x(0)=x _{ s } and input u=0.4 at time point t _{ j }.
Illustration of hybrid optimization method
are derived. The local sensitivities depend merely on derivatives of the vector field and can be computed without knowledge of an analytical expression of the steady state.
Illustration of simulationbased optimization method
with initial conditions θ _{1}(0)=θ _{1,0}, θ _{2}(0)=θ _{2,0} and x _{ s }(0)=x _{ s,0}. It can be verified that the objective function J is locally strictly convex in θ – the parameters are locally identifiable – and that the model (20) is asymptotically stable. Accordingly, system (25) converges to a local optimum of the constrained optimization problem (22) (see Additional file 1: Theorem 1).
Application example 1: NGFinduced ERK signaling in primary sensory neurons
To evaluate and compare existing and proposed local optimization methods for problems with steadystate constraints, we analyze NGFinduced ERK phosphorylation in primary sensory neurons. Primary sensory neurons are among others used to investigate pain sensitization in response to inflammation. During inflammation a cocktail of stimuli is present, including NGF. NGF binds to cellular receptors and induces the ERK phosphorylation [45]. This modulates neuronal activity by triggering ion channel phosphorylation and protein expression [46].
Growthfactor induced ERK signaling is a potential target for novel pain therapies [47] and therefore of high practical relevance. In addition, this application is wellsuited for the evaluation of the methods as NGF doseresponse curves at late time points have been recorded. These data provide multiple steadystate constraints for the thorough assessment of the methods. In the following, we will compare the performance of unconstrained, constrained, hybrid and simulationbased optimization in the presence of multiple steadystate constraints.
Experimental data for NGFinduced ERK phosphorylation
ERK phosphorylation in response to different concentrations of NGF was previously quantified using quantitative automated microscopy [45]. This technique provides singlecell data from which population average data can be derived. These population average data are highly reproducible and quantitative but provide merely the relative ERK phosphorylation in comparison to the control as no calibration curve is employed. The unknown scaling constant is denoted by s.
Mathematical model of NGFinduced ERK phosphorylation
to describe the activities of the NGF receptor TrkA (x _{1}) and ERK phosphorylation (x _{2}) in response to NGF stimulation, u= [NGF]_{0}. This model possesses a minimal number of model 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 Additional file 1: Section 2 and the original publication [34]. The experimental noise ε is assumed to be normally distributed with the unknown variance σ ^{2}, \(\epsilon \sim \mathcal {N}(0,\sigma ^{2})\).
This steady state exists for all positive parameters and is exponentially stable.
Parameter estimation problem with multiple steadystate constraints
In this study, the unknown parameters \(\theta \in \mathbb {R}_{+}^{7}\) and the states x _{ s,1} and x _{ s,2} for each considered input of NGF are inferred from published dose response data [45] using ML estimation. The dataset contains 6 different NGF doses, yielding an optimization problem with 7+2·6=19 optimization variables and 2·6=12 nonlinear equality constraints. This nonlinear optimization problem is solved using multistart local optimization. The local optimization is performed using unconstrained, constrained and hybrid optimization as well as simulationbased optimization using gradient and Newtontype descent. Bounds and scales for the parameters are provided in the Additional file 1: Table S1.
The convergence properties of unconstrained, hybrid and simulationbased optimization are comparable
To assess the convergence of the optimization method, we sort and visualize the objective function values achieved in the individual optimizer starts (Fig. 6 a). In addition, we determine the percentage of converged starts. A start is considered to be converged if the final point cannot be rejected compared to the ML estimate using the likelihood ratio test with a significance level of 0.05.
As expected, we find that the gold standard – the unconstrained optimization method – shows the best convergence properties. It converges in 75 % of the starts to the global optimum. A similar convergence is achieved by the proposed methods, hybrid optimization and simulationbased optimization using gradient descent. The third proposed method – simulationbased optimization using Newtontype descent – displays intermediate convergence properties (60 % of the starts converged to the global optimum). The stateoftheart method – constrained optimization – exhibits the poorest convergence. It converges in 45 % of the starts. Hence, the proposed optimization methods are superior to constrained optimization regarding convergence to the global optimum.
Beyond differences in the convergence to the global optimum, the convergence to local optima differs. The results of unconstrained, constrained and hybrid optimization reveal three local optima. The local optima with the worst objective function values are hardly found using simulationbased optimization, indicating altered regions of attraction.
Hybrid and simulationbased optimization provide reliable estimates of the steady states
The individual optimization methods enforce the steadystate constraints differently. What all methods have in common is that the steadystate constraint f(x _{ s },θ,u _{ s })=0 is relaxed to a constraint on the norm of the vector field, i.e., f(x _{ s },θ,u _{ s })_{2}<ε _{ f }. Accordingly, parameterstate pairs returned by the optimization methods usually do not fulfill steadystate constraints exactly. Different optimization methods might even achieve different accuracies. In addition, a bound for the difference of the estimated steady state x _{ s } for a parameter θ and the true steady state x _{ s }(θ), Δ x _{ s }=x _{ s }−x _{ s }(θ), is usually not available.
We studied the relation of the solver indicating convergence based on the vector field (f_{2}<ε _{ f }) and the difference of the estimated to the analytical steadystate being small (Δ x _{ s }_{2}<ε _{ x }) for the different optimization methods. In our opinion a good optimizer should achieve equivalence of the two criteria. This would mean that enforcing the constraint of the vector field ensures a good approximation of the steady state. The result is depicted in Fig. 6 b for a tolerance of 10^{−6}.
The unconstrained optimization uses an analytical expression of the steady state and therefore the two criteria are identical. Hybrid and simulationbased optimization also achieved a good agreement of both criteria, with ∼85 %. In ∼15 % of the cases, the solver indicates convergence based on the vector field constraint but the steadystate estimate is off (Δ x _{ s }_{2}>ε _{ x }). The precise percentage depends heavily on the retraction factor λ for the simulationbased optimization method. For the constrained optimization, all possible combinations are observed and the two criteria agree in merely 55 % of the runs. In summary, the results indicate that the proposed methods provide reliable estimates for the steady states while constrained optimization yields many inconsistent parameterstate pairs.
Hybrid and simulationbased optimization are faster than constrained optimization
A key performance metric for local optimizers is the average computation time per converged start. This metric summarizes convergence properties (Fig. 6 a and b) and computation times for individual starts (Fig. 6 c). It is computed by dividing the overall computation time for the multistart optimization by the number of converged starts. This measure of optimizer performance is depicted in Fig. 6 d.
Unconstrained optimization using a (usually not available) analytical expressions for steady states is most efficient. The individual runs are fast and the percentage of converged starts is high. Hybrid and simulationbased methods are roughly 10 times slower but these methods can be applied if analytical expressions for steady states are not available. Furthermore, these methods are 1.5 times faster than constrained optimization due to the improved convergence rate. Additionally, the fit to the data for the optimal parameters is convincing (Fig. 6 e). Accordingly, we conclude that hybrid and simulationbased optimization are promising approaches in the presence of multiple steadystate constraints.
Application 2: Raf/MEK/ERK signaling in HeLa cells after release from Sphase arrest
In this section, we study Raf/MEK/ERK signaling in HeLa cells after release from Sphase arrest. Experimental studies revealed that cellcycle is, among others, controlled by Raf/MEK/ERK signaling [50, 51]. The signaling dynamics in different cellcycle phases as well as the cellcycledependent relevance of feedback mechanisms [52] are however still not completely unraveled although a more thorough understanding could provide valuable insights into treatment resistance [52]. Using the new data and model selection we study the relevance of negative feedback from phosphoERK to Raf activation during G1/S phase transition.
In addition to its biological relevance, the Raf/MEK/ ERK pathway is wellsuited for the evaluation of the proposed optimization methods and the comparison to stateoftheart methods. The pathway is nonlinear, yielding a nonlinear and nonconvex optimization problem. Furthermore, we will consider a synchronized cell population which reached a steady state before the start of the experiment. Accordingly, a steadystate constraint has to be enforced and fitted along with timeresolved data for perturbation experiments.
Experimental data for Raf/MEK/ERK signaling after release from Sphase arrest
To study the Raf/MEK/ERK pathway, HeLa cells were synchronized at the G1/S border using an aphidicolin treatment. After synchronization was achieved, aphidicolin was removed and the dynamics of phosphoMEK and phosphoERK were quantified using Western blotting. This was repeated after treatment with Sorafenib and UO126 to explore the dynamic range of the pathway. Sorafenib is an inhibitor of Raf kinases [53] and UO126 is a highly selective inhibitor of MEK [54].
As Western blots are merely semiquantitive, they provide the relative activity of phosphoMEK and phosphoERK at different time points and under different conditions. The unknown scaling constants differ between blots and measured species. For a detailed discussion of characteristics of Western blot data we refer to [55].
Mathematical model for Raf/MEK/ERK signaling after release from Sphase arrest
(as proposed in the Data2Dynamics toolbox [58]). The effects of Sorafenib and UO126 are captured by a reduction in the kinase activity of pRaf and pMEK (R_{4} and R_{6}).

Inhibition of Raf phosphorylation by pERK: \(\xi (t) = \frac {K_{1}}{K_{1} + [\text {pERK}]}\)

No inhibition: ξ(t)=1
for the relative phosphorylation levels x _{1}= [pRaf]/[Raf]_{0}, x _{2}= [pMEK]/[MEK]_{0} and x _{3}= [pERK]/[ERK]_{0} and the input u=([sora],[UO126]). The model for the relative phosphorylation levels 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 nonidentifiabilities and reduces the number of parameters. Each Western blot, indexed by b=1,…,4, provides timeresolved data for y _{1,b } and y _{2,b } for a combination of different experimental conditions. The measurement noise is assumed to be normally distributed and its variance is estimated from the experimental data. As all parameters are nonnegative, a logparameterization is used for parameter estimation [15]. The states of the reformulated model are between 0 and 1. Details regarding parameters and initial conditions are provided in the Additional file 1: Table S2.
In addition to the kinetic, scaling and noise parameters, the initial conditions of the models for H1 and H2 are unknown. As the cells are however arrested in Sphase with k _{1,max}(0)=k _{1,0} and u=0, the initial conditions are the corresponding steady states. After significant manual preprocessing of the steadystate constraints, analytical expressions x _{ s }(θ) for the steadystates as a function of the other parameters could be calculated with symbolic math toolboxes. These analytical expressions are provided in the Additional file 1: Equation (10), (11).
Parameter estimation problem with multiple perturbation datasets
Hybrid and simulationbased optimization outperforms constrained optimization
Unconstrained optimization using the analytical expression for the steady state – the gold standard – converged in ∼50 % of the starts (Fig. 8 a). Hybrid and simulationbased optimization methods achieved a percentage of converged starts comparable to the gold standard (40–60 %), but without requiring an analytical expression for the steady state. Constrained optimization – the stateoftheart – converged in less than 10 % of the starts, resulting in a relatively large computation time per converged start (Fig. 8 b). Even though hybrid and simulationbased optimization were slower than the gold standard, they were more than 10 times faster than constrained optimization. Hence, the proposed optimization methods also outperform constrained optimization for this problem.
A detailed comparison of the proposed methods revealed that simulationbased optimization using gradient descent achieved the highest percentage of converged starts. Hybrid optimization required however fewer simulations of the perturbation experiments – the timeconsuming step – rendering this method computationally more efficient. Simulationbased optimization using Newtontype descent was the least efficient of the proposed methods. This might be related to the challenges in tuning the regularization parameters.
Model selection reveals importance of negative feedback
The model with negative feedback (H1) fits the experimental data (Fig. 8 c). It captures the transient phosphorylation of MEK and ERK after release from Sphase arrest, the reduced ERK phosphorylation in the presence of Sorafenib and UO126. Furthermore, the increased MEK phosphorylation after UO126 treatment is explained via a decrease in the strength of the negative feedback which is caused by the reduced abundance of pERK. The model without the negative feedback loop (H2) is not able to capture the difference between the control condition and the simulation with UO126. The value of the Bayesian Information Criterion (BIC) [59] is 278.4 for the model with negative feedback (H1) and 317.4 for the model without negative feedback (H2). The difference of 39.0 indicates a strong preference for H1 [60]. The same conclusion is reached using the Akaike Information Criterion (AIC) [61]. We conclude that Raf phosphorylation is inhibited by pERK during G1/S phase transition.
To conclude, in this section we illustrated the proposed hybrid and simulationbased optimization methods. The applicability of the methods was demonstrated by studying relevant biological problems. The comparison with stateoftheart methods revealed convergence and computational efficiency. The study of Raf/MEK/ERK signaling using the methods underlined the feedback regulation of ERK phosphorylation during cell cycle progression.
Discussion
Optimization problems with steady state constraints arise in many biology applications for a wide range of models. For some models an analytical expression for the steady state can be derived and used to eliminate the steadystate constraints [18]. While this is favorable, it is not always possible. In cases in which no analytical expressions are available, the vector of optimization variables contains the unknown parameters as well as the corresponding steady states. The optimizers have to evolve on the nonlinear manifold, the set of steady states. In this manuscript, we propose a hybrid optimization method and a continuous analogue to solve optimization problems with steadystate constraints more efficiently. This simulationbased method exploits the local geometry of the steadystate manifold for optimization.
The proposed hybrid and simulationbased optimization methods are evaluated using three models for biological processes. Following a simple illustration example, an application with multiple steadystate constraints and an application with timeresolved data for multiple perturbation conditions are considered. For this rich set of scenarios we find that the hybrid and simulationbased optimization methods possess improved convergence properties in comparison to standard constrained optimization methods implemented in the MATLAB routine fmincon. We expect that the proposed methods also outperform alternative optimization routines (e.g. IPOPT [62]), this, however, remains to be analyzed. The proposed optimization methods yield convergence properties comparable to those of unconstrained optimization methods exploiting an analytical expression for the steady state. However, if analytical expressions for the steady state can be determined using available methods [18, 25, 26], unconstrained optimization should be used as the computation time is lower. The proposed methods are also applicable to a broader class of problems for which no analytical expression for the steady state is available. Furthermore, the method directly allows for multiple steadystate constraints. Unlike methods based on sequential geometric programming [63, 64], steadystate and kinetic data can be handled.
Beyond the evaluation of the proposed methods, our experimental and computational analysis of novel data for Raf/MEK/ERK signaling after release for Sphase arrest provided new insights. Parameter estimation and model comparison indicate that the negative feedback from ERK to Raf phosphorylation is also active at the G1/S border. This complements previous knowledge of stimulus, contextdependence of this stimulus [49] and its relevance for the robustness of MAPK signaling in tumor cells [52].
The implementation of the hybrid optimization method employed in this study is a simulationbased retraction operator. Alternatively, efficient and accurate schemes combining simulation and local optimization could be employed to compute steady states and sensitivities [42]. This should improve the computational efficiency further. To relax the stability assumption for the steady state, conservation relations can be incorporated in the local optimization scheme.
For the simulationbased optimization method we established local asymptotic stability of optimal points using perturbation theory. This result is however restricted to the gradienttype descent and locally convex objective functions. The latter implies local practical identifiability. The theoretical properties of the Newtontype descent and the properties in the presence of practical and structural nonidentifiability remain to be analyzed. Preliminary results and the applications suggest that in the presence of nonidentifiabilities the simulationbased optimization method always yields a point on the nonidentifiable subspace. Furthermore, the available proof shows the retraction factor λ has to be chosen large enough to ensure convergence. However, as too large λ will result in a stiff system, an intelligent choice of λ is necessary.
The simulation example and the applications possess a unique exponentially stable steady state. However, preliminary results suggest that the methods also achieve good convergence for dynamical systems with multiple stable steady states and bifurcations [65] (see Additional file 1: Section 6). The theoretical analysis of the proposed methods and a detailed performance evaluation for dynamical systems with such properties remains to be addressed.
Beyond parameter estimation, the proposed optimization methods can also accelerate practical indentifiability analysis and uncertainty quantification by speeding up optimization runs in bootstrapping uncertainty analysis [66, 67] and profile likelihood calculation [68]. In addition, Bayesian uncertainty analysis using Markov chain Monte Carlo sampling [11] can profit from an efficient initial optimization prior to the sampling. This has been shown to reduce the warmup and to improve convergence [69]. For a more detailed discussion of identifiabilty and uncertainty analysis methods, we refer to recent comparative studies [70, 71].
Conclusion
In summary, the proposed optimization methods are promising alternatives to constrained optimization for optimization problems with steadystate constraints. They are applicable to a wide range of ODEconstrained optimization problems [34, 72] and can – unlike methods which rely on an analytical expression for the steady state – be extended to PDE constrained optimization problems [73]. The availability of the MATLAB code will facilitate the application and extension of the methods, as well as the integration in toolboxes such as Data2Dynamics [58] and COPASI [74]. Accordingly, our study has a strong potential influence on the analysis of optimization problems with steadystate constraints in practice.
Abbreviations
Not applicable.
Declarations
Acknowledgements
Not applicable.
Funding
This work was supported by the German Federal Ministry of Education and Research (BMBF) within the SYSStomach project [grant number 01ZX1310B] [J.H. and F.J.T.] and the Postdoctoral Fellowship Program (PFP) of the Helmholtz Zentrum München [J.H.].
Availability of data and materials
All data generated and analyzed during this study are included in this published article and its supplementary information files.
Authors’ contributions
AF, FJT and JH developed the proposed optimization methods. AF and JH implemented and evaluated the optimization methods. AH and SR designed and conducted the experiments for Raf/MEK/ERK signaling upon cell cycle release. AF and JH developed the model of Raf/MEK/ERK signaling upon cell cycle release and fitted it to experimental data. AF, AH and JH analyzed the estimation results. AF, FJT and JH wrote the manuscript. The final manuscript has been approved by all authors.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Klipp E, Herwig R, Kowald A, Wierling C, Lehrach H. Systems Biology in Practice. Weinheim: WileyVCH Verlag GmbH & Co. KGaA; 2005.View ArticleGoogle Scholar
 Kitano H. Systems biology: a brief overview. Science. 2002; 295(5560):1662–4.View ArticlePubMedGoogle Scholar
 Kauffman S, Peterson C, Samuelsson B, Troein C. Random boolean network models and the yeast transcriptional network. Proc Natl Acad Sci U S A. 2003; 100(25):14796–9. doi:10.1073/pnas.2036429100.View ArticlePubMedPubMed CentralGoogle Scholar
 Klamt S, Haus UU, Theis F. Hypergraphs and cellular networks. PLoS Comput Biol. 2009; 5(5):1000385. doi:10.1371/journal.pcbi.1000385.View ArticleGoogle Scholar
 Klann MT, Lapin A, Reuss M. Stochastic simulation of signal transduction: Impact of the cellular architecture on diffusion. Biophys J. 2009; 96(12):5122–9. doi:10.1016/j.bpj.2009.03.049.View ArticlePubMedPubMed CentralGoogle Scholar
 Hasenauer J, Jagiella N, Hross S, Theis FJ. Datadriven modelling of biological multiscale processes. J Coupled Syst Multiscale Dyn. 2015; 3(2):101–21. doi:10.1166/jcsmd.2015.1069.View ArticleGoogle Scholar
 Klipp E, Nordlander B, Krüger R, Gennemark P, Hohmann S. Integrative model of the response of yeast to osmotic shock. Nat Biotechnol. 2005; 23(8):975–82.View ArticlePubMedGoogle Scholar
 Schöberl B, Pace EA, Fitzgerald JB, Harms BD, Xu L, Nie L, Linggi B, Kalra A, Paragas V, Bukhalid R, Grantcharova V, Kohli N, West KA, Leszczyniecka M, Feldhaus MJ, Kudla AJ, Nielsen UB. Therapeutically targeting ErbB3: a key node in ligandinduced activation of the ErbB receptor–PI3K axis. Sci Signal. 2009; 2(77):31.Google Scholar
 Bachmann J, Raue A, Schilling M, Böhm ME, Kreutz C, Kaschek D, Busch H, Gretz N, Lehmann WD, Timmer J, Klingmüller U. Division of labor by dual feedback regulators controls JAK2/STAT5 signaling over broad ligand range. Mol Syst Biol. 2011;7:516.Google Scholar
 Tarantola A. Inverse problem theory and methods for model parameter estimation. Philadelphia: SIAM; 2005.View ArticleGoogle Scholar
 Xu TR, Vyshemirsky V, Gormand A, von Kriegsheim A, Girolami M, Baillie GS, Ketley D, Dunlop AJ, Milligan G, Houslay MD, Kolch W. Inferring signaling pathway topologies from multiple perturbation measurements of specific biochemical species. Sci Signal. 2010; 3(113):20. doi:10.1126/scisignal.2000517.View ArticleGoogle Scholar
 Molinelli EJ, Korkut A, Wang W, Miller ML, Gauthier NP, Jing X, Kaushik P, He Q, Mills G, Solit DB, Pratilas CA, Weigt M, Braunstein A, Pagnani A, Zecchina R, Sander C. Perturbation biology: Inferring signaling networks in cellular systems. PLoS Comput Biol. 2013; 9(12):1003290.View ArticleGoogle Scholar
 Banga JR. Optimization in computational systems biology. BMC Syst Biol. 2008; 2(47):1–7.Google Scholar
 Chou I, Voit EO. Recent developments in parameter estimation and structure identification of biochemical and genomic systems. Math Biosci. 2009; 219(2):57–83. doi:10.1016/j.mbs.2009.03.002.View ArticlePubMedPubMed CentralGoogle Scholar
 Raue A, Schilling M, Bachmann J, Matteson A, Schelke M, Kaschek D, Hug S, Kreutz C, Harms BD, Theis FJ, Klingmüller U, Timmer J. Lessons learned from quantitative dynamical modeling in systems biology. PLoS ONE. 2013; 8(9):74335. doi:10.1371/journal.pone.0074335.View ArticleGoogle Scholar
 Moles CG, Mendes P, Banga JR. Parameter estimation in biochemical pathways: a comparison of global optimization methods. Genome Res. 2003; 13:2467–74.View ArticlePubMedPubMed CentralGoogle Scholar
 Villaverde AF, Henriques D, Smallbone K, Bongard S, Schmid J, CicinSain D, Crombach A, SaezRodriguez J, Mauch K, BalsaCanto E, Mendes P, Jaeger J, Banga JR. BioPreDynbench: a suite of benchmark problems for dynamic modelling in systems biology. BMC Syst Biol. 2015;9(8). doi:10.1186/s1291801501444.
 Rosenblatt M, Timmer J, Kaschek D. Customized steadystate constraints for parameter estimation in nonlinear ordinary differential equation models. Front. 2016; 4:41.Google Scholar
 Vaz A, Vicente L. A particle swarm pattern search method for bound constrained global optimization. J Global Optim. 2007; 39(2):197–219. doi:10.1007/s1089800791335.View ArticleGoogle Scholar
 Egea JA, Henriques D, Cokelaer T, Villaverde AF, MacNamara A, Danciu DP, Banga JR, SaezRodriguez J. MEIGO: an opensource software suite based on metaheuristics for global optimization in systems biology and bioinformatics. BMC Bioinf. 2014; 15(136):1–9.Google Scholar
 King EL, Altman C. A schematic method of deriving the rate laws for enzymecatalyzed reactions. J Phys Chem. 1956; 60(10):1375–8. doi:10.1021/j150544a010.View ArticleGoogle Scholar
 CornishBowden A. An automatic method for deriving steadystate rate equations. Biochem J. 1977; 165(1):55–9.View ArticlePubMedPubMed CentralGoogle Scholar
 Chou KC. Applications of graph theory to enzyme kinetics and protein folding kinetics. Biophys Chem. 1990; 35(1):1–24. doi:10.1016/03014622(90)80056D.View ArticlePubMedGoogle Scholar
 Feliu E, Wiuf C. Variable elimination in chemical reaction networks with massaction kinetics. SIAM J Appl Math. 2012; 72(4):959–81. doi:10.1137/110847305.View ArticleGoogle Scholar
 Halasz A, Lai HJ, McCabe Pryor M, Radhakrishnan K, Edwards J. Analytical solution of steadystate equations for chemical reaction networks with bilinear rate laws. IEEE/ACM Trans Comput Biol Bioinformatics. 2013; 10(4):957–69. doi:10.1109/TCBB.2013.41.View ArticleGoogle Scholar
 Loriaux PM, Tesler G, Hoffmann A. Characterizing the relationship between steady state and response using analytical expressions for the steady states of mass action models. PLoS Comput Biol. 2013; 9(2):1002901. doi:10.1371/journal.pcbi.1002901.View ArticleGoogle Scholar
 Bertsekas DP. Nonlinear Programming, 2nd edn. Belmont: Athena Scientific; 1999.Google Scholar
 Absil PA, Mahony R, Sepulchre R. Optimization Algorithms on Matrix Manifolds. Princeton, New Jersey: Princeton University Press; 2007.Google Scholar
 Kose T. Solutions of saddle value problems by differential equations. Econometrica. 1956; 24(1):59–70.View ArticleGoogle Scholar
 Tanabe K. Global analysis of continuous analogues of the LevenbergMarquardt and NewtonRaphson methods for solving nonlinear equations. Ann Inst Statist Math. 1985; 37(Part B):189–203.View ArticleGoogle Scholar
 Tanabe K. Continuous NewtonRaphson method for solving an underdetermined system of nonlinear equations. Nonlinear Anal Theory Methods Appl. 1979; 3(4):495–503. doi:0.1016/0362546X(79)900646.View ArticleGoogle Scholar
 Dürr HB, Ebenbauer C. A smooth vector field for saddle point problems. In: Proceedings of the 50th Conference on Decision and Control (CDC 2011). Orlando, Florida, USA: IEEE: 2011. p. 4654–660.Google Scholar
 Khalil HK. Nonlinear Systems, 3rd edn. Upper Saddle River, New Jersey: Prentice Hall; 2002.Google Scholar
 Hasenauer J, Hasenauer C, Hucho T, Theis FJ. ODE constrained mixture modelling: a method for unraveling subpopulation structures and dynamics. PLoS Comput Biol. 2014; 10(7):1003686. doi:10.1371/journal.pcbi.1003686.View ArticleGoogle Scholar
 Schittkowski K. A robust implementation of a sequential quadratic programming algorithm with successive error restoration. Optim Lett. 2011; 5(2):283–96. doi:10.1007/s1159001002079.View ArticleGoogle Scholar
 Bäck T. Evolutionary algorithms in theory and practice: evolution strategies, evolutionary programming, genetic algorithms. New York and Oxford: Oxford University Press; 1996.Google Scholar
 Yang XS. Natureinspired Metaheuristic Algorithms, 2nd edn. Bristol, UK: Luniver Press; 2010.Google Scholar
 Kirkpatrick S, Gelatt Jr. CD, Vecchi MP. Optimization by simulated annealing. Science. 1983; 220(4598):671–80.View ArticlePubMedGoogle Scholar
 BalsaCanto E, Peifer M, Banga JR, Timmer J, Fleck C. Hybrid optimization method with general switching strategy for parameter estimation. BMC Syst Biol. 2008; 2(26):1–9.Google Scholar
 RodriguezFernandez M, Egea JA, Banga JR. Novel metaheuristic for parameter estimation in nonlinear dynamic biological systems. BMC Bioinf. 2006; 7(483):1–18.Google Scholar
 Raue A. Quantitative dynamic modeling: Theory and application to signal transduction in the erythropoietic system. Phd. thesis: AlbertLudwigsUniversität Freiburg im Breisgau; 2013.Google Scholar
 Shiraishi F, Yoshida E, Voit EO. An efficient and very accurate method for calculating steadystate sensitivities in metabolic reaction systems. IEE/ACM Trans Comp Biol Bioinf. 2014; 11(6):1077–86.View ArticleGoogle Scholar
 Faller D, Klingmüller U, Timmer J. Simulation methods for optimal experimental design in systems biology. Simul. 2003; 79(12):717–25. doi:10.1177/0037549703040937.View ArticleGoogle Scholar
 Hindmarsh AC, Brown PN, Grant KE, Lee SL, Serban R, Shumaker DE, Woodward CS. SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers. ACM T Math Softw. 2005; 31(3):363–96.View ArticleGoogle Scholar
 Andres C, Meyer S, Dina OA, Levine JD, Hucho T. Quantitative automated microscopy (QuAM) elucidates growth factor specific signalling in pain sensitization. Mol Pain. 2010; 6(98):1–16. doi:10.1186/17448069698.Google Scholar
 Nicol GD, Vasko MR. Unraveling the story of NGFmediated sensitization of nociceptive sensory neurons: ON or OFF the Trks?Mol Interv. 2007; 7(1):26–41.View ArticlePubMedGoogle Scholar
 Andres C, Hasenauer J, Ahn HS, Joseph EK, Theis FJ, Allgöwer F, Levine JD, DibHajj SD, Waxman SG, Hucho T. Wound healing growth factor, basic FGF, induces Erk1/2 dependent mechanical hyperalgesia. Pain. 2013; 154(10):2216–26. doi:10.1016/j.pain.2013.07.005.View ArticlePubMedGoogle Scholar
 Fujioka A, Terai K, Itoh RE, Aoki K, Nakamura T, Kuroda S, Nishida E, Matsuda M. Dynamics of the Ras/ERK MAPK cascade as monitored by fluorescent probes. J Biol Chem. 2006; 281(13):8917–26.View ArticlePubMedGoogle Scholar
 Santos SDM, Verveer PJ, Bastiaens PIH. Growth factorinduced MAPK network topology shapes Erk response determining PC12 cell fate. Nat Cell Biol. 2007; 9(3):324–30. doi:10.1038/ncb1543.View ArticlePubMedGoogle Scholar
 Zhang W, Liu HT. MAPK signal pathways in the regulation of cell proliferation in mammalian cells. Cell Res. 2002; 12(1):9–18. doi:10.1038/sj.cr.7290105.View ArticlePubMedGoogle Scholar
 Chambard JC, Lefloch R, Pouysségur J, Lenormand P. ERK implication in cell cycle regulation. Biochim Biophys Acta. 2007; 1773(8):1299–310.View ArticlePubMedGoogle Scholar
 FritscheGuenther R, Witzel F, Sieber A, Herr R, Schmidt N, Braun S, Brummer T, Sers C, Blüthgen N. Strong negative feedback from Erk to Raf confers robustness to MAPK signalling. Mol Syst Biol. 2011;7(489). doi:http://dx.doi.org/10.1038/msb.2011.27.
 Wilhelm SM, Adnane L, Newell P, Villanueva A, Llovet JM, Lynch M. Preclinical overview of sorafenib, a multikinase inhibitor that targets both Raf and VEGF and PDGF receptor tyrosine kinase signaling. Mol Cancer Ther. 2008; 7(10):3129–140. doi:10.1158/15357163.MCT080013.View ArticlePubMedGoogle Scholar
 Favata MF, Horiuchi KY, Manos EJ, Daulerio AJ, Stradley DA, Feeser WS, Van Dyk DE, Pitts WJ, Earl RA, Hobbs F, Copeland RA, Magolda RL, Scherle PA, Trzaskos JM. Identification of a novel inhibitor of mitogenactivated protein kinase kinase. J Biol Chem. 1998; 273(29):18623–32.View ArticlePubMedGoogle Scholar
 Kreutz C, Bartolome Rodriguez MM, Maiwald T, Seidl M, Blum HE, Mohr L, Timmer J. An error model for protein quantification. Bioinf. 2007; 23(20):2747–53.View ArticleGoogle Scholar
 Kholodenko BN. Untangling the signalling wires. Nat Cell Biol. 2007; 9(3):247–9. doi:10.1038/ncb0307247.View ArticlePubMedPubMed CentralGoogle Scholar
 Kholodenko BN. Negative feedback and ultrasensitivity can bring about oscillations in the mitogenactivated protein kinase cascades. Eur J Biochem. 2000; 267(6):1583–8.View ArticlePubMedGoogle Scholar
 Raue A, Steiert B, Schelker M, Kreutz C, Maiwald T, Hass H, Vanlier J, Tönsing C, Adlung L, Engesser R, Mader W, Heinemann T, Hasenauer J, Schilling M, Höfer T, Klipp E, Theis FJ, Klingmüller U, Schöberl B, Timmer J. Data2Dynamics: a modeling environment tailored to parameter estimation in dynamical systems. Bioinformatics. 2015. doi:10.1093/bioinformatics/btv405.
 Schwarz G. Estimating the dimension of a model. Ann Statist. 1978; 6(2):461–4.View ArticleGoogle Scholar
 Burnham KP, Anderson DR. Model selection and multimodel inference: a practical informationtheoretic approach, 2nd edn. New York, NY: Springer; 2002.Google Scholar
 Akaike H. On the likelihood of a time series model. The Statistician. 1978; 27(3/4):217–35.View ArticleGoogle Scholar
 Wächter A, Biegler LT. On the implementation of an interiorpoint filter linesearch algorithm for largescale nonlinear programming. Math Program. 2006; 106(1):25–57. doi:10.1007/s101070040559y.View ArticleGoogle Scholar
 Pozo C, MarínSanguino A, Alves R, GuillénGosálbez G, Jiménez L, Sorribas A. Steadystate global optimization of metabolic nonlinear dynamic models through recasting into powerlaw canonical models. BMC Syst Biol. 2011;5(137). doi:10.1186/175205095137.
 Xu G. Steadystate optimization of biochemical systems through geometric programming. Eur J Oper Res. 2013; 225(1):12–20. doi:10.1016/j.ejor.2012.07.026.View ArticleGoogle Scholar
 Müller J, Kuttler C. Methods and models in mathematical biology, 1st edn. Berlin / Heidelberg: Springer; 2015.View ArticleGoogle Scholar
 Joshi M, SeidelMorgenstern A, Kremling A. Exploiting the bootstrap method for quantifying parameter confidence intervals in dynamical systems. Metabolic Eng. 2006; 8:447–55.View ArticleGoogle Scholar
 Kirk PDW, Stumpf MPH. Gaussian process regression bootstrapping: exploring the effects of uncertainty in time course data. Bioinformatics. 2009; 25(10):1300–6.View ArticlePubMedPubMed CentralGoogle Scholar
 Raue A, Kreutz C, Maiwald T, Bachmann J, Schilling M, Klingmüller U, Timmer J. Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinf. 2009; 25(25):1923–9.View ArticleGoogle Scholar
 Hug S, Raue A, Hasenauer J, Bachmann J, Klingmüller U, Timmer J, Theis FJ. Highdimensional Bayesian parameter estimation: Case study for a model of JAK2/STAT5 signaling. Math Biosci. 2013; 246(2):293–304. doi:10.1016/j.mbs.2013.04.002.View ArticlePubMedGoogle Scholar
 Raue A, Kreutz C, Theis FJ, Timmer J. Joining forces of Bayesian and frequentist methodology: A study for inference in the presence of nonidentifiability. Phil Trans R Soc A. 2013;371(1984). doi:10.1098/rsta.2011.0544.
 Fröhlich F, Theis FJ, Hasenauer J. Uncertainty analysis for nonidentifiable dynamical systems: Profile likelihoods, bootstrapping and more In: Mendes P, Dada JO, Smallbone KO, editors. Proceedings of the 12th International Conference on Computational Methods in Systems Biology (CMSB 2014), Manchester, UK. Lecture Notes in Bioinformatics. Cham: Springer: 2014. p. 61–72.Google Scholar
 Zechner C, Ruess J, Krenn P, Pelet S, Peter M, Lygeros J, Koeppl H. Momentbased inference predicts bimodality in transient gene expression. Proc Natl Acad Sci U S A. 2012; 109(21):8340–5. doi:10.1073/pnas.1200161109.View ArticlePubMedPubMed CentralGoogle Scholar
 Hock S, Hasenauer J, Theis FJ. Modeling of 2D diffusion processes based on imaging data: Parameter estimation and practical identifiability analysis. BMC Bioinf. 2013;14(Suppl 10)(S7). doi:10.1186/1471210514S10S7.
 Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, Singhal M, Xu L, Mendes P, Kummer U. COPASI – a COmplex PAthway SImulator. Bioinf. 2006; 22:3067–74.View ArticleGoogle Scholar