A variational approach to parameter estimation in ordinary differential equations

Background 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. Results 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. Conclusions 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.


Background
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 nonautonomous 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 nonautonomous ODE. An example for such a cross-talk can be found in [1].
While reaction rates and initial reactant concentrations form a countable set of parameters, inputs correspond *Correspondence: daniel.kaschek@physik.uni-freiburg.de 1 Institute of Physics, Freiburg University, Freiburg, Germany Full list of author information is available at the end of the article 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 http://www.biomedcentral.com/1752-0509/ 6/99 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.

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 μ , μ = 1, . . . , n and reaction parameters p k , k = 1, . . . , r, comprises inputs x ν , ν = 1, . . . , m, 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 y ∈ R n and x ∈ R m . 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 . In addition, input measurements x D ν,i are compared with the input function values x ν (t i ). 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 [4] 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 probleṁ φ = ∇ y f φ with the matrix ∇ y f of first derivatives of f with respect to y and ∇ x fh constitutes the inhomogeneity. Furthermore, a weighted residual function is , analogously res x ν . 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 [4], 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 y D μ,i and x D ν,i have to be replaced by continuous and differentiable data representations by means of a mapping S : R N → C 1 (R) from N discrete values to a differentiable function. The resulting representations S y μ (t), S x ν (t) as well as S σ yμ (t), S σ xν (t) 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φ = −∇ y f φ which follows from being a fundamental system foṙ φ = +∇ y f φ. Together with eq. (1) we obtain: The right-hand sides of eqs. (6-7) 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, ∇ x f is independent of x and eq. (5) http://www.biomedcentral.com/1752-0509/6/99 can be directly solved for x, i.e. x = S x − S σ 2 x ∇ x f * u. However, even in the non-linear case, the implicit function theorem provides the possibility to check locally whether eq. (5) uniquely defines x(u, y). 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) = y 0 and u(T) = 0. This fact constitutes a remarkable difference to the original initial value problem.

Interpretation
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, y 0 ) = (y(t|p, y 0 ), x(t|p, y 0 )) minimise the χ 2 functional for given dynamic parameters p and initial values y 0 . However, there is still notable freedom in the choice of data and uncertainty representations, denoted by S y , S x and S σ , which decides about the meaning of the solution trajectories.
One possibility to define time-continuous data representations S y and S x 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 w(t) = 1 S σ 2 (t) 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 y 0 . Note that the time-discrete χ 2 function and the time-continuous χ 2 functional 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 [5].

Technical remarks
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, where g μν and g μ,0 do not depend on the input variables which have been transformed tox = ϕ(x, p) by a coordinate transformation ϕ. Examples could be ϕ(x) = x 2 or ϕ(x, K D ) = x K D +x for a bimolecular or an enzymatic reaction, respectively. The possibility of a change of variables covers a broad range of biologically relevant reaction networks. http://www.biomedcentral.com/1752-0509/6/99 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 A → B is mediated by x while the back reaction B → A is unaffected by the input x. According to eqs. (5-7), the augmented ODE system for A, B and x is given bẏ with the auxiliary state variables u A , u B , the data representations S A , S B and the uncertainty representations S σ 2 A and S σ 2

B
. The input x is related to the other state variables by x = S x + S σ 2 x k 1 A(u A − u B ). 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 [6]. 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.
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 [7] 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.
For both estimation approaches, confidence intervals of estimated parameters and initial values have been produced by means of the profile likelihood approach [8] 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 http://www.biomedcentral.com/1752-0509/6/99 the two approaches. Confidence intervals for k 1 and k 2 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 k 1 and k 2 , 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 timecontinuous objective functions.

Conclusion
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 [9] and eqs. (6-7) 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 http://www.biomedcentral.com/1752-0509/6/99 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 [10]. 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 timecontinuous 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.