Ordinary differential equations are widely-used in the field of systems biology and chemical engineering to model chemical reaction networks. Numerous techniques have been developed to estimate parameters like rate constants, initial conditions or steady state concentrations from time-resolved data. In contrast to this countable set of parameters, the estimation of entire courses of network components corresponds to an innumerable set of parameters.
The approach presented in this work is able to deal with course estimation for extrinsic system inputs or intrinsic reactants, both not being constrained by the reaction network itself. Our method is based on variational calculus which is carried out analytically to derive an augmented system of differential equations including the unconstrained components as ordinary state variables. Finally, conventional parameter estimation is applied to the augmented system resulting in a combined estimation of courses and parameters.
The combined estimation approach takes the uncertainty in input courses correctly into account. This leads to precise parameter estimates and correct confidence intervals. In particular this implies that small motifs of large reaction networks can be analysed independently of the rest. By the use of variational methods, elements from control theory and statistics are combined allowing for future transfer of methods between the two fields.
Keywords:Parameter estimation; Calculus of variations; Boundary value problem; Optimal control; Reaction networks; Ordinary differential equations; Statistical inference
Frequently, signalling pathways and chemical reaction networks in systems biology are modelled by ordinary differential equations (ODE). In many cases, the reaction networks are open systems comprising external inputs like drug stimuli. The system is then modelled by a non-autonomous ODE.
Similarly, modules of reaction networks are open systems. The nodes they have in common with the surrounding network are not or not entirely determined by the module species. They can be considered as intrinsic inputs and again the system can be modelled by a non-autonomous ODE. An example for such a cross-talk can be found in .
While reaction rates and initial reactant concentrations form a countable set of parameters, inputs correspond to an innumerable set of parameters since in general, every function of time is possible as input and each function value at each time point is a free parameter. Commonly, if measurements for the inputs are available, non-parametric estimates like smoothing splines are employed to describe the input data [2,3]. Given the input, an objective function depending on rate parameters and initial values is defined and its minimum is approached by numerical optimization methods. In this way, the problem of infinitely many parameters is avoided. As we will show, one problem associated to this approach is that it does not account for the uncertainty present in the input. As a consequence, estimated parameter confidence intervals do not cover the actual variability, i.e. they are too small.
Therefore, it is preferable to parametrize the input which is possible if certain knowledge about origin and processes underlying inputs is available. This enables a reasonable choice of basis functions and the parametrization becomes finite. Following this approach, the problem of erroneous confidence intervals is circumvented presuming that the input model is correct. However, this assumption is problematic if only sparse information about the inputs and few measurement points are available.
We propose to approach the problem of input parametrization by calculus of variations. In the Result section, the system’s objective function used for ordinary parameter estimation is extended to a functional to be minimised. The original non-autonomous ODE is transformed into an augmented autonomous ODE. The result is interpreted and applied to simulated data.
Results and discussion
Derivation of the augmented ODE system
In conventional parameter estimation, the objective function to be optimised is the likelihood function or the χ2 function. If a reaction network with species yμ, and reaction parameters pk, , comprises inputs xν, , the dynamics of the system is described by the model
with dynamic variables yμ and time-dependent input functions xν(t), each of them collected in the vectors and . In the following, the dependence on the whole course of x will be emphasized by the notation [x]. Furthermore, it is assumed that the input function x(t) is differentiable. Commonly used inputs like step functions or injections are rather distribution like than differentiable functions. However, it is assumed that on the physiological level the acting input is more accurately described by a differentiable function. The χ2 objective function
penalizes distances between species measurements and model prediction yμ(ti, [x], p,y0) at time points ti quadratically and weighted by the measurement uncertainties . In addition, input measurements are compared with the input function values xν(ti). In particular, χ2 is already a functional of [x]. In case of Gaussian noise, eq. (2) coincides with the maximum likelihood estimator.
Our aim is to find a unique input function which minimises the functional defined in eq. (2). To this purpose, we compute the first variation and check under which condition the first variation vanishes. See  for a general introduction to variational calculus as well as sections 1-2 in the Additional file 1. For the first variation we obtain
The trajectory variation δyh is derived by eq. (1) and is expressed by variation of constants: Φ(t) denotes the fundamental system of the homogeneous linear problem with the matrix ∇yf of first derivatives of f with respect to y and ∇xfh constitutes the inhomo- geneity. Furthermore, a weighted residual function is defined as , analogously . For a detailed derivation see sections 2-5 in the Additional file 1.
Next, h needs to be separated. Similarly to Euler-Lagrange’s equation , partial integration needs to be performed to extract h from the integral. However, therefore the sum in eq. (2) needs to be extended to an integral, all time-discrete measurement points and have to be replaced by continuous and differentiable data representations by means of a mapping from N discrete values to a differentiable function. The resulting representations , as well as , need to be defined at least on a finite interval [0, T where T denotes the latest time point to be considered. After partial integration, the first variation for the just defined time-continuous χ2 functional reads
with the auxiliary function u. The transpose is denoted by ∗. The integral, i.e. the first variation, vanishes for all choices of h if and only if the integrand is zero, leading to eq. (5). The auxiliary function u is equivalently expressed by its corresponding differential equation, eq. (6). Here, it is used that Φ−1 is a fundamental system for which follows from Φ being a fundamental system for . Together with eq. (1) we obtain:
The right-hand sides of eqs. (67) depend on state variables y, u, and x, the latter being constrained by eq. (5). Particularly, if the input enters linearly in the dynamics of the reaction network, ∇xf is independent of x and eq. (5) can be directly solved for x, i.e. . However, even in the non-linear case, the implicit function theorem provides the possibility to check locally whether eq. (5) uniquely defines x(uy). For the discussion of a global version of this statement, see section 6 of the Additional file 1.
From the definition of u it follows that u(T) = 0 needs to vanish at the final time point T. Hence, the augmented ODE system needs to satisfy two-way boundary conditions y(0) = y0 and u(T) = 0. This fact constitutes a remarkable difference to the original initial value problem.
Starting from a dynamic system with inputs and measurements for both, state variables and inputs, we have derived differential equations for both of them. The original initial value problem has been transformed into a boundary value problem which is to be solved numerically. The solution trajectories Y(t|p,y0) = (y(t|p,y0),x(t|p,y0)) minimise the χ2 functional for given dynamic parameters p and initial values y0. However, there is still notable freedom in the choice of data and uncertainty representations, denoted by Sy, Sx and Sσ, which decides about the meaning of the solution trajectories.
One possibility to define time-continuous data representations Sy and Sx is smoothing splines. They constitute prior knowledge for each component about shape and time-scale of changes based solely on the measurement points. Also Sσ needs to be chosen appropriately. Differences between model prediction and data prior are usually weighted by at each time point t. Especially if data sampling is sparse, the data prior has larger uncertainty when far away from measurement time points. In this case, a reasonable choice of w(t) is given by
i.e. a sum of Gaussians located around the measurement points. The parameter τ is a measure for the correlation length of the data prior.
Once data and uncertainty representations are chosen, the solution trajectories Y can be employed for conventional parameter estimation minimising
over the finite dimensional parameter space of p and y0. Note that the time-discrete χ2 function and the time-continuous χ2functional do not coincide exactly. Thereby, different measures of optimality are applied to input functions and parameters. This difference is resolved in the asymptotic case of infinitely many measurement points.
The distinction between parameter estimation and input reconstruction has further implications on the estimation of uncertainty bounds. Confidence intervals can only be assigned to the dynamic parameters and initial conditions. In contrast, the input becomes a usual state variable by construction. For state variables, the confidence region in parameter space needs to be mapped to state space by prediction, i.e. by evaluating the model for different parameter values within the confidence region. This can e.g. be realized by parameter sampling using MCMC methods. Alternatively, profile likelihoods can be employed .
It is important for the interpretation of x(t) as a species concentration that x(t) > 0 for all times t ∈[0,T]. This is not imposed a priori on the solution x(t). Rather, it needs to be enforced by construction, analogously to the state variables in the ODE of the dynamic system. This can be realized by the following extension of the dynamic system,
with a diagonal matrix of new inputs . By construction, x can not change sign over time. The choice SD(t) = 0 and for all t reflects a constant input prior with penalized first derivative and can serve as starting point. Besides enforcing positivity of the input, the extension by eq. (11) presents a workaround for dealing with non-linear inputs because the new input variables dν enter linearly and the old inputs xν become regular state variables.
If f depends linearly on x, eq. (5) can be solved for x explicitly. This ensures computational efficiency. In the non-linear case, matrix inversion has to be performed in each evaluation step of the ODE which might slow down the computation of the solution remarkably. Alternatively to the introduction of new input variables, eq. (11), the computationally intensive approach can be avoided by a change of variables. This is possible if state variables and input variables factorize, i.e. if
where gμν and gμ,0 do not depend on the input variables which have been transformed to by a coordinate transformation φ. Examples could be φ(x) = x2 or for a bimolecular or an enzymatic reaction, respectively. The possibility of a change of variables covers a broad range of biologically relevant reaction networks.
Although computation for linear input is remarkably faster than for non-linear input, it is still slower than solving an initial value problem. On the other hand, the solution of the boundary value problem is already optimal with regard to the input course. Therefore, computing time has to be compared to the time a parameter optimization algorithm takes to estimate the parametrized input course. The comparison will strongly depend on the number of parameters that are necessary to parametrize the partially unknown input. So far, there has not been a comprehensive study comparing the two methods.
Application to simulated data
The approach is applied to the following toy model:
The forward reaction is mediated by x while the back reaction is unaffected by the input x. According to eqs. (57), the augmented ODE system for A, B and x is given by
with the auxiliary state variables uA, uB, the data representations SA, SB and the uncertainty representations and . The input x is related to the other state variables by . Several input functions x have been chosen for data generation, among them an exponential decay, x ∼ e−αt, an activation dynamics with a slow decay after a fast increase, x ∼ e−αt − e−βt with α > β, and a Gaussian input, . The example is numerically implemented in C and in R . Optimization is performed by a Gauss-Newton algorithm for nonlinear least-squares estimation.
The purpose of this section is to compare parameter estimation for the variational and the fixed input approach. The input data prior, i.e. the smoothing spline through the simulated input data points, is employed as input function for the fixed input approach.
Examples with Gaussian input are depicted in Figures 1 and 2. All components, A, B, and x depicted in Figure 1A-B have been measured at 20 time points. In this case of dense sampling, the data priors, charted as dotted lines in Figure 1A, come close to the estimated time-courses, charted as dotted lines in Figure 1B. This is reflected in the distributions of the parameter estimates in Figure 1C: for the same set-up, 1000 noise realizations have been generated and the variational approach has been used for parameter estimation. In order to compare the result with the fixed input approach, the data prior of x has been employed as input and conventional parameter estimation has been performed. Hence, in the setting of dense sampling and small noise, both estimation approaches perform equally in terms of accuracy and precision.
Figure 1. Rich input measurement. (A-B) Simulated data for the species A, B and the input x. True time courses are denoted by continuous lines. Data points are subject to Gaussian noise with σ = 0.1. (A) Data representations are indicated as dashed lines, (B) solution trajectories after parameter estimation are shown as dashed lines. (C) Comparison of parameter distributions obtained from 1000 repetitions of data generation and parameter estimation for the variational and the fixed input approach.
Figure 2. Poor input measurement. (A-B) Input reconstruction – simulated data for the species A, B and the input x. True time courses are denoted by continuous lines. Data points are subject to Gaussian noise with σ = 0.1. For the input, only 4 data points are provided. (A) Data representations are indicated as dashed lines, (B) solution trajectories after parameter estimation are shown as dashed lines. (C) Comparison of parameter distributions obtained from 1000 repetitions of data generation and parameter estimation for the variational and the fixed input approach.
A rather different situation is depicted in Figure 2A-B. The input x is measured at four time points only, leading to a poor data prior shown as green dotted line in Figure 2A. Like before, the species A and B have been measured at 20 time points. Most of the information about the dynamics of the input is encoded in these measurements. The correlation time τ has been chosen to be much smaller than the distance between time points allowing for much interstitial variability. The resulting trajectories Y after parameter estimation are shown as dotted lines in Figure 2B. The true input curve is reconstructed almost entirely. The noticeable fluctuations are caused by coincidental noise correlations between species A and B: simultaneous deviations from the true course in opposite directions lead to immediate breakouts of the reconstructed input.
Also for this set-up, 1000 noise realizations have been generated for the comparison of the variational and the fixed input approach. The parameter and initial value distributions for both approaches are shown in Figure 2C. Since the true input can be reconstructed, the variational approach is able to estimate all parameters accurately. In contrast, when the input is fixed to the apparent input data prior, parameter estimation leads to biased parameter estimates.
Finally, we investigated the coverage probability  of the confidence intervals derived from the variational and the fixed input approach: for each simulated data set, parameter estimation is performed, confidence intervals are computed and the information if the true parameter value is situated within the 68%/90% confidence interval is collected. This information is cumulated over many runs of data generation.
Figures 3A and B show the results for Gaussian input with 20 input measurements and 4 input measurements respectively. In each case, 20 measurement points have been provided for each of the species A and B.
Figure 3. Coverage. Coverage for variational approach (blue) and fixed input approach (red). Continuous lines correspond to a 90% coverage probability, dashed lines to 68%. Both probabilities are indicated as black dashed horizontal lines. (A) Example with 20 input measurement points, (B) example with 4 input measurement points.
For both estimation approaches, confidence intervals of estimated parameters and initial values have been produced by means of the profile likelihood approach  with respect to eq. (9).
For the set-up with 20 input measurement points, both estimation approaches provide accurate estimators with similar variances as confirmed by Figure 1. However, as Figure 3 shows, the coverage differs significantly between the two approaches. Confidence intervals for k1 and k2 are systematically underestimated by the fixed input approach. The variational approach in contrast is able to correctly take the degrees of freedom in the input into account. Thus, the coverage is close to the expected values.
For the set-up with 4 input measurement points, the variational approach performs significantly better than the fixed input approach with respect to coverage. However, also the variational approach produces confidence intervals that are slightly too small for the dynamic parameters k1 and k2, Figure 3B left, and too small for the estimated initial values, Figure 3B right. The reason for this behaviour is a combination of the small correlation length τ and the objective function given by eq. (9). Short values of τallow that the input function has fast fluctuations. Especially around the input measurement points, function values tend to punctually approach the measured values, favoured by the time-discrete objective function. Since these fluctuations occur at a short time scale, it has little influence on the course of A and B and thereby, estimation of the dynamic parameters is almost unaffected.
This case shows that τneeds to be chosen appropriately for the problem: small for comprehensive input reconstruction and larger for propagation of uncertainties. A second possibility would be to adapt statistical results for conventional parameter estimation to the case of time-continuous objective functions.
In many applications, it is difficult to guess a proper input model because input data is not available or too noisy. Instead of parametrizing the input, we employed variational calculus to transform the ODE into an augmented system of ODEs describing the original and the input components. The solution of this system minimises the χ2 functional which plays a central role and is directly associated to the objective function of the original estimation problem. Since the extension of the χ2 function to the χ2 functional is not unique, the new functions, i.e. continuous data and uncertainty representations, need to be chosen intentionally. To this end we propose smoothing splines that have a concrete interpretation as data priors. Especially in the case of sparse sampling we propose to use weighting functions for the uncertainties. By this means, existing measurement points are taken into account and the course between time points is not excessively constrained by the data prior.
In the field of control theory and optimal control, so called cost functionals take the role of our χ2 functional. Once chosen the appropriate χ2 functional, our approach to input estimation can be embedded in the general framework of Pontryagin’s minimum principle  and eqs. (67) can be identified as a Hamiltonian system.
We showed that our combined variational approach to parameter estimation enables the assembly of all information present in species and input measurements. By this means, it accounts properly for variability in the input due to measurement uncertainties and produces correct confidence bounds. Depending on the situation, the combination of all information leads to comprehensive reconstruction of the input curves. Information about the dynamics of the input can be concentrated in the species measurements like Figure 2 shows. In such cases our approach clearly outperforms conventional approaches. The variational method is even applicable if no input measurements are available or if species are partially unobserved. A prominent example where the presented method could be applied is the PI3K/AKT/mTOR pathway . Even though various mTOR complexes and their phosphorylated states can be measured, it is not clear how they mediate AKT activation. By applying the variational method to AKT data, it would be possible to reconstruct the required mediator and subsequently relate it to mTOR complex measurements.
A completely different field of application is network modularization. The entire network can be dissected preferably at nodes where measurements are available. These nodes are then treated as independent inputs thus disentangling the network. In this way, the number of equations the variational approach has to deal with is kept small and computational efficiency is ensured.
A further step after the introduction of a time-continuous objective function would be to use the same function for parameter estimation. The time-continuous version of the objective function is closely related to the original function. Therefore, we are confident that it is possible to endow the time-continuous objective function with statistical meaning. This would not only allow for employing the same objective for parameter estimation and input reconstruction in our application. It would also enable the transfer of many more results from control theory and make it suitable for statistical inference.
The authors declare that they have no competing interests.
DK developed the methodology, wrote the software, designed the study and wrote the manuscript. JT supervised the study and critiqued the manuscript. All authors read and approved the final manuscript.
The authors thank Raphael Engesser, Clemens Kreutz and Jan Hasenauer for their advice and valuable discussions. This work has been supported by the German Federal Ministry for Education and Research programme Medical Systems Biology SARA (0315394E).
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.
Bioinformatics 2009, 25(15):1923-1929.PubMed Abstract | Publisher Full Text
Proc Nat Acad Sci 2003, 100(3):1028-1033.PubMed Abstract | Publisher Full Text | PubMed Central Full Text
ArXiv Preprint 2011.
J Am Stat Assoc 2000, 95(450):449-465.Publisher Full Text
Dalle Pezze P, Sonntag AG, Thien A, Prentzell MT, Gdel M, Fischer S, Neumann-Haefelin E, Huber TB, Baumeister R, Shanley DP, Thedieck K: A Dynamic Network Model of mTOR Signaling Reveals TSC-Independent mTORC2 Regulation.
Sci Signaling 2012, 5(217):ra25. Publisher Full Text