Bioinformatics Advance Access originally published online on November 29, 2005
Bioinformatics 2006 22(3):346-353; doi:10.1093/bioinformatics/bti800
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Conservation analysis of large biochemical networks
Keck Graduate Institute 535 Watson Drive, Claremont, CA 91711, USA
*To whom correspondence should be addressed.
 |
ABSTRACT |
|---|
TOP
ABSTRACT
1 INTRODUCTIONMotivation: Large biochemical networks pose a unique challenge from the point of view of evaluating conservation laws. The computational problem in most cases exceeds the capability of available software tools, often resulting in inaccurate computation of the number and form of conserved cycles. Such errors have profound effects on subsequent calculations, particularly in the evaluation of the Jacobian which is a critical quantity in many other calculations. The goal of this paper is to outline a new algorithm that is computationally efficient and robust at extracting the correct conservation laws for very large biochemical networks.
Results: We show that our algorithm can perform the conservation analysis of large biochemical networks, and can evaluate the correct conserved cycles when compared with other similar software tools. Biochemical simulators such as Jarnac and COPASI are successful at extracting only a subset of the conservation laws that our algorithm can. This is illustrated with examples for some large networks which show the advantages of our method.
Availability: The software is available as part of the latest release of Systems Biology Workbench (SBW version 2.5.0) and can be downloaded from http://www.sys-bio.org. The software is licensed under the BSD open source license and is freely available at sourceforge.
Contact: rrao{at}kgi.edu
|
1 INTRODUCTION |
|---|
One of the characteristics of biological networks is the conservation of certain molecular subgroups, termed moieties (Reich and Selkov, 1981). A typical example of a conserved group is the conservation of the adenine nucleotide moiety, i.e. the total amount of ATP, ADP and AMP is constant during the evolution of the system. Other more common examples include the conservation of protein between phosphorylated and unphosphorylated states. Over long time periods the conservations do not hold as other slower processes will intervene. However, in many studies conservation laws represent important invariants which can greatly influence dynamical behavior (Sauro and Ingalls, 2004; Koshland et al., 1982).
In addition to effects on behavior, conservation laws in a model can affect the way we carry out certain numerical procedures. The Jacobian matrix in particular is a central quantity that is required in the computation of many analyzes. Examples include, Bifurcation Analysis (Chickarmane et al., 2005), Frequency Analysis, Metabolic Control Analysis (Reder, 1988) and more traditional approaches such as the numerical integration of stiff differential equations, steady-state analysis and certain optimization approaches. All these analyzes depend on the computation of a non-singular Jacobian matrix. When a model contains conserved moieties, the Jacobian becomes singular and such analyzes are no longer possible. Evaluating the conservation laws is therefore critical. Finally, taking into account conservation laws also permits a simulator to reduce the size of a model and thereby increase the performance of the software.
The determination of conservation laws in biochemical network models has been a subject of discussion for many years. The interest in this topic can be traced back to the study of stoichiometric networks by the chemical engineering community in the 1960s (for a review, see Sauro and Ingalls, 2004). Studies of relevance include works by Aris (1965) to be followed by Horn and Jackson (1972), Clarke (1980) and Feinberg (1989). Many of the early simulators, such as SCAMP (Sauro, 1993) and Gepasi (Mendes, 1993) specifically incorporated algorithms to compute the conservation laws both as a means to reduce the model size and more importantly as a means to compute the reduced Jacobian. None of the approaches taken by these early simulators were robust for large networks, indeed as we will show, a number of contemporary simulators suffer the same problem.1
Traditionally, most simulators have implemented Gaussian Elimination (Hofmeyr, 1986; Holstein and Greenshaw, 1994; Heinrich and Schuster, 1996; Cornish-Bowden and Hofmeyr, 2002; Sauro and Ingalls, 2004) to extract the conservation laws. Gaussian Elimination is an easy and fast algorithm to implement and is well suited for small networks. However, it is not applicable to large networks as errors tend to accumulate in the computation of pivots and the subsequent elimination of rows. An alternative to avoid round-off errors involves the use of integer arithmetic. This approach, while being accurate can be computationally slow for large systems.
The development of a robust method of computing the conservation laws for biochemical networks is important for large systems. It is in this regard that the US Department of Energy sponsored GTL (see http://doegenomestolife.org/) program has invested in the development of robust means of handling large biochemical systems. This is a very critical step for the GTL program that envisages among its goals the development of advanced computational methods and capabilities to further the understanding of complex biological systems and predict their behavior. This paper therefore addresses the problem of reliable computation of conservation laws for large networks by presenting a new method based on Householder QR factorization (Householder, 1958). We illustrate this method with a simple example and compare its ability to extract the correct conservation laws for large systems with other similar tools. Various applications of this method are also addressed.
|
2 CONSERVATION ANALYSIS |
|---|
In the following, we adopt the nomenclature pertaining to biochemical networks established by Reder (1988). A biochemical network involves reactions between various species whose time evolution can be described by
 | (1) |
Sm(t)]T is the vector of time-dependent species concentrations, N is the Stoichiometric matrix relating the species to the reactions they participate in and v(t) = [v1(t)v2(t) vn(t)]T is the vector of rates of the reactions that comprise the network. Consequently N is a matrix with m rows and n columns, where each column corresponds to a particular reaction and each row corresponds to a particular species. Equation (1) therefore describes the dynamics of the network species in relation to the reactions. In many networks however, some of the reactions are formed such that certain conservation relations naturally follow. These relations can be interpreted as dependencies among species or alternatively, as dependencies in the rows of the N matrix. This means the rank of the N matrix is less than what it would be without the conservation relations. Let us denote this rank by m0. The goal of any algorithm that identifies these conservation relations should be to identify and separate independent and dependent rows of the stoichiometry matrix, or equivalently partition the species participating in the biochemical network into dependent and independent entities. It follows that there are m0 independent species and m – m0 dependent species. Let us denote these as Si(t) and Sd(t). The time evolution of the biochemical network described in Equation (1) can then be expressed in terms of dependent and independent species as
 | (2) |
 | (3) |
 | (4) |
