Identification of neutral biochemical network models from time series data

Mathematical models in modern molecular biology have become attractive as compactors of the massive amounts of multidimensional data produced by high-throughput techniques, thus following similar ideas that previously led from reductionism to quantitative inroads into physiology and ecology. In the smaller-dimensional world described by the model structure and its parameters, new experiments are easier to conceive, hypotheses can be tested with greater clarity, and knowledge can be extended with inexpensive computational effort [1].

Generally, mathematical models are implemented with a set of parameters, which give them the flexibility of mapping a range of behaviors into a unifying mathematical framework. Except for singular cases where parameters are directly measured experimentally, parameter estimation from experimental data is an inevitable step in the process of constructing models [2]. A good, first-tier compromise between the need for a closed-form, computable representation of the biological process and the risk of ignoring meaningful parameters of mechanistic, hypothesis-driven reductions can be found in the use of generic, "canonical" modeling frameworks. Specifically, within the framework of biochemical systems theory (BST) [3–6], S-system models (Equation 1) offer a particularly convenient solution because their parameters more or less directly describe the interactions between the components of the system of interest [6].

(1)

In the general S-system form (Equation 1), the time variation of the concentration or amount of each component X _i of the system is given by the difference between production and degradation terms. The constant rates α _i and β _i represent the turnover rates of the production and degradation fluxes and the kinetic rates g _ij and h _ij quantitatively characterize the influence of the component X _j on the production and degradation term of the system component X _i , respectively [6, 7]. Thus, the network structure and the nature of the interactions driving the phenomenon under investigation are mapped essentially one-to-one onto the parameter values of the model. This modeling framework has been successfully applied to many biochemical systems [6] and can generally be considered a good first default for representing complex biological systems, especially if the governing mechanisms are not well characterized [8]. Automation of the estimation of the high-dimensional parameter set of an S-system model from multivariate time series data has therefore become a widely pursued computational challenge that has been addressed by a wide variety of optimization techniques: from relatively slow global heuristic optimization techniques like genetic algorithms and simulated annealing [9–12] to fast local optimization algorithms such as alternating regression and eigenvector optimization [13, 14] among others [15, 16]. Most of these optimization algorithms share the strategy of decoupling the differential equation system into a larger, nonlinear algebraic system [6, 11, 17, 18]. This strategy eliminates the need for numerical system integration at each step of the optimization process, which is expensive because S-systems can be numerically stiff, just like most other nonlinear models.

Ironically, the difficulty of finding any numerically integrated S-system model that fits a given set of experimental time series data well is accompanied by the "opposite" problem: many recent publications have pointed out that multiparametric models tend to have the capacity of accommodating whole ranges of parameter values without much affecting the system dynamics [19–28]. Furthermore, it was found that there are typically well-defined directions in the parameter space to which the system dynamics is insensitive, a phenomenon that was termed "sloppiness" [22–25]. Since redundancy appears to be a wide-spread design feature of biological processes, exploring the admissible parameter space of a model and a dataset of interest has relevance that reaches well beyond typical sensitivity analyses of model parameters.

Although the concept of sloppiness has been discussed quite intensely, little attention has been paid to the question of whether or not sloppiness can be translated into the structure of the biological system itself. In other words, is the biological system in reality more or less uniquely parameterized or is there such significant inter-individual variation that we could in principle find large "clouds" of parameter manifestations if we were able to determine the parameters in individual cells or organisms. The answer to this question is not without consequence, because it would affect the definition of what it means to have a good model fit to a given set of data.

In the case of S-systems, the question has further implications. If an admissible parameter cloud, defined by a sufficiently accurate overall fit of some data, permits some kinetic order parameter to be positive, zero, or negative, the interpretation of the estimated model becomes distinctly different. In the first case, the effect is activating, in the second it is negligible, and in the third it is inhibiting. Is such a parameter cloud a computational artifact that would disappear if more data were available, or is it possible that a real biological system would actually allow such different effects of a variable on the system? To reformulate the question, do natural systems only allow for sloppiness in parameter values or also for sloppiness in structure? While the question itself is not new (e.g., [29]), it can presently not be answered in generality and with reliability.

In this report, we address the intrinsic redundancies in the interactions between biological system components by proposing an embedding method for S-system parameter space estimation based on a Monte Carlo process that is combined with an improved version of our previous parameter optimization algorithm. We apply this methodology to experimental time series data characterizing the glycolytic pathway in the bacterium Lactococcus lactis [30–32]. In the same context we explore the concept of sloppiness in S-systems by studying the implications of admissible ensembles of models that dynamically represent the data well but lead to different interpretations.

