Parameter estimation in modulated, unbranched reaction chains within biochemical systems

Abstract

Modern biology is increasingly developing techniques for measuring time series of global gene expression and of many simultaneous proteins or metabolites. These data contain valuable information on the dynamics of cells, which has to be extracted with computational means. Given a suitable mathematical model, this extraction is in principle a straightforward regression task, but the complexity and nonlinearity of the differential equations that describe biological systems cause severe difficulties when the systems are of realistic size. We propose a method of stepwise regression that can be applied effectively to linear portions of pathways. The method may be combined with other estimation methods and either directly yields reasonable parameter estimates or at least provides appropriate start values for subsequent nonlinear search algorithms. We illustrate the method with the analysis of in vivo NMR data describing the dynamics of glycolytic metabolites in Lactococcus lactis.

Introduction

Mathematical models of biochemical systems have traditionally been constructed from the bottom up. Individual enzyme-catalyzed steps were formulated as rate laws and their parameters were estimated from experimental results in vitro, for instance, in the form of Lineweaver-Burk plots (e.g., Segel, 1991). Subsequently all rate functions were integrated into a comprehensive model describing the dynamics of select metabolite pools of interest. The integrative behavior was compared to observed system responses, if these were available, which was rarely the case, and discrepancies were resolved through secondary adjustments of parameter values.

Modern biology is generating new classes of data that permit a complementary approach. These data are of high density, capturing simultaneously the expression of hundreds or thousands of genes or proteins, or the concentrations of many metabolites. As a direct extension, the same experimental techniques allow the construction of time series of these measurements. These time series are very interesting, because they illuminate the dynamics of cellular responses in unprecedented ways. It is easy to imagine that these data have the potential of offering truly novel insights into the functioning of biological systems, if efficacious methods are available for their analysis. Indeed, several groups around the world have begun establishing methods that estimate model parameters from high-density time series (e.g., Voit and Savageau, 1982, Voit and Sands, 1996, Voit, 2000, Maki et al., 2002, Almeida and Voit, 2003, Kikuchi et al., 2003, Voit and Almeida, 2003, Voit and Almeida, 2004, Veflingstad et al., 2004, Kimura et al., 2005, Tsai and Wang, 2005, Naval et al., 2005; see also Vance et al., 2002, Torralba et al., 2003).

The estimation of system parameters from time series data is in principle a straightforward regression task, which, however, is obstructed by severe challenges in implementation. One class of issues pertains to the choice of the mathematical model that is supposed to capture the observed dynamics. Clearly, there are no a priori rules or laws that identify one specific type of model as optimal for a given biological system. Thus, one has to make assumptions and either use functions that have worked well in the past, or one employs generic representations that are supported by mathematical theory. The former approach might use modulated Michaelis-Menten or Hill functions, while the latter might use power-law representations. The second class of issues is of a technical nature and, for instance, includes the lack of convergence of nonlinear estimation algorithms, which is aggravated by the combinatorial explosion in the number of parameters to be estimated if the system grows in size.

In this article, we discuss a technique that, in combination with existing methods, helps tame the combinatorial explosion in the size of the parameter search space. This technique deals efficiently with linear portions of complex pathways, which may or may not be regulated by activators or inhibitors. The method does not provide advantages over other methods at branch points, but we show how linear parts may be separated out and thus render the proposed techniques admissible. This is important, because many parts of biochemical systems are indeed linear. As just one example, the conversion of PRPP into IMP consists of an unbranched chain of 10 steps, the first of which is subject to several inhibitory signals (cf. Stanbury et al., 1983: Chapter 50).

The method begins with the decoupling of the system equations, as it was proposed in the recent literature (Voit and Almeida, 2004; see also Voit and Savageau, 1982, Volgin et al., 2003, Clements et al., 2004); we will review this method in the next section. The decoupling renders it possible that the parameters of the entire system be estimated one equation at a time. Specifically, one considers the measured data quasi as forcing functions for the equation presently under consideration. Upon decoupling, and using power-law representations for the underlying model, we show that linear segments of pathways may be estimated with linear regression techniques that are extremely efficient.

We illustrate the method with time series data that Neves et al. (1999) obtained from Lactococcus lactis with methods of in vivo nuclear magnetic resonance.