Further simplification of Equation (4) yields the relation between time evolution of dependent and independent species as
 | (5) |
It follows upon integrating Equation (5) that Sd and Si are related by a constant vector 
such that
 | (6) |
Introducing a matrix 
in place of [–L0 I] in Equation (6), we can deduce that S = T. The (m – m0) x n matrix is called the conservation matrix, as it relates the species vector S to the vector of conserved moieties T. These represent molecular subgroups that are conserved during the evolution of the network (Reich and Selkov, 1981). Each row of the conservation matrix represents a distinct conserved cycle of the network. The values of the vector T can in practice be obtained by substituting the initial conditions of the species into the relation T = Sd(0) – L0Si(0).
Another aspect of interest in the study of biochemical networks relates to the evaluation of the Jacobian Matrix. The dynamics of the network described by Equation (1) are most often non-linear, and hence a linearized system can be obtained by studying the perturbation around a steady-state value. Ignoring second-order terms, a linear equivalent of Equation (1) can be obtained as
 | (7) |
S represents infinitesimal changes in S. The dependencies in the network between species are reflected in the Jacobian matrix JF being singular. The full Jacobian matrix can therefore be built using the relationship between dependent and independent species, given by L and NR matrices. This can be shown after modifying an expression given by Heinrich and Schuster (1996, p. 40) to be
 | (8) |
is the Elasticity coefficient matrix, with ij = 
vi/Sj. A non-singular reduced Jacobian matrix, JR, can be constructed using conservation analysis. Indeed, it can be shown, using Equations (2) and (7), that JR is given by the relation JR = NR L. |
3 DEALING WITH LARGE SYSTEMS |
|---|
The efficacy of a number of methods with varying degrees of computational complexity for conservation analysis have been reviewed by Sauro and Ingalls (2004). These include computing the Null space of N, Gauss–Jordan method, reduction to row echelon form and singular value decomposition (SVD). However, these methods are prone to errors when handling large networks. Sauro and Ingalls (2004) propose the use of SVD for such systems. However, even the SVD is inaccurate for large systems such as those used as examples later in this paper (see Fig. 1 for a metabolic network of E. coli). The reason for this is that the condition number of the system, defined as the ratio of the largest to the smallest singular value, is very large for nearly singular matrices. That implies that the NR matrix would not be of full rank, as one or more dependent rows has not been eliminated. This results in one of the solution vectors being closely aligned in the direction of the independent vectors, leading to errors when solutions are computed.
|
The objective of our paper is therefore to build a robust means of performing conservation analysis for large biochemical networks. We propose that a numerical tool with such capability can be built using the Householder QR method (Householder, 1958; Hansen, 1992). The QR method factorizes a given matrix A into an orthogonal matrix Q and an upper triangular matrix R. The potential of this method to retain high numerical stability makes it an ideal candidate for conservation analysis. The applicability of the Householder QR method to biochemical networks has in fact been suggested earlier (Holstein and Greenshaw, 1994) but was never explored or developed. While the QR factorization can be obtained using a number of algorithms, including Givens Rotations, Gram–Schmidt method as well as the Householder method, the Householder method is more numerically stable than the other approaches since it uses orthogonal similarity transforms (Householder and Bauer, 1959). For a general m x n matrix, Gaussian elimination requires about 2n3/3 floating point operations (FLOPS). The Householder QR method on the other hand requires about 2n2(m – n/3) FLOPS. This is also better than the QR factorization based on Givens rotations which requires nearly twice as many FLOPS due to additional computation of square roots. In practice, the Givens rotations method is useful only when a relatively few off diagonal elements need to be zeroed, as opposed to the whole lower triangular part of the matrix in Householder method.
|
4 QR DECOMPOSITION |
|---|
A general non-symmetric m x n matrix A can be written as the product of two matrices Q and R as
 | (9) |