After preliminary, successful tests with simulated data (see supplementary material), we applied the proposed optimization algorithm to actual time series. These data consist of metabolite profiles from the glycolytic pathway of the bacterium Lactococcus lactis, which were obtained with in vivo NMR experiments [30, 31]. For modeling purposes, the concentrations of the metabolites were coded as follows: glucose – X₁; glucose 6-phosphate (G6P) – X₂; fructose 1, 6-bisphosphate (FBP) – X₃; phosphoenolpyruvate (PEP) – X₄; lactate – X₅; acetate – X₆. As a pre-processing step for the parameter optimization, the time series were smoothed/denoised and their slopes (numerical first derivatives) were estimated as shown in previous work on non-stationary noise filtering [36].

Initially, the proposed algorithm exclusively used the time series of all metabolites in the parameter optimization, representing the case where no knowledge about the network connectivity is at hand. The optimal parameter set of each metabolite was separately translated into the eigenspace of the solution and subsequently fed into the MC process. The parameters were then perturbed and an ensemble of models was created as described in the Methods section. The parameter sets were selected based on a behavioral class [37], which was defined by a residual less than 5 (see Equation 15). As an example, Figure 1 shows all parameter sets for the G6P production term that fall within the defined behavioral class. After ordering, these parameters have the same distribution pattern (Figure 1B). Similar results were found for all other metabolites.

As an alternative, we performed a MC random walk in the original parameter space, as opposed to the eigenspace. Interestingly, although not surprising in the end, the results did not present as wide a range as was found in the exploration of the parameter set in the eigenspace. The difference in outcomes is explained by the compensation of errors between production and degradation terms: perturbations in the eigenspace of the matrix W affect both production and degradation terms while perturbing only one parameter does not always maintain the balance between the two terms. In other words, given one of the S-system's terms (production or degradation) and the vector of slopes, the complementary term can be obtained by multiple linear regression, which has been shown not to be sloppy [24]).

As expected, eigenvalues of the Hessian matrix fall within a sparse range [22–25, 38], thereby elucidating the stiff and sloppy directions (see Additional file1). The region of the parameter space that produces similar dynamical model behaviors can be approximated as an ellipsoid whose main axes are given by the direction of the eigenvectors of the Hessian matrix (see Figure 2) and whose width is inversely proportional to the squared root of the corresponding eigenvalue [22–25, 38]. A projection of the ellipsoid into three-dimensional space for the acetate production parameters g₆₄, g₆₅ and g₆₆ is shown in Figure 2. Revisiting reasons discussed for other modeling frameworks [22–25, 27], sloppiness can be explained in the proposed approach by the neutral space of solutions for Equation 6 and consequently Equation 8. Specifically, given the "right" linear combination of the eigenspace of the matrix W, any vector resulting from stretching or shrinking this combination is also a solution of Equation 6, however with different parameter values.

Ideally, all combinations of parameter sets found for the individual metabolites would generate recoupled models that fit the data upon numerical integration, within some error bound. Given the large number of parameters (at the order of 10¹⁴ combinations), a full exploration of this statement is not possible. In order to ameliorate the expensive combinatorial issue of assessing the trajectories of the recoupled systems, the parameter set for each metabolite generated from the MC process was randomly selected and the resulting system numerically integrated. This process was repeated 500 times and the results are presented in Figure 3.

As can be seen from Figure 3, the uncertainties and bounds in the prediction for most metabolites are actually close to the observed data. A notable exception is lactate, where the measured data contain little noise. Because PEP is the main precursor of lactate in the model [31] (pyruvate is not explicitly modeled), the class of predictions of the PEP dynamics was re-sampled for residuals less than 1. Furthermore, the newly sampled systems were integrated with glucose supplied in three different initial concentrations, namely 20, 30 and 40 mM. The results of these predictions are shown in Figure 4. For all three different scenarios the new ensemble model predictions provide accurate descriptions of the observed dynamics for the concentrations of lactate and the other metabolites in the model, except for the noticeable undershooting of PEP (see [30, 31] for information on the responses of the system to different initial glucose concentrations).

Although the models found by the proposed procedure accurately describe the dynamics observed in the experimental data, none of the parameter sets match well with the network topologies found in the literature, suggesting that distinctly different parameter combinations, and not just sloppy versions of some parameter set, area able to match these data. A variety of techniques can be applied to the ensemble in order to compare different models [37, 39], to cluster models by means of transformation groups, or even reduce them to a smaller subset [20]. One possible avenue for further analysis of this vast parameter space is by creating groups defined by different behavior classes based on biological and dynamical information [37, 39].

