(173a) 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.