The orthogonal matrix Q and the upper trapezoidal matrix R contain information regarding the dependencies in the rows of A. This is of interest to the factorization of the stoichiometry matrix that defines a biochemical network. The key to conservation analysis is the identification of dependent and independent species. Since the rows in the stoichiometry matrix correspond to species, it can be noted that changing the order of the species is equivalent to changing the order of the rows. This indicates that the rows can be permuted until those corresponding to the independent species occupy the first m0 rows yielding NR, while those corresponding to the dependent species will comprise the bottom part of the matrix, yielding N0. It should be noted that while a permutation of the rows of A results in a similarly row-permuted Q, the sign of the values of R in the permuted row is changed, which can be traced back to the way the QR algorithm is implemented. This indicates that if the Q and R matrices are partitioned into m0 and (m – m0) blocks such that
 | (10) |
 | (11) |
From the relationship N0 = L0NR and Equation (11), it follows that L0, Q11 and Q21 are related by 
. The matrix L0 will exist if Q11 is non-singular and hence invertible. A similar observation was made for the case of Gaussian Elimination by Holstein and Greensfaw (1994), who also briefly discussed the merits of the Householder QR factorization for large systems. In the case of the Householder QR however, we found that obtaining an invertible Q11 starting with a general stoichiometric matrix is a slow process. This is due to the fact that after each row permutation, the singular values of Q11 have to be obtained until there are no more zero singular values. The computational cost of repetitive calls to a SVD algorithm renders this approach impractical. In the following, we outline a new approach based on Householder QR for conservation analysis that is robust and that does not face the problem of repetitive singular value computation.
|
5 HOUSEHOLDER QR FACTORIZATION |
|---|
The Householder method was first described by Householder (1958) as an algorithm to triangularize a general non-symmetric matrix. A geometric interpretation of this method entails that we see the network reactions as vectors in the space of species involved. The Householder method is a series of projections of the reactions onto a coordinate system constructed by sequential Householder reflections. These reflections have the effect of creating an upper trazoidal matrix that can then be further reduced to obtain the conservation relations. Hansen (1992) describes this geometrical interpretation of the Householder method in more detail.
A Householder reflection can be defined as a transformation that takes any vector and reflects it about a plane. This plane can be constructed in a manner such that the reflected vector has certain desired properties. This has a simple geometrical interpretation—the method maps the original matrix, each column of which describes a vector—to another matrix where the elements below the diagonal in the column of interest have been made zero. This is done by replacing the column with a vector that reflects it in a particular plane.
The matrix, which in our case, is the stoichiometric matrix for the biochemical network consists of n columns of length m each. These can be treated as n vectors in an m-dimensional space. The Householder algorithm begins by reflecting the first vector onto the first axis of the m-dimensional coordinate system. This has the effect of setting all the elements of the first column of the stoichiometric matrix below the first element to zero. The reflection plane is perpendicular to a given vector v, which is used to construct the Householder matrix H. This matrix is given by H = I – 2vvT, where v is a normalized vector (||v||2 = vTv = 1). The Householder matrices are symmetric (HT = H) and orthogonal (H–1 = HT). Therefore H is a non-singular matrix which is its own inverse (H–1 = H). Each householder matrix eliminates the zeros below the diagonal for each column of the original matrix. Let us denote the first Householder matrix as H1. This matrix can eliminate the elements below the diagonal in the first column of the matrix.
It can therefore be shown that for an m x n stoichiometric matrix, n – 1 such Householder matrices are required to create an R matrix with zeros below the main diagonal. It follows that the product of the Householder matrices is the orthogonal matrix Q in Equation (9) is given as Q = H1H2 Hn–2Hn–1, where Hi is the Householder matrix that can annihilate elements below the diagonal in the i-th column of the stoichiometric matrix.
|
6 CONSERVATION MATRIX EVALUATION |
|---|
The QR decomposition using Householder reflections is available as a part of the LAPACK library (Anderson et al., 1995) which is in the public domain. We used the CLAPACK version to build a module that would connect to the Systems Biology Workbench (SBW, Sauro et al., 2003). This allowed us to develop the conservation analysis module that first builds the Stoichiometric matrix N in full form for a chosen model. The key to obtaining the conservation relations in the model lies in using NT to carry out a QR decomposition. This is because the permutation matrix P contains the new order of the columns of NT that is needed to generate Q and R in a numerically stable manner. The new order of the columns in the transpose is equivalent to the new order of the rows in the original stoichiometric matrix. This can be expressed as NTPt = QtRt where the subscript t indicates that the matrices Qt, Rt and Pt correspond to those obtained from the QR decomposition of the transpose of the original stoichiometric matrix. Since Qt is an orthogonal matrix, multiplying both sides by
leads to the following relation:
 | (12) |