Regardless of specific follow-up analyses, one of the purposes of this work is to demonstrate that the proposed optimization algorithm can be effectively applied to the parameter identification of specific networks, for instance, by taking kinetic parameters of some component X _j out of the optimization process of X _i . Conversely, previous knowledge can be used to restrict the values of the g and h interaction parameters. For example, existing knowledge about Lactococcus lactis primary metabolism [30–32, 40] provides precise clues about which interactions are reasonable. To analyze the benefit of information outside the time series data, we applied the proposed optimization algorithm to the same Lactococcus lactis data, this time however constraining some parameter values with the pre-existing information. The simulation results with the resulting system (Figure 5; see equations in Additional file1) not only describe the data but also agree with the double pulse of glucose described in [32].

Exactly as described before, the system resulting from this combined approach could be used as the starting point of a MC random walk process, and the same analysis of sloppiness and behavior classes could be performed. In order to assess the accuracy of the solutions for large systems, the optimization algorithm was also tested with a symbolic genetic network model consisting of 30 components [41] (see Additional file1). Because the algorithm performs the parameter optimization with the decoupled form of the system, its complexity is linear with the number of the system's components [14]. Thus, rather than system dimension, the real limitation of the optimization algorithm is the time series dynamics. Poor dynamical variability (components close to the steady-state) and collinear time series will result in a conditioning deficiency of the matrix L, causing numerical problems with the inverse operations and misleading the convergence. Despite the successful retrieval of the 30-dimensional system dynamics, the algorithm's limitation becomes more evident with large systems (higher possibility of having time series portions with collinear components, resulting in ill-conditioned blocks in the matrix L [42]). This drawback is partially resolved with the regularization technique presented in the Methods section. Also, problems of this nature can be prevented by removing the collinear components from the matrix L or placing them as independent variables (variables that do not take place on the parameter optimization [6]). Of course, this issue is drastically diminished when a chosen network topology constraints the matrix L.

Within the new field of systems biology, an extraordinary effort has been devoted to kinetic modeling with the aim of understanding biological processes better [1, 43, 44]. A very significant part of this effort has been directly related to the identification of model parameters [2, 45–47]. Until recently, the quality of the estimated parameter values was judged by the fit of experimental data. However, new analyses have pointed out the importance of other criteria, such as extrapolability and error compensation among terms [48], which is in some sense related to the sloppiness of admissible parameter sets and the fact that different parameter values can generate similar dynamical behaviors in nonlinear biological systems models [22, 23, 25–27]. These observations have a direct impact on the robustness of models, which may translate into robustness of the biological system itself, which may become apparent in tolerance to mutations, changes in gene expression, and insensitivity to modest changes in environmental conditions. This need for physiological robustness implies that biochemical networks should be able to preserve their dynamical properties within reasonable ranges of their kinetic parameters [49].

In this report, we present an extension of our previous optimization algorithm for S-system parameter identification from time series data [14]. The proposed method turned out to be faster and more accurate than its predecessor (see Additional file1) and was used here in combination with a Monte Carlo random walk technique to explore the space of admissible parameter sets of S-system models. This strategy allowed us to explore the concept of sloppy models. The results indicate that both, a fully integrated and a decoupled model, can be sloppy. We also reanalyzed time series of the concentrations of six metabolites within the glycolytic pathway of the bacterium Lactococcus lactis and demonstrated how a sloppiness analysis can elucidate the admissible parameter space and ultimately lead to more reliable estimates.

