Tight Interval Bounds for the Parametric Solutions of Sytems Biology Models

Biochemical reaction kinetics models and compartmental models formulated as systems of ordinary differential equations (ODEs) often do not have a known analytical solution. In the absence of an analytical solution, it is difficult to explore the behavior of the model solution with respect to variations in parameters, such as rate constants and initial conditions. Unfortunately, rate constants are often not known accurately, especially for complex biological systems where the influence of individual model parameters may be very difficult to discern experimentally. While sensitivity analysis provides local information about the parametric dependence of the model solution, it has limited applicability to the global behavior of the solution subject to a range of possible parameters and/or initial conditions. The method presented here provides time-varying upper and lower bounds on the solutions of systems biology models over a specified range of parameter values and /or initial conditions. Such bounds are useful for propagating uncertainty or measurement errors in parameters and initial conditions through a model, and for globally solving dynamic optimization problems involving systems biology models, which includes the important case of parameter estimation problems.

Many authors have addressed the problem of generating time-varying bounds on the parametric solutions of general ODEs in the literature on dynamic optimization, process safety verification and validated numerical integration. Many of these techniques are only effective for a very restrictive class of ODEs, which excludes all but the simplest systems biology models, and result in bounds which dramatically overestimate the parametric solution for ODEs outside of this class. The remaining techniques are based on Taylor expansions and can potentially provide tighter bounds, but often only through the use of high-order expansions which are very computationally expensive. None of these techniques are intended for systems biology models specifically, and therefore do not take advantage of the special structure of these models. In particular, the method proposed here identifies and exploits affine invariants such as reaction invariants and metabolic pools, which are known to exist for ODE models which can be written in terms of a matrix pre-multiplying a vector of rates, as is the case for biochemical reaction kinetics models and compartmental models. These affine invariants may correspond to the conservation of mass in closed subsystems, or pools formed by the reaction network topology, or to stoichiometric proportionality relations between species, and are particularly prevalent in systems biology models. For example, an enzyme which may exist in one of many excited or complexed states is often represented by several distinct species in a model, but the sum of the molecules of these species is invariant as the reactions proceed. Accordingly, such species need not be treated as entirely independent when generating bounds.

The bounding technique proposed here is an extension of the work in [1], where bounds are derived for the solutions of general ODEs based on an application of a standard results from the study of differential inequalities along with basic interval arithmetic methods. This technique is computationally inexpensive but may produce very weak bounds in the general case.  However, for systems biology models which obey affine invariants, it is not difficult to see that the solutions of the model must lie in an affine subspace of the full space of species concentrations.  Until now, no method has been proposed for enforcing these affine constraints when using bounding techniques involving interval arithmetic.  In the method proposed here, this complication is circumvented by using the invariants to define an affine mapping from the space of species concentrations to a lower-dimensional space, and to define an equivalent system of ODEs there.  In this reduced space, established bounding techniques [1] are applied directly.  Bounds on the systems biology model solutions can then be calculated through a suitable inverse mapping.  This method typically produces much tighter bounds than those computed without using the invariants, while the additional computational cost is negligible.

Using this technique, it is possible to produce accurate, meaningful bounds on the solutions of fairly large systems biology models over substantial ranges of parameters and initial conditions. At present, bounding techniques in the literature are only capable of accurately bounding ODEs with three or four state variables and fewer uncertain parameters, and often results for these systems are only presented over very narrow parameter ranges. Accordingly, these methods are of little use for large biological networks where many rate constants are not accurately known. In contrast, preliminary results using the method proposed here show very accurate bounds for models with roughly fifteen state variables and ten uncertain parameters varying over orders of magnitude. Thus, this technology can provide tight bounds on the parametric solutions of biologically relevant models under significant parametric uncertainty. Furthermore, this method potentially facilitates the global solution of parameter estimation problems for many such systems.

REFERENCES

1. A. B. Singer, P. I. Barton, Bounding the solutions of parameter dependent nonlinear ordinary differential equations," SIAM J. Sci. Comput., vol. 27, pp.2167-2182, 2006.