The above form has a matrix Rt on the right-hand side that is upper trapezoidal. Furthermore, if the stoichiometry matrix has rows that are not linearly independent, then Rt can be partitioned into four submatrices that are central to obtaining the conservation laws. The form of Equation (12) can be made simpler by scaling each row of Rt to have a unity on the main diagonal and then performing a Gauss–Jordan elimination to zero out the above diagonal elements. This can be restricted only to the non-zero rows of Rt, yielding
 | (13) |
can then be built from the relation = [–L0 I]. The Link matrix can similarly be obtained using the relation in Equation (3).
The order of the species that satisfies the conservation relation S = T can be obtained using the m x m column permutation matrix Pt. As observed earlier, Pt arises due to rearrangement of the columns of NT, or equivalently rearrangement of the rows of N. Hence the new order of the species can be built by noting the indices of rows of N that have been interchanged. The matrices NR and N0 can then be obtained from the reordered stoichiometry matrix.
The QR decomposition for the stoichiometric matrix results in a very robust numerical scheme to extract the conserved cycles in the biochemical network. This can be seen from the results for some large network models in Table 1. These models were obtained from the website for in silico organisms at the Systems Biology Research Group at University of San Diego, California (see http://gcrg.ucsd.edu/). These models were built from genomic studies and are available as Microsoft Excel spreadsheets that we have converted to SBML (Hucka et al., 2003) format. We compared our implementation of QR and LU with three of the most well-known and established tools, Jarnac (Sauro, 2000), PySCeS (Olivier et al., 2005) and COPASI (http://www.copasi.org).
|
Computation of conservation relations requires only one call to the Householder QR factorization routine, and the permutation matrix returned by the routine can be used to reorder the original species into independent and dependent species lists.
Another point of interest is that the R matrix in the QR decomposition of a general m x n matrix A can be related to A by constructing the matrix product ATA, which is a square matrix. We can then proceed to factorize this square matrix using Cholesky decomposition into the product of an upper triangular matrix and its transpose. It can then be shown that this upper triangular matrix is in fact R. However this approach is not preferred as dependencies in the rows of A are reflected in the ill-conditionality of the product ATA, making it difficult to construct R. It is hence desirable to compute R using QR decomposition instead of using ATA.
It is noted in brief that the procedure described in this section to obtain L0 can also be used to obtain the Kernel matrix K0 and the Null Space (Olivier, 2005; Hofmeyr, 2000). However, instead of NT we should use N to compute these matrices. The K matrix is representative of the number of independent fluxes (Reder, 1988; Heinrich and Schuster, 1996; Klamt and Stelling, 2003).
|
7 AN EXAMPLE NETWORK |
|---|
The Householder QR method described in this paper is now illustrated by means of an example to compute the QR factorization of the network as shown in Figure 2. The original order of species for this network is [ES, E, S1, S2]. The reactions [v1, v2, v3] are
 | (14) |
 | (15) |
|
We now show how the conservation relations for this network can be extracted using the Householder QR approach. The factorization for NT is given by NTPt = QtRt, where Pt and Qt are
 |
 | (16) |
It should be kept in mind that Qt is built out of a product of Householder matrices, each of which results in zeroing elements below the main diagonal of Rt. It can be noted that the rank of the stoichiometric matrix m0 is the number of rows with non-zero elements, in this case, m0 = 2. All the rows below the first m0 non-zero rows will contain only zeros and reflect the dependencies in the network. The objective of our algorithm is to generate L0 matrix, which is done by row operations on R matrix to first convert all diagonal elements to unity by dividing each row with a non-zero diagonal element by its value yielding a modified matrix R*. Performing a Gauss–Jordan reduction to eliminate non-zero values above the diagonal of R* gives 
as a submatrix of R*. This form can be matched to that described by Equation (13), where M corresponds to the transpose of L0. The L0 matrix is then obtained as the transpose of M
 |
The new order of species can be obtained by permuting the original species order with the permutation matrix P, and noting that the first m0 species will comprise the independent species while the remaining ones will be dependent species. In our example, this new order is [ES, S1, E, S2], with ES and S1 being independent species and E and S2 being dependent species. The conserved cycles can now be obtained as a product of the conservation matrix with the reordered species list as
 | (17) |
|
8 COMPARISON TO LU DECOMPOSITION |
|---|
The LU decomposition is a very well-known and standard numerical scheme employed to separate a matrix A into a lower triangular matrix L and an upper triangular matrix U. The LU decomposition is in fact used by PySCeS (Olivier et al., 2005) and COPASI (http://www.copasi.org) to factorize the stoichiometry matrix for structural and conservation analysis. LAPACK (Anderson et al., 1995) contains methods for obtaining the LU decomposition using both partial or full pivoting. Availability of these methods makes it feasible to apply the LU decomposition to conservation analysis in a straightforward manner, particularly for the method that uses partial pivoting. However, the method that employs full pivoting can be used only with square matrices. To overcome this restriction, we modify the stoichiometric matrix that is used as an input to the full pivoting method by adding additional zero rows or columns to convert it into a square matrix. We implemented both these methods using the CLAPACK library to compare the robustness and performance of the Householder QR factorization method. In the following, we briefly outline the implementation of the LU decomposition using partial as well as full pivoting.
The LU decomposition using partial pivoting of a general m x n matrix A is given by A = PLU, where L is an m x min(m, n) lower triangular matrix, U is a min(m, n) x n upper triangular matrix and P is an m x m matrix of row interchanges. While the partial pivoting method generates the L and U matrices, linear dependencies in the rows of A matrix are reflected by some of the diagonal elements of U being zero. The primary shortcoming of partial pivoting is that when LAPACK encounters a singular value on the diagonal of U, the algorithm is stopped and has to be restarted after the zero is exchanged with a non-zero value. An efficient search locates the positions of the zero diagonal elements and rearranges the columns such that all the non-zero diagonal elements are moved to the top of the U matrix. It should be noted that this reordering has to be applied to columns of the transpose of the stoichiometric matrix as well as the order of the species. Once this has been done, the LU decomposition method has to be called again to complete the process that was interrupted by the existence of a zero on the diagonal element. When applied to the transpose of the stoichiometric matrix, the equivalent LU decomposition is given by NT = PtLtUt, where as before, the subscript t is indicative of matrices tracing their origin to the transpose of the stoichiometric matrix. If we denote by (NT)* the stoichiometric matrix whose columns have been rearranged as described above, the second LU decomposition can be written as 
where 
is the new lower triangular matrix, 
is the new upper triangular matrix and 
is the new row reordering matrix. The matrix should now be reducible by Gauss–Jordan elimination of the elements above the diagonal to the form of Equation (13). The L0 matrix can then be extracted as a submatrix of the Gauss–Jordan reduced matrix. It can be seen from Table 1 that the implementation of LU decomposition using partial pivoting generates nearly as many conservation laws as the QR factorization and is faster than the Householder QR, as can be seen from Table 2. However, for large network models (iJR904 and iND750), these conservation laws fail the validity tests described in the following section indicating the lack of reliability in using LU decomposition using partial pivoting on such large systems.
|
In the case of LU decomposition with full pivoting, an additional permutation matrix Q containing column exchanges is used to increase the numerical stability of the method. The resulting decomposition can be written as A = PLUQ. The advantage of using the full pivoting lies in the fact that the computation need be done only once. Further, the structure of the decomposition renders the method to a form that is similar to the Householder method, in that the column permutation matrix Q in LU decomposition with full pivoting plays the same role as the transpose of the permutation matrix P in the Householder method. We also infer that this gives the Householder method a marginal advantage in computational terms, since only one permutation matrix has to be computed. The remaining steps that are needed to compute the L0 matrix from the U matrix are similar to that followed for reducing the Rt matrix in Equation (12) to the Gauss–Jordan reduced form in Equation (13). The L0 matrix can then be obtained as L0 = MT and the conservation laws constructed using the relation
= [–L0I]. We can indeed show that the conservation laws obtained by the Householder decomposition and the LU decomposition with full pivoting are equivalent, in that both pass the validation tests detailed in Section 9. Consequently, we now have two independent methods to verify the conservation laws for large biochemical networks, with the Householder QR method being slightly faster than the full pivoted LU decomposition method. |
9 TESTING CONSERVATION LAWS |
|---|
The conservation laws generated by the algorithm must be tested to ensure that they are valid. Our algorithm therefore includes five tests to validate the conservation laws. The first test checks if the conserved sums are indeed constant values. This is done by observing that the conserved sums are given by T =
S. Their rate of change, dT/dt must be zero if T is a constant vector. This is equivalent to
 | (18) |
The second test checks if the rank obtained by QR factorization of the stoichiometry matrix N is the same as that obtained by counting the non-zero singular values obtained from an SVD of N. If the rank computed by the two methods is different, it is evident that the number of conserved cycles has been miscalculated.
The third test checks the rank of NR matrix. Since the NR matrix corresponds to the contribution of all the independent species in the network, it must have a full rank. Therefore, one can compare the rank of NR using QR factorization to that obtained by SVD. If these ranks are full and match, then the conservation laws are correct.
The fourth test involves the computation of the eigenvalues of a submatrix obtained from the QR factorization of N. We note that the orthogonal matrix Q obtained after factorizing the reordered stoichiometric matrix can be partitioned into four submatrices as described in Equation (10). If the stoichiometric matrix N has been reordered such that the first m0 rows correspond to independent species, and the remaining (m – m0) rows correspond to dependent species, it can be shown that Q11 is a non-singular and invertible matrix. We can therefore compute its rank by counting the number of non-zero eigenvalues of Q11. If Q11 is of full rank, i.e. if m0 = rank(Q11), the conservation laws must be correct.
The fifth test evaluates the L0 matrix as the product of the submatrices of Q. As observed earlier, . The matrix obtained by multiplying Q21 and the inverse of Q11 are compared with the L0 matrix obtained directly by the Householder method. The test can then be deemed successful if both matrices match.
The results for these validity tests depend on the size of the model and the method employed to generate the conservation laws. The numerical errors inherent in the algorithm are highest in the case of LU decomposition with partial pivoting, leading to the generation of conservation laws that do not pass the validity tests. This can be seen in Table 1. On the other hand, the Householder QR method developed in this paper and LU decomposition with full pivoting (used by PySCeS) generate the same number of conservation laws, both of which pass the validity tests.
|
10 SOFTWARE IMPLEMENTATION |
|---|
The software developed for the purpose of demonstrating the capabilities of the Householder QR method is a part of the latest Systems Biology Workbench (Sauro et al., 2003), release 2.5.0 (http://www.sys-bio.org). The conservation analysis software is written in C++, with numerical analysis capability provided by CLAPACK (see http://www.netlib.org/clapack/). An SBW interface to essential routines such as matrix inversion, QR decomposition, LU Decomposition and Eigenvalue computation is available via SBW_CLAPACK module, which is also part of SBW release 2.5.0. A separate module (Structural Analysis Module) provides an interface to higher level routines. The software also includes a graphical interface tool for Structural Analysis developed in C# for the Windows platform, that should in principle run on other platforms as well. This tool allows users to load multiple models to compute their conservation laws and relevant matrices and to store and view them as required. It can be called from other biochemical network modeling tools such as JDesigner (Sauro, 2001, http://sbw.kgi.edu/software/jdesigner/htm; Sauro et al., 2003), and CellDesigner (Funahashi and Kitano, 2003), or any tool that can interface with the Systems Biology Workbench. More detailed information regarding this tool can be obtained from the author’s website at http://public.kgi.edu/~rrao. We are also in the process of developing a new high performance simulator using these tools.
|
11 APPLICATION TO SBML TRANSLATION |
|---|
Many biochemical models are now available in SBML (Hucka et al., 2003) format, which makes it feasible to share and disseminate information regarding new developments (http://www.biomodels.net). However, from a simulation point of view, SBML has to be parsed to obtain information about the species, reactions and rate-laws in order to create a model that can be simulated. Currently existing software tools make it feasible to translate SBML into MATLAB (http://www.mathworks.com/) code that can be run in the MATLAB simulation environment. For small models, this is quite convenient as the time evolution of various species can be graphed easily. In the case of large models however, MATLAB turns out to be very slow as the number of species and reactions increase many fold.
Faced with the slow simulation in MATLAB, we have developed an alternative to this by developing tools that can translate SBML into higher languages such as C, C# and Java, where such simulation can be carried out very quickly for large systems. For instance, a 20-fold speedup in simulation time over MATLAB was achieved using the SBML to C Translator. The translators also provide methods that allow computation of the reduced Jacobian of the system that is very useful for carrying out other analyzes.
The conservation analysis approach presented in this paper has now been integrated into all our SBML translators. This allows the users to translate models in SBML to a more simulation friendly language such as C or C# or Java. The SBML translators call the Structural Analysis methods to obtain a list of the independent and dependent species. The dependent species are then eliminated from the differential equations that have to be simulated using the conservation laws that are generated by the Structural Analysis tool. This also allows for computation of the reduced Jacobian that is non-singular. The users can therefore easily convert existing model information contained in SBML format to other languages. This has been illustrated in Figure 3.
|
Another example of the application of the Structural Analysis module is that of Oscill8 (http://sourceforge.net/projects/oscill8), a Bifurcation Analysis package that utilizes AUTO (Doedel et al., 1991). This package interfaces with Systems Biology Workbench to provide a seamless interface that allows users to study the dynamics of their models. This is done by translating the biochemical model described in SBML to XPP format, which defines the reduced system in terms of ordinary differential equations and their associated initial conditions. This XPP translated model is then used by Oscill8 to carry out bifurcation studies.
|
12 RESULTS AND DISCUSSION |
|---|
A Householder method-based approach for conservation Analysis has been presented in this paper. We have shown that it is a robust method capable of finding the correct conserved cycles for some very large biochemical models. This method uses the CLAPACK SBW library for dense matrices, and a sparse matrix-based solution is under development for very large systems.
The Householder-based method along with the LU decomposition using full pivoting appear to be the most promising of the currently known methods for conservation analysis of large biochemical networks. In addition to extracting the conservation laws, the software developed for this purpose also provides five tests to verify if the conservation laws satisfy required conditions.
Biochemical models involve reactions between various species, and often the action of many species is restricted to just a few reactions. Similarly, some reactions involve very few reactants and products. As a result, the stoichiometric matrix tends to be sparsely populated, with the sparseness increasing with the size of the network. Most analysis methods tend to treat the stoichiometric matrix as a dense matrix where the computations are carried out on all elements of the matrix, even if they are zero. This makes the algorithm slow and at the same time occupies computer memory. A sparse matrix scheme for representing stoichiometric matrices, along with numerical routines capable of performing analysis on sparse matrices will be beneficial in increasing the computational speed and efficiency in future.
Typical sparse matrix representation methods keep a list of the non-zero values, along with their row and column index locations. This is a compact representation that saves memory, but is not efficient from a computational perspective. For most matrix operations that involve rowwise or columnwise operations, a lot of time would be expended on locating the next non-zero row or column element. We are therefore implementing a sparse matrix scheme that not only keeps a list of non-zero values and their column indices, but also assigns pointers to the previous and next non-zero elements, and the previous and next column non-zero elements. This is called a doubly linked data structure and allows for faster element search. We believe this approach will lead to the development of sparse matrix-based analysis tools that will make it possible to study the dynamics of very large biochemical networks.
In this paper, we have not addressed the issue of computing conservation laws where the number of negative terms is minimized (Sauro, 1993). However, we refer the reader to the paper by Kholodenko et al. (1995) which briefly discusses ordering rules to minimize the number of negative terms.
The numerical method developed herein makes it feasible to attempt the more complex problem of simulation of biochemical systems with varied time scales. Typically, these large systems are stiff due to the occurrence of several fast processes which require small step sizes for integrating the solution. Park (1974) first addressed the issue of developing simulators that took into account the temporal hierarchical structure of metabolic networks by splitting the system into slow and fast variables (Heinrich and Schuster, 1996). Such a split assumes that the system moves quickly onto a subset of the allowed phase space, in which its subsequent evolution is slow. Hence, if such a manifold is found, then the task of integrating the systems of differential equations is much simpler. Recently, efficient methods have been developed, which use these intrinsic low-dimensional manifolds (ILDM) to solve large chemical systems, for example in the field of combustion (Mass and Pope, 1992; Nafe and Mass, 2002). In another publication, we will describe the implementation of these methods in conjunction with the conservation analysis approach that we have elaborated upon in this paper.
|
Acknowledgments |
|---|
R.V. would like to acknowledge Brett Olivier for assistance with converting large SBML models to PySCeS format, Emery Conrad for testing the conservation Analysis tool with Oscill8, Anastasia Deckard for generating displays of large biochemical network models and Frank Bergmann for assistance with software development. We would also like to thank the reviewers of this paper for bringing to our attention the fact that PySCeS employs LU decomposition using full pivoting. R.R.V. and V.C. are grateful to the generous support from the US DOE GTL Program and the DARPA/IPTO BioCOMP program, contract number MIPR 03-M296-01 respectively.
Conflict of Interest: none declared.
|
FOOTNOTES |
|---|
Associate Editor: Martin Bishop
1It should be noted that the majority of simulators for systems biology do not attempt a conservation analysis. 
Received on September 23, 2005; revised on November 22, 2005; accepted on November 23, 2005
|
REFERENCES |
|---|
Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A., Sorensen, D. LAPACK User's Guide, (1995) , Philadelphia SIAM.
Aris, R. (1965) Prolegomena to the rational analysis of systems of chemical reactions. Arch. Rational Mech. Anal, . 19, 81–99.
Chickarmane, V., et al. (2005) Bifurcation discovery tool. Bioinformatics, 21, 3688–3690
Clarke, B.L. (1980) Stability of complex reaction networks. In Prigogine, I. and Rice, S.A. (Eds.). Advances in Chemical Physics, , NY Wiley.
Cornish-Bowden, A. and Hofmeyr, J.S. (2002) The role of stoichiometric analysis in studies of metabolism: an example. J. Theor. Biol, . 216, 179–191[CrossRef][Web of Science][Medline].
Cox, A.J. and Higham, N.J. (1997) Stability of Householder QR factorization for weighted least squares problems. Numerical Analysis, , UK TR 301, Department of Mathematics, University of Manchester.
Doedel, E.J., et al. (1991) Numerical analysis and control of bifurcation problems. Part I. Int. J. Bifurcat. Chaos, 1, 493–520[CrossRef].
Duarte, N.C., et al. (2004) Reconstruction and validation of Saccharomyces cerevisiae iND750, a fully compartmentalized genome-scale metabolic model. Genome Res, . 14, 1298–1309
Feinberg, M. (1989) Necessary and sufficient conditions for detailed balancing in mass action systems of arbitary complexity. Chem. Eng. Sci, . 44, 1819–1827[CrossRef].
Funahashi, A. and Kitano, H. (2003) CellDesigner: a process diagram editor for gene-regulatory and biochemical networks. BioSilico, 1, 159–162[CrossRef].
Golub, G.H. (1965) Numerical methods for solving linear least-squares problems. Numer. Math, . 7, 206–216[CrossRef].
Golub, G.H. and Van Loan, C.F. Matrix Computations, (1996) 3rd edn , Baltimore, MD Johns Hopkins University Press.
Hansen, P.B. (1992) Householder reduction of linear equations. ACM Computing Surveys, 24, 185–194.
Heinrich, R. and Schuster, S. The Regulation of Cellular Systems, (1996) Chapman & Hall.
Higham, N.J. (2000) QR factorization with complete pivoting and accurate computation of the SVD. Linear Algebra Appl, . 209, 153–174.
Hofmeyr, J.-H.S. (1986) Steady-state modelling of metabolic pathways: a guide for the prospective simulator. Comput. Appl. Biosci, . 2, 5–11
Hofmeyr, J.-H.S. (2000) Metabolic control analysis in a nutshell. Proceedings of the International Conference on Systems BiologyPasadena, California , pp. 291–300.
Holstein, H. and Greenshaw, C.P. (1994) A numerical treatment of metabolic control methods. In Westerhoff, H.V. (Ed.). Biothermokinetics, Intercept Ltd., pp. 293–299.
Horn, F. and Jackson, R. (1972) General mass action kinetics. Arch. Rational Mech. Anal, . 81–116.
Householder, A.S. (1958) Unitary triangularization of a nonsymmetric matrix. J. ACM, 5, 339–342[CrossRef].
Householder, A.S. and Bauer, F.L. (1959) On certain methods for exapanding the characteristic polynomial. Numerische Math, . 1, 29–37.
Hucka, M., et al. (2003) The Systems Biology Markup Language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics, 19, 524–531
Kholodenko, B.N., et al. (1995) Control in channelled pathways. A matrix method calculating the enzyme control coefficients. Biophys. Chem, . 53, 247–258[Medline].
Klamt, S. and Stelling, J. (2003) Stoichiometric analysis of metabolic networks. proceedings of the Tutorial at the 4th International Conference on Systems BiologyHeidelberg, Germany.
Koshland, D.E., et al. (1982) Amplification and Adaptation in Regulatory and Sensory Systems. Science, 217, , pp. 220–225
Maas, U.A. and Pope, S.B. (1992) Simplifying chemical kinetics: intrinsic low-dimensional manifolds in composition space. Combustion Flame, 88, 239–264.
Mendes, P. (1993) GEPASI: a software package for modelling the dynamics, steady states and control of biochemical and other systems. Comp. Appl. Biosci, . 9, 563–571[Medline].
Nafe, J. and Mass, U. (2002) A general algorithm for improving ILDMs. Combustion Theor. Modeling, 6, 697–709[CrossRef].
Olivier, B.G. (2005) Simulation and database software for computational systems biology: PySCeS and JWS Online. PhD Thesis, Stellenbosch University.
Olivier, B.G., et al. (2005) Modelling cellular systems with PySCeS. Bioinformatics, 21, 560–561
Park, D.J.M. (1974) The hierarchical structure of metabolic networks and the construction of efficient metabolic simulators. J. Theor. Biol, . 46, 31–74[CrossRef][Web of Science][Medline].
Reder, C. (1988) Metabolic control theory: a structural approach. J. Theor. Biol, . 135, 175–201[CrossRef][Web of Science][Medline].
Reed, J.L., et al. (2003) An expanded genome-scale model of Escherichia coli K-12 (iJR904 GSM/GPR). Genome Biol, . 4, R54.1–R54.12.
Reich, J.G. and Selkov, E.E. (1981) Energy metabolism of the cell. , London Academic Press.
Sauro, H.M. (1993) SCAMP: a general-purpose simulator and metabolic control analysis program. Comput. Appl. Biosci, . 9, 441–450
Sauro, H.M. (2000) Jarnac: a system for interactive metabolic analysis. In Hofmeyr, J.-H.S., Rohwer, J.M., Snoep, J.L. (Eds.). Animating the Cellular Map. Intl. BioThermoKinetics Meeting, Stellenbosch Univeristy Press. Chapter 33, pp. 221–228.
Sauro, H.M. (2001) JDesigner: a simple biochemical network designer.
Sauro, H.M., et al. (2003) Next generation simulation tools: the Systems Biology Workbench and BioSPICE integration. OMICS Winter, 7, 355–372.
Sauro, H.M. and Ingalls, B. (2004) Conservation analysis in biochemical networks: computational issues for software writers. Biophys. Chem, . 109, 1–15[CrossRef][Web of Science][Medline].
Schilling, C.H., et al. (2002) Genome-scale metabolic model of Helicobacter pylori 26695. J. Bacteriol, . 184, 4582–4593
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  S. Hoops, S. Sahle, R. Gauges, C. Lee, J. Pahle, N. Simus, M. Singhal, L. Xu, P. Mendes, and U. Kummer COPASI--a COmplex PAthway SImulator Bioinformatics, December 15, 2006; 22(24): 3067 - 3074. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 F. J Doyle III and J. Stelling Systems interface biology J R Soc Interface, October 22, 2006; 3(10): 603 - 616. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||