The central result reported here is that a diversity of parameter sets may produce quasi-isomorphic dynamics for S-system models. Most of the parameter variations extended to both positive and negative parts of the parameter space (Figure 1B). This result is interesting, because it could be the consequence of two distinctly different scenarios. First, it could reflect redundancy or sloppiness caused by insufficient data. In other words, the data are not informative enough to distinguish between alternative models that fit equally well. In the past, such situations have often been "resolved" by setting as many parameters to zero as possible with the set of admissible solutions, borrowing arguments of parsimony or Ockham's razor. The important feature of this scenario is that further experimental data, maybe obtained under similar yet sufficiently different conditions, would declare one of the candidate models the (sole) winner. The second possible root cause of distinct, well-fitting parameter sets is the actual natural co-existence of different regulatory signals between components (inhibition versus activation) or different regulatory networks (some with and some without particular regulatory signals). Thus, otherwise similar cells or organisms would function under slightly different regulatory regimens. The difference between these scenarios is conceptually similar to the distinction between uncertainty and variability, which was discussed intensely in the 1980s and 1990s within the fields of risk assessment and exposure analysis [50]. The former case of uncertainty (due to insufficient data) is only valid if the model is, in principle, uniquely identifiable and structurally distinguishable [19, 51] from among all feasible S-systems. This aspect raises the possibility that a sloppy model could be unidentifiable and that the sloppy directions, which are given by the possible parameter combinations [22] could be a measure of non-identifiability [51]. For S-systems this argumentation can be extended. Our results show that an ensemble of S-systems could be interpreted as a collection of structurally indistinguishable and unidentifiable models [19]. This conjecture would explain the range of variability (negative, zero, positive) of the parameter sets that were found. Furthermore, the identifiability characteristic of the mathematical framework could be associated with the robustness of the biological system to environmental changes. Analogously, the structure distinguishability characteristic could be, in some sense, associated with the robustness of the biological system to mutations that change the network of interaction among components. A more rigorous study of S-system identifiability and distinguishability will be needed to reveal more concrete conclusions and mathematical implications [52, 53].

Whether the distinction between parameter sloppiness and structural distinguishability in biological models is an important issue will have to await further investigation. Nevertheless, it is clear that extrapolability in the generic network identification problem from time series data is a more complex task than the computational fitting of a model. One should expect that this observation is true for any model structure, but it was shown here that S-system models allow its exploration in a most translucent fashion.

The assessment of sloppiness provides valuable information that can be extracted from the prediction of the ensemble of models and through the investigation of behavioral classes that differ in their dynamical features (e.g. overshoots or response time). These classes may be further reduced by biological knowledge such as biochemical, physico-chemical, or thermodynamic constraints [54, 55]. The proposed techniques could also serve as a powerful exploration tool for the testing of hypotheses and the design of new experiments. Moreover, the union of the proposed optimization algorithm with statistical methods could also result in a robust network inference implementation [56].

Maybe the most important value of this report is the clear definition of a framework to explore sloppiness, robustness and evolutionary innovation [38, 49, 57], where neutral parameter spaces (system with essentially the same dynamical behavior) are merged with neutral networks [49, 57] into one unique structure. This approach could reveal insights into how different metabolic networks could possibly have been changed during the evolutionary process.

The method described here is an extension of our previous work in S-systems parameter estimation [14]. Exactly as its predecessor, the proposed algorithm exploits the advantages of the numerical decoupling of S-systems, which allows that each component has its own optimization process [6, 11]. However, the optimization is now performed in a transformed equation of the decoupled form. This transformation enhances the ability of the algorithm to reach local optimal points. Given a time series collection X _m (t _i ) of the concentrations of species or metabolites m = [1,..., M] in the time interval t _i = [t₁,..., t _N ] and their respective slope values S _m (t _i ) estimated in the same time interval, one can numerically decouple the S-system equations and write a semi-linearization (only the left-hand side is linearized) form for each species m as follows

(2)

where

(3)

and Vp _m = [log α _m g_{m 1}g_{m 2}⋯ g _mM ] is the parameter vector for the production term. In Equation (2), the time point index t _i is suppressed in order to keep a simple notation. If a multiple regression step is applied to the system in Equation (2), the production parameter vector can be written as

(4)

The matrix L⁺ is the Moore-Penrose inverse of the matrix L, given by L⁺ = (L ^T L)^-1L ^T . Substituting Equation (4) into Equation (2), we obtain the following eigenvector problem [14, 58]:

(5)

In this equation, the matrix W = LL⁺ is idempotent with M+1 unit eigenvalues, while the remaining eigenvalues are zero. It is clear from Equation (5) that not only but also are eigenvectors of W corresponding to eigenvalue 1. If Equations (2), (4) and (5) are rewritten in order to isolate the degradation parameter vector instead of the production parameter vector, the vector also belongs to the eigenspace. Standard routines for eigenvalue analysis can easily calculate the eigenvectors of matrix W. However, this does not imply that any of these eigenvectors will map onto the "correct" parameter set. In fact, these vectors form the eigenspace of W corresponding to eigenvalue 1, leading to the conclusion that the true parameter set can be formulated as a linear combination of these eigenvectors. Thus, let EigS be a matrix where the columns are the eigenvectors with correspondent eigenvalues 1. Then the following vectors can be defined:

(6)

or

(7)

In Equations (6) and (7), ψ and δ are (M+ 1)-dimensional vectors that represent an arbitrary linear combination of the eigenvectors of W. We can use Equation (7) to write the decoupled form of the S-system as

(8)

Equation (8), once written as a function of the production and degradation parameter vectors, can now be seen as a function of the vectors ψ and δ. Using the uniqueness of the mapping between the vectors ψ and δ,

(9)

we can rewrite Equaiton (8) as a function of only the vector δ. For the estimation of vector δ and consequently of vector ψ, we define a cost function as the sum of squared residuals between the two sides of Eq. (8).

(10)

where is the right-hand side of the Equation (8). Thus, Equation (10) leads to the following gradient equations

(11)

or

(12)

where

(13)

The symbols and represent the Hadamard and Shanghai Jiao Tong (SJT) product respectively [42, 59], defined as

(14)

In Equations (11) and (12), υ is the argument of the logarithm in Equation (10). We use the Levenberg-Marquardt routine [60] for the minimization of F with the following nonlinear constraints

(14a)

In these inequality constraints (see also Equation (7)), the pairs lcb-ucb and lkb-ukb are the intervals (lower and upper bounds) for the rate constants and kinetic orders respectively. Preliminary tests with a simulated systems show the potential and accuracy of the algorithm (see Additional file1). Without change in the gradient equations, a Tikhonov regularization [61] was also implemented in order to overcome potential problems with the ill-conditioning of the matrix L, resulting in the following matrix:

(15)

where λ is the regularization parameter and I the identity matrix.

Decoupled ensembles of S-system models

Recently, a Monte Carlo random walk method [33] in the parameter space of multiparametric system biology models was suggested in order to assess uncertainties in model parameters and model structure [22, 23, 25, 26, 38]. In this context, S-system models are particularly convenient because they map any network topology directly onto its parameter values (g and h) with no change in model structure [6]. Therefore, we propose a walk in the parameter space of the S-system and proffer that it corresponds to a walk in the 'topology space'. This idea is explored here by combining of the Monte Carlo approach with a nonlinear optimization algorithm. In contrast to relevant work in the literature, but in accordance with our optimization method, the random walk is performed separately with the decoupled S-system format for each system component. This strategy avoids difficulties encountered in the numerical integration of stiff equations and poses no significant loss of information, under the reasonable assumption that there are no bifurcation points within the considered parameter ranges. Thus, we define the following cost function, which quantifies variations in the first derivatives of a variable with respect to variations in the parameter set Δδ:

(16)

Here, δ* is the optimal parameter set found using the proposed optimization algorithm. Evaluation of the Hessian matrix,

(17)

allows us to analyze the sensitivity of the first derivative of each component of the system in relation to its parameters. This Hessian is also used to guide the Monte Carlo random walk, which selects the next perturbations within the parameter set. Thus, the Monte Carlo random walk starts with the optimal vector δ* as initial condition, randomly selects one component of this vector and perturbs it proportional to the Hessian eigenvector (see [25, 35]), and uses the new vector as initial condition for the proposed optimization algorithm, which now has an early stop criterion (number of iterations). The new optimal vector Δδ (i.e., the outcome of the optimization algorithm) is accepted as the next step of the process with probability

(18)

[26], where

(19)

is the probability distribution of the parameter set δ given the model for the slopes δ m(δ*). In Equation (19), k is a normalization factor (that vanishes in the probability of acceptance [26], Eq. (18)) and σ is the standard deviation of the slopes S _m . If a new vector δ is accepted for the next step, the hessian matrix is recalculated and the process follows as described above. This strategy permits a more detailed exploration of the parameter space where the MC steps are taken in the sloppy directions. A small range for the number of iterations of the optimization algorithm (early stop criteria) was tested and no significant differences were observed. Relative large numbers were avoided in order to prevent the algorithm convergence to the same neighborhood. After the optimization process, the vectors δ are mapped onto the parameter degradation vector Vd _m = [log β _m h_{m 1}h_{m 2}⋯ h _mM ] ^T and consequently onto the parameter production term Vp _m = [log α _m g_{m 1}g_{m 2}⋯ g _mM ] ^T , restating the decoupled system in its original parameters. All the numerical integrations presented in this report were carried out by the MATLAB^® ode23s routine. The rate constants were optimized within the range [0.1, 300] while the kinetic order were optimized within the range [-2, 2]. All the results shown in this report can be reproduced with the freely available MATLAB^® scripts (see Additional file2).